text:
using Base: has_offset_axes

## (GE) general matrices, solvers with factorization, solver and inverse
# for (gels, gesv, getrs, getri, elty) in
#     ((:dgels,:dgesv,:dgetrs,:dgetri,:Float64),
#      (:sgels,:sgesv,:sgetrs,:sgetri,:Float32),
#      (:zgels,:zgesv,:zgetrs,:zgetri,:ComplexF64),
#      (:cgels,:cgesv,:cgetrs,:cgetri,:ComplexF32))
const function_list = ("gels", "gesv", "getrs", "getri", "getrf")
for type in magmaTypeList
    # create the symbols for element types
    elty = Symbol(type)
    # generate the symbol variables for our wrappers
    for func_name in function_list
        @eval $(Symbol(func_name)) = (Symbol(magmaTypeDict[$type], $func_name))
    end

    @eval begin
        #      SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,INFO)
        # *     .. Scalar Arguments ..
        #       CHARACTER          TRANS
        #       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
        function magma_gels!(trans::magma_trans_t, A::CuArray{$elty}, B::CuArray{$elty})
            m, n  = size(A)
            info  = Ref{Int}()
            work  = Vector{$elty}(undef, 1)
            lwork = Int(-1)

            for i = 1:2  # first call returns lwork as work[1]
                func = eval(@magmafunc_gpu($gels))
                func(
                    MagmaNoTrans, m, n, size(B,2),
                    A, max(1,stride(A,2)), B, max(1,stride(B,2)),
                    work, lwork, info
                )
                # ccall((@magmafunc_gpu($gels), libmagma), Int,
                #       (Int, Int, Int, Int,
                #        PtrOrCuPtr{$elty}, Int, PtrOrCuPtr{$elty}, Int,
                #        Ptr{$elty}, Int, Ptr{Int}),
                #       MagmaNoTrans, m, n, size(B,2),
                #       A, max(1,stride(A,2)), B, max(1,stride(B,2)),
                #       work, lwork, info)

                if i == 1
                    lwork = ceil(Int, real(work[1]))
                    resize!(work, lwork)
                end
            end

            k   = min(m, n)
            F   = m < n ? tril(A[1:k, 1:k]) : triu(A[1:k, 1:k])
            ssr = Vector{$elty}(undef, size(B, 2))
            for i = 1:size(B,2)
                x = zero($elty)
                for j = k+1:size(B,1)
                    x += abs2(B[j,i])
                end
                ssr[i] = x
            end
            F, subsetrows(B, B, k), ssr
        end
        function magma_gels!(trans::magma_trans_t, A::Array{$elty}, B::Array{$elty})
            m, n  = size(A)
            info  = Ref{Int}()
            work  = Vector{$elty}(undef, 1)
            lwork = Int(-1)

            for i = 1:2  # first call returns lwork as work[1]
                func = eval(@magmafunc($gels))
                func(
                    MagmaNoTrans, m, n, size(B,2),
                    A, max(1,stride(A,2)), B, max(1,stride(B,2)),
                    work, lwork, info
                )
                # ccall((@magmafunc($gels), libmagma), Int,
                #       (Int, Int, Int, Int,
                #        PtrOrCuPtr{$elty}, Int, PtrOrCuPtr{$elty}, Int,
                #        Ptr{$elty}, Int, Ptr{Int}),
                #       MagmaNoTrans, m, n, size(B,2),
                #       A, max(1,stride(A,2)), B, max(1,stride(B,2)),
                #       work, lwork, info)

                if i == 1
                    lwork = ceil(Int, real(work[1]))
                    resize!(work, lwork)
                end
            end

            k   = min(m, n)
            F   = m < n ? tril(A[1:k, 1:k]) : triu(A[1:k, 1:k])
            ssr = Vector{$elty}(undef, size(B, 2))
            for i = 1:size(B,2)
                x = zero($elty)
                for j = k+1:size(B,1)
                    x += abs2(B[j,i])
                end
                ssr[i] = x
            end
            F, subsetrows(B, B, k), ssr
        end

        #     Parameters
        # [in]	     n	    INTEGER The order of the matrix A. N >= 0.
        # [in]	     nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
        # [in,out]	 A	    COMPLEX_16 array, dimension (LDA,N). On entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A = P*L*U; the unit diagonal elements of L are not stored.
        # [in]	     lda	INTEGER The leading dimension of the array A. LDA >= max(1,N).
        # [out]	     ipiv	INTEGER array, dimension (min(M,N)) The pivot indices; for 1 <= i <= min(M,N), row i of the matrix was interchanged with row IPIV(i).
        # [in,out]	 B	    COMPLEX_16 array, dimension (LDB,NRHS) On entry, the right hand side matrix B. On exit, the solution matrix X.
        # [in]	     ldb	INTEGER The leading dimension of the array B. LDB >= max(1,N).
        # [out]	     info	INTEGER
        #
        #               = 0: successful exit
        #               < 0: if INFO = -i, the i-th argument had an illegal value
        function magma_gesv!(A::CuArray{$elty}, B::CuArray{$elty})
            # require_one_based_indexing(A, B)
            # chkstride1(A, B)
            n = checksquare(A)
            isGPU = (A isa CuArray && B isa CuArray)
            # if size(B,1) != n
            #     throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n"))
            # end
            lda  = max(1,stride(A,2))
            ldb  = max(1,stride(B,2))
            ipiv = similar(Matrix(A), Int, n)
            info = Ref{Int}()

            func = eval(@magmafunc_gpu($gesv))
            func(
                n, size(B,2),
                A, lda, ipiv,
                B, ldb, info
            )
            # ccall((@magmafunc_gpu($gesv), libmagma), Int,
            #       (Int, Int,
            #        PtrOrCuPtr{$elty}, Int, Ptr{Int},
            #        PtrOrCuPtr{$elty}, Int, Ptr{Int}),
            #        n, size(B,2),
            #        A, lda, ipiv,
            #        B, ldb, info)
            # chkmagmaerror(info[])
            B, A, ipiv
        end

        function magma_gesv!(A::Array{$elty}, B::Array{$elty})
            # require_one_based_indexing(A, B)
            # chkstride1(A, B)
            n = checksquare(A)
            isGPU = (A isa CuArray && B isa CuArray)
            # if size(B,1) != n
            #     throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n"))
            # end
            lda  = max(1,stride(A,2))
            ldb  = max(1,stride(B,2))
            ipiv = similar(Matrix(A), Int, n)
            info = Ref{Int}()

            func = eval(@magmafunc($gesv))
            func(
                n, size(B,2),
                A, lda, ipiv,
                B, ldb, info
            )
            # ccall((@magmafunc($gesv), libmagma), Int,
            #       (Int, Int,
            #        PtrOrCuPtr{$elty}, Int, Ptr{Int},
            #        PtrOrCuPtr{$elty}, Int, Ptr{Int}),
            #        n, size(B,2),
            #        A, lda, ipiv,
            #        B, ldb, info)

            # chkmagmaerror(info[])
            B, A, ipiv
        end

        #     SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
        #*     .. Scalar Arguments ..
        #      CHARACTER          TRANS
        #      INTEGER            INFO, LDA, LDB, N, NRHS
        #     .. Array Arguments ..
        #      INTEGER            IPIV( * )
        #      DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
        function magma_getrs!(trans::magma_trans_t, A::CuMatrix{$elty}, ipiv::Array{Int}, B::CuArray{$elty})
            # require_one_based_indexing(A, ipiv, B)
            # chktrans(trans)
            # chkstride1(A, B, ipiv)
            @assert !has_offset_axes(A, ipiv, B)
            n = checksquare(A)
            if n != size(B, 1)
                throw(DimensionMismatch("B has leading dimension $(size(B,1)), but needs $n"))
            end
            nrhs = size(B, 2)
            lda  = max(1,stride(A,2))
            ldb  = max(1,stride(B,2))
            info = Ref{Int}()

            func = eval(@magmafunc_gpu($getrs))
            func(trans, n, nrhs, A, lda, ipiv, B, ldb, info)
            B
            # println("\nThe info variable is: ", info, "\n")
        end
        function magma_getrs!(trans::magma_trans_t, A::Matrix{$elty}, ipiv::AbstractVector{Int}, B::Array{$elty})
            magma_getrs!(trans, cu(A), ipiv, cu(B))
        end

        #*# Parameters of getri
        #*# [in]	n	INTEGER The order of the matrix A. N >= 0.
        #*# [in,out]	dA	COMPLEX array on the GPU, dimension (LDDA,N) On entry, the factors L and U from the factorization A = P*L*U as computed by CGETRF_GPU. On exit, if INFO = 0, the inverse of the original matrix A.
        #*# [in]	ldda	INTEGER The leading dimension of the array A. LDDA >= max(1,N).
        #*# [in]	ipiv	INTEGER array, dimension (N) The pivot indices from CGETRF; for 1 <= i <= N, row i of the matrix was interchanged with row IPIV(i).
        #*# [out]	dwork	(workspace) COMPLEX array on the GPU, dimension (MAX(1,LWORK))
        #*# [in]	lwork	INTEGER The dimension of the array DWORK. LWORK >= N*NB, where NB is the optimal blocksize returned by magma_get_cgetri_nb(n).
        #*# Unlike LAPACK, this version does not currently support a workspace query, because the workspace is on the GPU.
        #*# [out]	info	INTEGER
        #*#         = 0: successful exit
        #*#         < 0: if INFO = -i, the i-th argument had an illegal value
        #*#         > 0: if INFO = i, U(i,i) is exactly zero; the matrix is singular and its cannot be computed.
        function magma_getri!(A::CuMatrix{$elty}, ipiv::Array{Int})
            # require_one_based_indexing(A, ipiv)
            # chkstride1(A, ipiv)
            n = checksquare(A)
            if n != length(ipiv)
                throw(DimensionMismatch("ipiv has length $(length(ipiv)), but needs $n"))
            end
            lda = max(1,stride(A, 2))

            func_nb = eval(@magmafunc_nb($getri))
            nb    = func_nb(n)
            lwork = ceil(Int, real(n * nb))
            work  = cu(Vector{$elty}(undef, max(1, lwork)))
            info  = Ref{Int}()
            # for i = 1:2  # first call returns lwork as work[1]
                # ccall((@magmafunc_gpu($getri), libmagma), Cvoid, # ! MAGMA has no native cpu interface for getri
                #       (Ref{Int}, PtrOrCuPtr{$elty}, Ref{Int}, PtrOrCuPtr{Int},
                #        PtrOrCuPtr{$elty}, Ref{Int}, PtrOrCuPtr{Int}),
                #       n, A, lda, ipiv, work, lwork, info)
                # # chkmagmaerror(info[])
            func = eval(@magmafunc_gpu($getri))
            func(n, A, lda, ipiv, work, lwork, info)
            #     if i == 1
            #         lwork = ceil(Int, real(work[1]))
            #         resize!(work, lwork)
            #     end
            # end
            A
        end
        function magma_getri!(A::Matrix{$elty}, ipiv::Array{Int})
            magma_getri!(cu(A), ipiv)
        end

        function magma_getrf!(A::AbstractMatrix{$elty})
            @assert !has_offset_axes(A)
            chkstride1(A)
            m, n = size(A)
            lda = max(1, stride(A, 2))
            ipiv = similar(A, Int, min(m, n))
            info = Ref{Int}()
            func = eval(@magmafunc($getrf))
            func(m, n, A, lda, ipiv, info)
            A, ipiv, info[]
        end
    end
end


text:
## ------------------------------------------------------------------
# Meta

get_book(d::RBDoc) = get_parent(d)
set_book!(d::RBDoc, book::RBook) = set_parent!(d, book)

get_title(d::RBDoc) = getproperty!(d, :title, "")
set_title!(d::RBDoc, title::String) = setproperty!(d, :title, strip(title))

get_doi(d::RBDoc) = getproperty!(d, :doi, "")
set_doi!(d::RBDoc, doi::String) = setproperty!(d, :doi, _doi_to_url(strip(doi)))

get_year(d::RBDoc) = getproperty!(d, :year, "")
set_year!(d::RBDoc, year::String) = setproperty!(d, :year, strip(year))

get_bibkey(d::RBDoc) = getproperty!(d, :bibkey, "")
set_bibkey!(d::RBDoc, bibkey::String) = setproperty!(d, :bibkey, strip(bibkey))

get_abstract(d::RBDoc) = getproperty!(d, :abstract, "")
set_abstract!(d::RBDoc, abstract::String) = setproperty!(d, :abstract, strip(abstract))

get_ctime(d::RBDoc) = getproperty!(d, :ctime, "")
set_ctime!(d::RBDoc, ctime::String) = setproperty!(d, :ctime, strip(ctime))

get_author(d::RBDoc) = getproperty!(d, :authors, Set{String}())
set_author!(d::RBDoc, author::String, authors::String...) = (_push_csv!(empty!(get_author(d)), author, authors...); d)

## ------------------------------------------------------------------
# Data
sections(d::RBDoc) = getproperty!(() -> OrderedDict{String, RBSection}(), d, :secs)
hassec(d::RBDoc, label::String) = haskey(sections(d), strip(label))

## ------------------------------------------------------------------
# OrderedDict
Base.length(d::RBDoc) = length(sections(d))
Base.getindex(d::RBDoc, key::String) = getindex(sections(d), key)
Base.get(d::RBDoc, key::String, default) = get(sections(d), key, default)
Base.get!(d::RBDoc, key::String, default) = get!(sections(d), key, default)

## ------------------------------------------------------------------
# Array
Base.getindex(d::RBDoc, idx) = (secs = sections(d); getindex(_values(secs), idx))
Base.lastindex(d::RBDoc) = lastindex(_values(sections(d)))
Base.firstindex(d::RBDoc) = firstindex(_values(sections(d)))
Base.iterate(d::RBDoc) = iterate(_values(sections(d)))
Base.iterate(d::RBDoc, state) = iterate(_values(sections(d)), state)

## ------------------------------------------------------------------
# show
function Base.show(io::IO, d::RBDoc)
    nsecs = length(d)
    println(io, "RBDoc(\"", get_label(d), "\") with ", nsecs, " section(s)")
    print(io, " - path: '", localrelpath(d), "'")

    # meta
    for meta in [:title, :doi]
        str = getproperty(d, meta, "")
        if !isempty(str)
            print(io, "\n - ", meta, ": \"", _preview(io, str), "\"")
        end
    end

    # data
    if nsecs > 0
        print(io, "\nsection(s):")
        _show_data_preview(io, d) do sec
            string("[", get_label(sec), "] ", get_title(sec))
        end
    end
    print(io, "\n", _preview(io, "-"^70))
    return nothing
end

## ------------------------------------------------------------------
# API
function eachobj(doc::RBDoc)
    ch = Channel{RBObject}(1) do ch_
        for sec in doc
            put!(ch_, sec)
            for obj in sec
                put!(ch_, obj)
            end
        end
    end

    return (obj for obj in ch)
end

# This should be the only method to add nre documents
# It updates all sort of stuff
function add_sec!(doc::RBDoc, sec::RBSection)

    seclabel = get_label(sec)

    # Update stuff
    # obj labels
    add_objlabel!(sec)

    # add doc
    secs = sections(doc)
    push!(secs, seclabel => sec)

    return doc
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

struct IntSet <: AbstractSet{Int}
    bits::BitVector
    IntSet() = new(falses(256))
end

"""
    IntSet([itr])

Construct a sorted set of positive `Int`s generated by the given iterable object, or an
empty set. Implemented as a bit string, and therefore designed for dense integer sets. Only
`Int`s greater than 0 can be stored. If the set will be sparse (for example holding a few
very large integers), use [`Set`](@ref) instead.
"""
IntSet(itr) = union!(IntSet(), itr)

eltype(::Type{IntSet}) = Int
similar(s::IntSet) = IntSet()
copy(s1::IntSet) = copy!(IntSet(), s1)
function copy!(dest::IntSet, src::IntSet)
    resize!(dest.bits, length(src.bits))
    copy!(dest.bits, src.bits)
    dest
end
eltype(s::IntSet) = Int
sizehint!(s::IntSet, n::Integer) = (_resize0!(s.bits, max(n, length(s.bits))); s)

# An internal function for setting the inclusion bit for a given integer n >= 0
@inline function _setint!(s::IntSet, idx::Integer, b::Bool)
    if idx > length(s.bits)
        b || return s # setting a bit to zero outside the set's bits is a no-op
        newlen = idx + idx>>1 # This operation may overflow; we want saturation
        _resize0!(s.bits, ifelse(newlen<0, typemax(Int), newlen))
    end
    @inbounds s.bits[idx] = b
    s
end

# An internal function to resize a bitarray and ensure the newly allocated
# elements are zeroed (will become unnecessary if this behavior changes)
@inline function _resize0!(b::BitVector, newlen::Integer)
    len = length(b)
    resize!(b, newlen)
    len < newlen && @inbounds b[len+1:newlen] = false # resize! gives dirty memory
    b
end

# An internal function that takes a pure function `f` and maps across two BitArrays
# allowing the lengths to be different and altering b1 with the result
function _matched_map!(f, b1::BitArray, b2::BitArray)
    l1, l2 = length(b1), length(b2)
    if l1 == l2
        map!(f, b1, b1, b2)
    elseif l1 < l2
        _resize0!(b1, l2)
        map!(f, b1, b1, b2)
    elseif l1 > l2
        if f(false, false) == f(true, false) == false
            # We don't need to worry about the trailing bits — they're all false
            resize!(b1, l2)
            map!(f, b1, b1, b2)
        else
            # We transiently extend b2 — as IntSet internal storage this is unobservable
            _resize0!(b2, l1)
            map!(f, b1, b1, b2)
            resize!(b2, l2)
        end
    end
    b1
end

@noinline _throw_intset_bounds_err() = throw(ArgumentError("elements of IntSet must be between 1 and typemax(Int)"))
@noinline _throw_keyerror(n) = throw(KeyError(n))

@inline function push!(s::IntSet, n::Integer)
    0 < n <= typemax(Int) || _throw_intset_bounds_err()
    _setint!(s, n, true)
end
push!(s::IntSet, ns::Integer...) = (for n in ns; push!(s, n); end; s)

@inline function pop!(s::IntSet)
    pop!(s, last(s))
end
@inline function pop!(s::IntSet, n::Integer)
    n in s ? (_delete!(s, n); n) : _throw_keyerror(n)
end
@inline function pop!(s::IntSet, n::Integer, default)
    n in s ? (_delete!(s, n); n) : default
end
@inline _delete!(s::IntSet, n::Integer) = _setint!(s, n, false)
@inline delete!(s::IntSet, n::Integer) = n < 0 ? s : _delete!(s, n)
shift!(s::IntSet) = pop!(s, first(s))

empty!(s::IntSet) = (fill!(s.bits, false); s)
isempty(s::IntSet) = !any(s.bits)

# Mathematical set functions: union!, intersect!, setdiff!, symdiff!

union(s::IntSet) = copy(s)
union(s1::IntSet, s2::IntSet) = union!(copy(s1), s2)
union(s1::IntSet, ss::IntSet...) = union(s1, union(ss...))
union(s::IntSet, ns) = union!(copy(s), ns)
union!(s::IntSet, ns) = (for n in ns; push!(s, n); end; s)
function union!(s1::IntSet, s2::IntSet)
    _matched_map!(|, s1.bits, s2.bits)
    s1
end

intersect(s1::IntSet) = copy(s1)
intersect(s1::IntSet, ss::IntSet...) = intersect(s1, intersect(ss...))
function intersect(s1::IntSet, ns)
    s = IntSet()
    for n in ns
        n in s1 && push!(s, n)
    end
    s
end
intersect(s1::IntSet, s2::IntSet) =
    (length(s1.bits) >= length(s2.bits) ? intersect!(copy(s1), s2) : intersect!(copy(s2), s1))
"""
    intersect!(s1::IntSet, s2::IntSet)

Intersects sets `s1` and `s2` and overwrites the set `s1` with the result. If needed, `s1`
will be expanded to the size of `s2`.
"""
function intersect!(s1::IntSet, s2::IntSet)
    _matched_map!(&, s1.bits, s2.bits)
    s1
end

setdiff(s::IntSet, ns) = setdiff!(copy(s), ns)
setdiff!(s::IntSet, ns) = (for n in ns; _delete!(s, n); end; s)
function setdiff!(s1::IntSet, s2::IntSet)
    _matched_map!(>, s1.bits, s2.bits)
    s1
end

symdiff(s::IntSet, ns) = symdiff!(copy(s), ns)
"""
    symdiff!(s, itr)

For each element in `itr`, destructively toggle its inclusion in set `s`.
"""
symdiff!(s::IntSet, ns) = (for n in ns; symdiff!(s, n); end; s)
"""
    symdiff!(s, n)

The set `s` is destructively modified to toggle the inclusion of integer `n`.
"""
function symdiff!(s::IntSet, n::Integer)
    0 < n < typemax(Int) || _throw_intset_bounds_err()
    val = !(n in s)
    _setint!(s, n, val)
    s
end
function symdiff!(s1::IntSet, s2::IntSet)
    _matched_map!(xor, s1.bits, s2.bits)
    s1
end

@inline function in(n::Integer, s::IntSet)
    if 1 <= n <= length(s.bits)
        @inbounds b = s.bits[n]
    else
        b = false
    end
    b
end

# Use the next-set index as the state to prevent looking it up again in done
start(s::IntSet) = next(s, 0)[2]
function next(s::IntSet, i)
    nextidx = i == typemax(Int) ? 0 : findnext(s.bits, i+1)
    (i, nextidx)
end
done(s::IntSet, i) = i <= 0

@noinline _throw_intset_notempty_error() = throw(ArgumentError("collection must be non-empty"))

function first(s::IntSet)
    idx = findfirst(s.bits)
    idx == 0 ? _throw_intset_notempty_error() : idx
end

function last(s::IntSet)
    idx = findprev(s.bits, length(s.bits))
    idx == 0 ? _throw_intset_notempty_error() : idx
end

length(s::IntSet) = sum(s.bits)

function show(io::IO, s::IntSet)
    print(io, "IntSet([")
    first = true
    for n in s
        !first && print(io, ", ")
        print(io, n)
        first = false
    end
    print(io, "])")
end

function ==(s1::IntSet, s2::IntSet)
    l1 = length(s1.bits)
    l2 = length(s2.bits)
    # If the lengths are the same, simply punt to bitarray comparison
    l1 == l2 && return s1.bits == s2.bits

    # Swap so s1 is always longer
    if l1 < l2
        s2, s1 = s1, s2
        l2, l1 = l1, l2
    end
    # Iteratively check the chunks of the bitarrays
    c1 = s1.bits.chunks
    c2 = s2.bits.chunks
    @inbounds for i in 1:length(c2)
        c1[i] == c2[i] || return false
    end
    # Ensure remaining chunks are zero
    @inbounds for i in length(c2)+1:length(c1)
        c1[i] == UInt64(0) || return false
    end
    return true
end

issubset(a::IntSet, b::IntSet) = isequal(a, intersect(a,b))
<(a::IntSet, b::IntSet) = (a<=b) && !isequal(a,b)
<=(a::IntSet, b::IntSet) = issubset(a, b)

const hashis_seed = UInt === UInt64 ? 0x88989f1fc7dea67d : 0xc7dea67d
function hash(s::IntSet, h::UInt)
    h ⊻= hashis_seed
    bc = s.bits.chunks
    i = length(bc)
    while i > 0 && bc[i] == UInt64(0)
        # Skip trailing empty bytes to prevent extra space from changing the hash
        i -= 1
    end
    while i > 0
        h = hash(bc[i], h)
        i -= 1
    end
    h
end

minimum(s::IntSet) = first(s)
maximum(s::IntSet) = last(s)
extrema(s::IntSet) = (first(s), last(s))
issorted(s::IntSet) = true

max_stars_repo_name: jdlangs/julia

max_issues_repo_name: jdlangs/julia

max_forks_repo_name: jdlangs/julia


text:
function inv(m,n)
    a = m
    b = n
    u_a = 1
    v_b = 1
    u_b = 0
    v_a = 0

    if NOD(a,b)!=1
        return Nothing
    end

    while b != 0
        k, r = divrem(a,b)
        a, b = b, r
        u, v = u_a, v_a
        u_a, v_a = u_b, u_a
        u_b, v_b = u-k*u_b, v-k*v_b
    end
    while u_a<0
        u_a+=n
    end
    return u_a
end

function NOD(a,b)
    m = a
    n = b
    while n != 0
        m, n = n, m % n
    end
    return m
end

Base.inv(a::ResidueRing{m}) where m = ResidueRing{m}.inv(a,m)


text:
"""
A SYK model solver based on a forward fixed-point iteration method.

 The self-energy of the SYK model is given by,

    Σ(τ) = J² * G(τ) * G(τ) * G(β-τ),

 where Green's function of the SYK model is given by the Dyson equation,

    G(iωₙ) = -1/(iωₙ -μ + Σ(iωₙ))

 We solve the Dyson equation self-consistently by a weighted fixed point iteration,
 with weight `mix` assigned to the new iterate and weight `1-mix` assigned to the previous iterate.

 The self-energy is evaluated in the imaginary time domain,
 and the Dyson equation is solved in the Matsubara frequency domain.

 The SYK Green's function has particle-hole symmetry when μ=0.
 You may enforce such symmetry by setting `symmetry = :ph` when initialize the DLR grids.
 A symmetrized solver tends to be more robust than a unsymmetrized one.
"""

using Lehmann
using Printf

diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b

conformal_tau(τ, β) = π^(1 / 4) / sqrt(2β) * 1 / sqrt(sin(π * τ / β))

function syk_sigma_dlr(d, G_x, J = 1.0; sumrule = nothing, verbose = false)

    tau_k = d.τ # DLR imaginary time nodes
    tau_k_rev = d.β .- tau_k # Reversed imaginary time nodes

    G_x_rev = tau2tau(d, G_x, tau_k_rev, sumrule = sumrule, verbose = verbose) # G at beta - tau_k

    Sigma_x = J .^ 2 .* G_x .^ 2 .* G_x_rev # SYK self-energy in imaginary time

    return Sigma_x
end

function dyson(d, sigma_q, mu)
    if d.symmetry == :ph #symmetrized G
        @assert mu ≈ 0.0 "Only the case μ=0 enjoys the particle-hole symmetry."
        return 1im * imag.(-1 ./ (d.ωn * 1im .- mu .+ sigma_q))
    elseif d.symmetry == :none
        return -1 ./ (d.ωn * 1im .- mu .+ sigma_q)
    else
        error("Not implemented!")
    end
end

function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 5000, G_x = zeros(ComplexF64, length(d)), sumrule = nothing, verbose = true)

    for iter in 1:maxiter

        Sigma_x = syk_sigma_dlr(d, G_x, sumrule = sumrule, verbose = verbose)

        G_q_new = dyson(d, tau2matfreq(d, Sigma_x), mu)

        G_x_new = matfreq2tau(d, G_q_new, sumrule = sumrule, verbose = verbose)

        if verbose
            if iter % (maxiter / 10) == 0
                println("round $iter: change $(diff(G_x_new, G_x))")
            end
        end
        if maximum(abs.(G_x_new .- G_x)) < tol && iter > 10
            break
        end

        G_x = mix * G_x_new + (1 - mix) * G_x # Linear mixing
    end
    return G_x
end

function printG(d, G_x)
    @printf("%15s%40s%40s%40s\n", "τ", "DLR imag", "DLR real", "asymtotically exact")
    for i in 1:d.size
        if d.τ[i] <= d.β / 2
            @printf("%15.8f%40.15f%40.15f%40.15f\n", d.τ[i], imag(G_x[i]), real(G_x[i]), conformal_tau(d.τ[i], d.β))
        end
    end
    println()
end

verbose = false

printstyled("=====    Prepare the expected Green's function of the SYK model     =======\n", color = :yellow)
dsym_correct = DLRGrid(Euv = 5.0, β = 10000.0, isFermi = true, rtol = 1e-14, symmetry = :ph) # Initialize DLR object
G_x_correct = solve_syk_with_fixpoint_iter(dsym_correct, 0.00, mix = 0.1, verbose = false)
printG(dsym_correct, G_x_correct)

printstyled("=====    Test Symmetrized and Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)

@printf("%30s%30s%30s%30s%20s\n", "Euv", "sym_solver", "unsym_solver", "unsym_solver+sum_rule", "good or bad")
for Euv in LinRange(5.0, 10.0, 50)

    rtol = 1e-10
    β = 10000.0
    # printstyled("=====     Symmetrized DLR solver for SYK model     =======\n", color = :yellow)
    mix = 0.01
    dsym = DLRGrid(Euv = Euv, β = β, isFermi = true, rtol = rtol, symmetry = :ph, rebuild = true, verbose = false) # Initialize DLR object
    G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix, sumrule = nothing, verbose = verbose)

    # printstyled("=====     Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)
    mix = 0.01
    dnone = DLRGrid(Euv = Euv, β = β, isFermi = true, rtol = rtol, symmetry = :none, rebuild = true, verbose = false) # Initialize DLR object
    G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, sumrule = nothing, verbose = verbose)

    # printstyled("=====     Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)
    mix = 0.01
    G_x_none_sumrule = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, sumrule = 1.0, verbose = verbose)
    # printG(dnone, G_x_none)

    # printstyled("=====     Unsymmetrized versus Symmetrized DLR solver    =======\n", color = :yellow)
    # @printf("%15s%40s%40s%40s\n", "τ", "sym DLR (interpolated)", "unsym DLR", "difference")
    # G_x_interp = tau2tau(dsym_correct, G_x_correct, dnone.τ)
    # for i in 1:dnone.size
    #     if dnone.τ[i] <= dnone.β / 2
    #         @printf("%15.8f%40.15f%40.15f%40.15f\n", dnone.τ[i], real(G_x_interp[i]), real(G_x_none[i]), abs(real(G_x_interp[i] - G_x_none[i])))
    #     end
    # end

    G_x_interp_ph = tau2tau(dsym_correct, G_x_correct, dsym.τ)
    G_x_interp_none = tau2tau(dsym_correct, G_x_correct, dnone.τ)
    G_x_interp_none_sumrule = tau2tau(dsym_correct, G_x_correct, dnone.τ)
    d_ph = diff(G_x_interp_ph, G_x_ph)
    d_none = diff(G_x_interp_none, G_x_none)
    d_none_sumrule = diff(G_x_interp_none_sumrule, G_x_none_sumrule)
    flag = (d_ph < 100rtol) && (d_none < 100rtol) && (d_none_sumrule < 100rtol) ? "good" : "bad"

    @printf("%30.15f%30.15e%30.15e%30.15e%20s\n", Euv, d_ph, d_none, d_none_sumrule, flag)
    # println("symmetric Euv = $Euv maximumal difference: ", diff(G_x_interp, G_x_ph))
    # println("non symmetric Euv = $Euv maximumal difference: ", diff(G_x_interp, G_x_none))

end


text:
using SpecialFunctions, QuadGK, ForwardDiff

##### Diffusion approximation for 2-layered medium with matched index of refraction #####

"""
    fluence_spatial_freq(z, s, μa, μsp, l)

Compute the spatial frequency components of the fluence for 2-layered model

# Arguments
- `z`: depth (z=0 at the surface, z=l at interface of two layers).
- `s`: spatial frequency component s^2 = sx^2 + sy^2
- `μa`: absorption coefficient of both layers [μa1, μa2]
- `μsp`: reduced scattering coefficient of both layers [μsp1, μsp2]
- `l`: thickness of first layer

# Examples
```jldoctest
julia> fluence_spatial_freq(0., 1,  [0.005,0.1], [10,100.], 1)
0.8647044203305079
```
"""
function fluence_spatial_freq(z, s, μa, μsp, l)

    #get numerical instabillity at large s
    D = @. 1/(3(μsp + μa))
    α = @. sqrt((D*s^2 + μa)/D)
    A = 0.493
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    ## equation 11
    if z < z0

        # have to account for exp blowing up to inf
        if (α[1]*z0 - α[1]*z) > 15
            ϕ = exp(α[1]*(z0 - 2*l)) - exp(α[1]*(-z0 - 2*l - zb)) - exp(α[1]*(z0 - 2*l - 2*zb)) + exp(α[1]*(-z0 - 2*l - 4*zb))
            ϕ = ϕ*(D[1]*α[1] - D[2]*α[2])/(D[1]*α[1] + D[2]*α[2])
            ϕ = ϕ + exp(-α[1]*z0) - exp(-α[1]*(2*zb + z0))
            ϕ = ϕ/(2*D[1]*α[1])

        else

            ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
            ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
            ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
            ϕ = ϕ*exp(α[1]*z0 - α[1]*z)
            ϕ = ϕ - ((exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z))))/(2*D[1]*α[1]))

        end

    ## equation 12
    elseif z > z0 && z < l
        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

    #equation 13
    elseif z >= l

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    return ϕ
end

function fluence_spatial_complexfreq(z, s, μa, μsp, l, f, nmed)

    #get numerical instabillity at large s
    D = @. 1/(3(μsp + μa))
    ν = 29.9792345/nmed
    ω = 2*π*f
    α = @. sqrt((D*s^2 + μa + ω*im/ν)/D)
    A = 0.493
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    ## equation 11
    if z < z0

        # have to account for exp blowing up to inf
        if abs2(α[1]*z0 - α[1]*z) > 15
            ϕ = exp(α[1]*(z0 - 2*l)) - exp(α[1]*(-z0 - 2*l - zb)) - exp(α[1]*(z0 - 2*l - 2*zb)) + exp(α[1]*(-z0 - 2*l - 4*zb))
            ϕ = ϕ*(D[1]*α[1] - D[2]*α[2])/(D[1]*α[1] + D[2]*α[2])
            ϕ = ϕ + exp(-α[1]*z0) - exp(-α[1]*(2*zb + z0))
            ϕ = ϕ/(2*D[1]*α[1])

        else

            ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
            ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
            ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
            ϕ = ϕ*exp(α[1]*z0 - α[1]*z)
            ϕ = ϕ - ((exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z))))/(2*D[1]*α[1]))

        end

    ## equation 12
    elseif z > z0 && z < l
        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

    #equation 13
    elseif z >= l

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    return ϕ
end

function compute_complexintegral(ρ, z, s, μa, μsp, l, f, nmed)
    return fluence_spatial_complexfreq(z, s, μa, μsp, l, f, nmed)*s*besselj0(ρ*s)
end

function fluence_FD(ρ, z, μa, μsp, l, f, nmed)
    return (quadgk((s) -> compute_complexintegral(ρ, z, s, μa, μsp, l, f, nmed), 0, Inf, rtol=1e-6)[1])/2π
end

function fluence_TD(ρ, z, μa, μsp, l, nmed)
    fr = 0.19531 .* Vector(1:512)
	Rt = Array{Complex{Float64},1}(undef, length(fr))
    map!((f) -> fluence_FD(ρ, z, μa, μsp, l, f, nmed), Rt, fr)
end

# compute equation 14 integrand
function compute_integral(ρ, z, s, μa, μsp, l)
    return fluence_spatial_freq(z, s, μa, μsp, l)*s*besselj0(ρ*s)
end

# numerically integrate equation 14
function fluence_steadystate(ρ, z, μa, μsp, l)
    return (quadgk((s) -> compute_integral(ρ, z, s, μa, μsp, l), 0, Inf, rtol=1e-6)[1])/2π
 end

# Eqn 16 to compute reflectance (use auto diff)
function refl_steadystate(ρ, μa, μsp, l)
    f(z) = fluence_steadystate(ρ, z, μa, μsp, l)
    g(z) = ForwardDiff.derivative(z -> f(z), z)
    return 0.118*fluence_steadystate(ρ, 0, μa, μsp, l) + 0.306/(3*(μa[1] + μsp[1]))*g(0)
end

# Eqn 3 Kienle 97 for spatial fluence semi-infinite (use as comparison for l -> ∞)
function ss_compare(ρ, z, μa, μsp)
    D = 1/(3(μsp + μa))
    A = 0.493
    z0 = 1/(μsp + μa)
    zb = 2*A*D
    μe = sqrt(3*μa*(μa + μsp))

    ϕ = exp(-μe*sqrt((z - z0)^2 + ρ^2))/sqrt((z - z0)^2 + ρ^2)
    ϕ = ϕ - exp(-μe*sqrt((z + z0 + 2*zb)^2 + ρ^2))/sqrt((z + z0 + 2*zb)^2 + ρ^2)
    ϕ = ϕ/(4*π*D)
end

function refl_steadystate1(ρ, μa, μsp, l)
    f(z) = ss_compare(ρ, z, μa, μsp)
    g(z) = ForwardDiff.derivative(z -> f(z), z)
    return 0.118*ss_compare(ρ, 0, μa, μsp) + 0.306/(3*(μa[1] + μsp[1]))*g(0)
end

##reproduce Figure 2 Kienle 98 *seem to be slightly off by some factor
a = map((ρ) -> refl_steadystate(ρ, [0.02,0.01], [1.3,1.2], 2), 0.5:0.5:20)
b = map((ρ) -> refl_steadystate(ρ, [0.02,0.01], [1.3,0.7], 2), 0.5:0.5:20)

plot(0.5:0.5:20, a, yscale=:log10, lw = 2, label = "crosses", xlabel = "distance [mm]", ylabel = "reflectance")
plot!(0.5:0.5:20, b, yscale=:log10, lw = 2, label = "circles")

##reproduce Figure 3 Kienle 98
b = map((ρ) -> refl_steadystate(ρ, [0.005,0.01], [1.3,1.0], 6), 0.5:0.5:20)
a = map((ρ) -> refl_steadystate(ρ, [0.005,0.022], [1.3,1.0], 6), 0.5:0.5:20)

plot(0.5:0.5:20, a, yscale=:log10, lw = 2, label = "crosses", xlabel = "distance [mm]", ylabel = "reflectance")
plot!(0.5:0.5:20, b, yscale=:log10, lw = 2, label = "circles")

##compare to semi-infinite model from Kienle 97 as l gets large
a = map((ρ) -> fluence_steadystate(ρ, 0, [0.02,0.02], [1.3,1.3], 10), 0.5:0.5:20)
b = map((ρ) -> ss_compare(ρ, 0, 0.02, 1.3), 0.5:0.5:20)

plot(0.5:0.5:20, a, yscale=:log10, lw = 2, label = "2-layer", xlabel = "distance [mm]", ylabel = "reflectance")
plot!(0.5:0.5:20, b, yscale=:log10, lw = 2, label = "semi-inf", alpha = 0.7)

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
# Ост - "Е6" - восток, Норд - "N8" - север, Запад - "W4" - запад, Юг - "S2" - юг.
function main(r::Robot)
    path = go_to_left_bottom_corner_and_get_path(r)
    for i ∈ (Ost, Nord, West, Sud)
        mark_direction(r, i)
    end
    go_back(r, path)
end

function go_to_left_bottom_corner_and_get_path(r::Robot) ::Array
    path=[]
    while !isborder(r,Sud) || isborder(r,West)
        if !isborder(r,Sud)
            move!(r,Sud)
            push!(path,Nord)
        end
        if !isborder(r,West)
            move!(r,West)
            push!(path,Ost)
        end
    end
    return path
end

function mark_direction(r,side::HorizonSide)
    while !isborder(r,side)
        putmarker!(r)
        move!(r,side)
    end
    putmarker!(r)
end

function go_back(r,path::Array)
    n=length(path)
    while n>0
        move!(r,path[n])
        n=n-1
    end
end


text:
#const background_img = zeros(UInt8,640*480)

#const background=Ref{WT_Image}(WT_Image(1,640,480,C_NULL,pointer(background_img)))

#1st dim = w
#2nd dim = h
function JT_trace(iFrame::Int,image_data::AbstractArray{UInt8,2},background_img=zeros(UInt8,length(image_data)),background=Ref{WT_Image}(WT_Image(1,size(image_data,1),size(image_data,2),C_NULL,pointer(background_img))))

    img=Ref{WT_Image}(WT_Image(1,size(image_data,1),size(image_data,2),C_NULL,pointer(image_data)))

    pnseg=Ref{Int32}(0)

    data=ccall((:find_segments,libwhisk_path),Ptr{Whisker2},(Int32,Ref{WT_Image},Ref{WT_Image},Ref{Int32}),iFrame,img,background,pnseg)

    wts=Array{Whisker1,1}()

    for i=1:pnseg[]
        ww=Whisker1(unsafe_load(data,i))
        push!(wts,deepcopy(ww))
        #push!(wts,Whisker1(unsafe_load(data,i)))
    end

    ccall((:Free_Whisker_Seg_Vec,libwhisk_path),Void,(Ptr{Whisker2}, Int32),data,pnseg[])

    wts
end

function WT_trace(iFrame::Int,image_data::AbstractArray{UInt8,2},min_length::Int,pad_pos::Tuple{Float32,Float32},mask::BitArray{2})

    whiskers=JT_trace(iFrame,image_data)

    WT_length_constraint(whiskers,min_length)

    WT_reorder_whisker(whiskers,pad_pos)

    #Apply mask
    apply_mask(whiskers,mask,min_length)

    apply_roi(whiskers,pad_pos)

    eliminate_redundant(whiskers)

    whiskers
end

function get_JT_params()
    ccall((:Params,libwhisk_path),Ptr{JT_Params},())
end

function change_JT_param(f_name,value)
    xx=get_JT_params()
    jt=unsafe_load(xx,1)

    setfield!(jt,f_name,value)

    unsafe_store!(xx,jt,1)

    nothing
end

#returns measurements for single
function JT_measure(wt::Tracker,frame_id,whisker_num,face_axis='x')

    facex=round(Int32,wt.pad_pos[1])
    facey=round(Int32,wt.pad_pos[2])
    whisk_num = convert(Int32,1)
    whisk_array = Ref{Whisker2}(Whisker2(wt.all_whiskers[frame_id][whisker_num]))

    mm=ccall((:Whisker_Segments_Measure,libwhisk_path),Ptr{JT_Measurements},
    (Ref{Whisker2},Int32,Int32,Int32,Cuchar),whisk_array,whisk_num,facex,facey,face_axis)

    unsafe_load(mm,1)
end

function JT_measure(w::Whisker1,facex::Float32,facey::Float32,face_axis='x')
    JT_measure(w,round(Int32,facex),round(Int32,facey),face_axis)
end

function JT_measure(w::Whisker1,facex::Int32,facey::Int32,face_axis='x')
    #This should be an option depending on how the whisker is oriented

    whisk_num = convert(Int32,1)

    whisk_array = Ref{Whisker2}(Whisker2(w))

    mm=ccall((:Whisker_Segments_Measure,libwhisk_path),Ptr{JT_Measurements},
    (Ref{Whisker2},Int32,Int32,Int32,Cuchar),whisk_array,whisk_num,facex,facey,face_axis)

    unsafe_load(mm,1)
end

function JT_measure(wt::Tracker,w::Whisker1,face_axis='x')

    facex=round(Int32,wt.pad_pos[1])
    facey=round(Int32,wt.pad_pos[2])

    JT_measure(w,facex,facey,face_axis)
end

function get_JT_measurements(wt::Tracker,face_axis='x')

    mymeas=Array{JT_Measurements,1}(0)

    for i=1:length(wt.all_whiskers)
        for j=1:length(wt.all_whiskers[i])
            xx=JT_measure(wt,i,j,face_axis)
            push!(mymeas,xx)
        end
    end

    mymeas
end

#Finds the best length threshold to discriminate whiskers and non whiskers for a known number
#of whiskers to track
function JT_Best_Threshold(wt,table,num_whisk)

    n_rows=convert(Int32,size(wt.vid,3))
    n_cols=convert(Int32,0)

    #Ranges for length threshold (I'm basically making these giant and ignoring them)
    low=1.0
    high=600.0

    is_gt=convert(Int32,1)

    #Number of whiskers that should be in frame
    target_count=convert(Int32,num_whisk)

    ccall((:Measurements_Table_Estimate_Best_Threshold_For_Known_Count,WhiskerTracking.libwhisk_path),Float64,
        (Ref{WhiskerTracking.JT_Measurements},Int32,Int32,Float64,Float64,Int32,Int32),
        table,n_rows,n_cols,low,high,is_gt,target_count)
end

function JT_find_segments(img,h,th,s,facemask)

    ccall(Libdl.dlsym(libwhisk,:compute_seed_from_point_field_on_grid),Void,(Ref{WT_Image},Int32,Int32,Int32,Float32,Float32,Ref{WT_Image},Ref{WT_Image},Ref{WT_Image}),
    img,50,4,5,0.0f0,0.0f0,h,th,s)

    mystride = img[].height

    s_j=unsafe_wrap(Array,convert(Ptr{Float32},s[].array),(length(img),))
    h_j=unsafe_wrap(Array,h[].array,(length(img),))
    th_j=unsafe_wrap(Array,convert(Ptr{Float32},th[].array),(length(img),))
    for i=1:length(h_j)
        if h_j[i]>0
            th_j[i] = th_j[i] / convert(Float32,h_j[i])
        end
    end

    mask=falses(length(img))

    nseeds=0
    for i=1:length(h_j)
        if (s_j[i] > 0.80f0)&(!facemask[i]) #SEED_THRESH
            nseeds+=1
            mask[i]=true
        end
    end

    myline=Ref{Line_Params}(Line_Params(0.5f0,convert(Float32,0.0f0),2.0f0,0.0f0))

    pi_cycle = pi / 4.0f0 / 18.0f0
    myangle=0.0f0

    scores=zeros(Float32,nseeds)
    inds=find(mask)
    j=1

    for i=inds

        offset=0.5f0
        xdir=100.0f0 * convert(Float32,cos(th_j[i]))
        ydir=100.0f0 * convert(Float32,sin(th_j[i]))

        if xdir < 0.0f0
            myangle = round(atan(-1.0f0 * ydir, -1.0f0 * xdir) / pi_cycle) * pi_cycle
        else
            myangle = round(atan(ydir,xdir) / 18.0f0) * pi_cycle
        end

        width = 2.0f0

        myline[].angle=myangle

        scores[j]=ccall(Libdl.dlsym(libwhisk,:eval_line),Float32,(Ref{Line_Params},Ref{WT_Image},Int32),myline,img,convert(Int32,i-1))
        j+=1
    end

    score_ord=sortperm(scores)

    nsegs=0
    ss=Ref{Seed}(Seed(0,0,0,0))
    for i=score_ord
        ss[].xpnt=(inds[i]-1) % mystride
        ss[].ypnt=div(inds[i]-1,mystride)
        ss[].xdir=round(Int32,100 * cos(th_j[inds[i]]))
        ss[].ydir=round(Int32,100 * sin(th_j[inds[i]]))
        w=ccall(Libdl.dlsym(libwhisk,:trace_whisker),Ptr{Whisker2},(Ref{Seed},Ref{WT_Image}),ss,img)

        if w == C_NULL
            ss[].ydir=round(Int32,100 * cos(th_j[inds[i]]))
            ss[].xdir=round(Int32,100 * sin(th_j[inds[i]]))
            w=ccall(Libdl.dlsym(libwhisk,:trace_whisker),Ptr{Whisker2},(Ref{Seed},Ref{WT_Image}),ss,img)
        end

        if w != C_NULL
            nsegs+=1
        end
    end

    nsegs
end


text:
# BLAS/LAPACK dependency for HSL.jl

libblas = library_dependency("libblas")
liblapack = library_dependency("liblapack")

@static if Sys.isapple()
  provides(Homebrew.HB, "homebrew/core/openblas", libblas, os=:Darwin)
  provides(Homebrew.HB, "homebrew/core/openblas", liblapack, os=:Darwin)
end

provides(AptGet, "libblas-dev", libblas, os=:Linux)
provides(AptGet, "liblapack-dev", liblapack, os=:Linux)

# try to guess what BinDeps installed
blas_lib = "-lblas"
lapack_lib = "-llapack"
@static if Sys.isapple()
  hb_libopenblas_dir = joinpath(Homebrew.prefix("openblas"), "lib")
  if isdir(hb_libopenblas_dir)
    blas_lib = "-L$hb_libopenblas_dir -lopenblas"
    lapack_lib = "-L$hb_libopenblas_dir -lopenblas"
  end
end


text:
# Command: Get Element Property
"""
    element_property(element::Element, property::AbstractString)

Get Element Property.
"""
function element_property(element::Element, property::AbstractString)
    @unpack addr, id = element.session
    element_id = element.id
    response = HTTP.get(
        "$addr/session/$id/element/$element_id/property/$property",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end


text:
module TangentMps

using LinearAlgebra, KrylovKit, TensorOperations
const AbstractTensor3 = AbstractArray{T,3} where {T}
const AbstractTensor4 = AbstractArray{T,4} where {T}
const Tensor3 = Array{T,3} where {T}
const Tensor4 = Array{T,4} where {T}

export ibcleft, ibcright, vumpsstep, tdvpstep
export al_from_ac_and_c, ar_from_ac_and_c
export vl_from_al, vr_from_ar
export enlargestep, twosite_variance
export transfer_from_left, transfer_from_right
export mul_matrix_from_left, mul_matrix_from_right
export mul_operator_onsite, mul_operator_with_left, mul_operator_with_right
export sx2, isy2, sz2, sp2, sm2, id2
export sx3, isy3, sz3, sp3, sm3, id3
export blbq
export mixedcanonical

const χ = nothing

include("algorithms.jl")

include("transfer_matrix.jl")
include("effective_hamiltonian.jl")
include("array_functions.jl")
include("default_operators.jl")

include("bond_enlarge.jl")

end # module


text:
module SAASGP4

using Compat

# include("deps.jl")

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("SAASGP4 not properly installed. Please run Pkg.build(\"SAASGP4\")")
end

function __init__()
    mainHandle = mainInit()

    sgp4SetLicFilePath(normpath(joinpath(@__DIR__, "../deps/usr/lib/")))

    envInit(mainHandle)
    astroFuncInit(mainHandle)
    timeFuncInit(mainHandle)
    tleInit(mainHandle)
    sgp4Init(mainHandle)
end

# Main

"""
    mainInit()

Returns a handle which can be used to access the static global data set needed by the Standardized Astrodynamic Algorithms DLLs to communicate among themselves
"""
function mainInit()
    ccall((:DllMainInit, dllmain), Clonglong, ())
end
const dllMainInit = mainInit

"""
    getLastErrMsg()

Returns a character string describing the last error that occurred
"""
function getLastErrMsg()
    msg = repeat(" ", 128)
    ccall((:GetLastErrMsg, dllmain), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    getLastInfoMsg()

Returns a character string describing the last informational message that was recorded
"""
function getLastInfoMsg()
    msg = repeat(" ", 128)
    ccall((:GetLastInfoMsg, dllmain), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    mainGetInfo()

Returns information about the DllMain DLL
"""
function mainGetInfo()
    msg = repeat(" ", 128)
    ccall((:DllMainGetInfo, dllmain), Nothing, (Cstring,), msg)
    strip(msg)
end
const dllMainGetInfo = mainGetInfo

"""
    getInitDllNames()

Returns a list of names of the Standardized Astrodynamic Algorithms DLLs that were initialized successfully
"""
function getInitDllNames()
    msg = repeat(" ", 512)
    ccall((:GetInitDllNames, dllmain), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    openLogFile(path)

Opens a log file and enables the writing of diagnostic information into it
"""
function openLogFile(path)
    ccall((:OpenLogFile, dllmain), Cint, (Cstring,), path)
end

"""
    closeLogFile()

Closes the currently open log file
"""
function closeLogFile()
    ccall((:CloseLogFile, dllmain), Nothing, ())
end

# EnvConst

"""
    envInit(mainHandle::Clonglong)

Initializes the EnvInit DLL for use in the program
"""
function envInit(mainHandle::Clonglong)
    errCode = ccall((:EnvInit, envconst), Clonglong, (Clonglong,), mainHandle)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    envGetInfo()

Returns information about the EnvConst DLL
"""
function envGetInfo()
    msg = repeat(" ", 128)
    ccall((:EnvGetInfo, envconst), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    envGetGeoStr()

Returns the name of the current Earth constants (GEO) model
"""
function envGetGeoStr()
    msg = repeat(" ", 6)
    ccall((:EnvGetGeoStr, envconst),
        Nothing,
        (Cstring,),
        msg)
    strip(msg)
end

"""
    envSetGeoStr(geoStr)

Changes the Earth constants (GEO) setting to the model specified by a string literal
"""
function envSetGeoStr(geoStr)
    ccall((:EnvSetGeoStr, envconst),
        Nothing,
        (Cstring,),
        geoStr)
end

"Returns the current Earth constants (GEO) setting"
function envGetGeoIdx()
    ccall((:EnvGetGeoIdx, envconst),
        Cint,
        ())
end

"""
    envSetGeoIdx(id)

Changes the Earth constants (GEO) setting to the specified value
"""
function envSetGeoIdx(id)
    ccall((:EnvSetGeoIdx, envconst),
        Nothing,
        (Cint,),
        id)
end

"""
    envGetGeoConst(id)

Retrieves the value of one of the constants from the current Earth constants (GEO) model
"""
function envGetGeoConst(id)
    ccall((:EnvGetGeoConst, envconst),
        Cdouble,
        (Cint,),
        id)
end

"Returns the current fundamental catalogue (FK) setting"
function envGetFkIdx()
    ccall((:EnvGetFkIdx, envconst),
        Cint,
        ())
end

"""
    envSetFkIdx(id)

Changes the fundamental catalogue (FK) setting to the specified value
"""
function envSetFkIdx(id)
    ccall((:EnvSetFkIdx, envconst),
        Nothing,
        (Cint,),
        id)
end

"""
    envGetFkConst(id)

Retrieves the value of one of the constants from the current fundamental catalogue (FK) model
"""
function envGetFkConst(id)
    ccall((:EnvGetFkConst, envconst),
        Cdouble,
        (Cint,),
        id)
end

"""
    envGetFkPtr()

Returns a handle that can be used to access the fundamental catalogue (FK) data structure
"""
function envGetFkPtr()
    ccall((:EnvGetFkPtr, envconst),
        Clonglong,
        ())
end

"""
    envLoadFile(path)

Reads Earth constants (GEO) model and fundamental catalogue (FK) model settings from a file
"""
function envLoadFile(path)
    retval = ccall((:EnvLoadFile, envconst),
        Cint,
        (Cstring,),
        path)
    if retval != 0
        error(getLastErrMsg())
    end
end

"""
    envSaveFile(path, append = 0, format = 0)

Saves the current Earth constants (GEO) model and fundamental catalogue (FK) model settings to a file
"""
function envSaveFile(path, append = 0, format = 0)
    retval = ccall((:EnvSaveFile, envconst),
        Cint,
        (Cstring, Cint, Cint),
        path,
        append,
        format)
    if retval != 0
        error(getLastErrMsg())
    end
end

# AstroFunc

"""
    astroFuncInit(mainHandle)

Initializes AstroFunc DLL for use in the program
"""
function astroFuncInit(mainHandle)
    errCode = ccall((:AstroFuncInit, astrofunc), Clonglong, (Clonglong,), mainHandle)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    astroFuncGetInfo()

Returns information about the AstroFunc DLL
"""
function astroFuncGetInfo()
    msg = repeat(" ", 128)
    ccall((:AstroFuncGetInfo, astrofunc), Nothing, (Cstring,), msg)
    strip(msg)
end

# TimeFunc

"""
    timeFuncInit(mainHandle)

Initializes the TimeFunc DLL for use in the program
"""
function timeFuncInit(mainHandle)
    errCode = ccall((:TimeFuncInit, timefunc), Clonglong, (Clonglong,), mainHandle)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    timeFuncGetInfo()

Returns information about the TimeFunc DLL
"""
function timeFuncGetInfo()
    msg = repeat(" ", 128)
    ccall((:TimeFuncGetInfo, timefunc), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    DTGToUTC(dtg)

Converts a time in one of the DTG formats to a time in ds50UTC. DTG15, DTG17, DTG19, and DTG20 formats are accepted
"""
function DTGToUTC(dtg)
    ccall((:DTGToUTC, timefunc), Cdouble, (Cstring,), dtg)
end

# TLE

"""
    tleInit(mainHandle)

Initializes Tle DLL for use in the program
"""
function tleInit(mainHandle)
    errCode = ccall((:TleInit, tle), Clonglong, (Clonglong,), mainHandle)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    tleGetInfo()

Returns information about the Tle DLL
"""
function tleGetInfo()
    msg = repeat(" ", 128)
    ccall((:TleGetInfo, tle), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    tleAddSatFrLines(line1, line2)

Adds a TLE (satellite), using its directly specified first and second lines
"""
function tleAddSatFrLines(line1, line2)
    satkey = ccall((:TleAddSatFrLines, tle),
        Clonglong,
        (Cstring, Cstring),
        line1, line2)
end

"""
    tleRemoveSat(satkey)

Removes a TLE represented by the satKey from memory
"""
function tleRemoveSat(satkey)
    errCode = ccall((:TleRemoveSat, tle),
        Cint,
        (Clonglong,),
        satkey)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    tleRemoveAllSats()

Removes all the TLEs from memory
"""
function tleRemoveAllSats()
    errCode = ccall((:TleRemoveAllSats, tle),
        Cint,
        ())
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    tleLoadFile(path)

Loads TLEs (satellites) contained in a text file into the TLE DLL's binary tree
"""
function tleLoadFile(path)

    errCode = ccall((:TleLoadFile, tle),
        Cint,
        (Cstring,),
        path)
    if errCode != 0
        error(getLastErrMsg())
    end

end

"""
    tleGetCount()

Returns the number of TLEs currently loaded
"""
function tleGetCount()
    ccall((:TleGetCount, tle),
        Cint,
        ())
end

# SGP4Prop

"""
    sgp4Init(mainHandle)

Initializes the Sgp4 DLL for use in the program
"""
function sgp4Init(mainHandle)
    errCode = ccall((:Sgp4Init, sgp4prop), Clonglong, (Clonglong,), mainHandle)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    sgp4GetInfo()

Returns information about the Sgp4 DLL
"""
function sgp4GetInfo()
    msg = repeat(" ", 128)
    ccall((:Sgp4GetInfo, sgp4prop), Nothing, (Cstring,), msg)
    strip(msg)
end

"""
    sgp4SetLicFilePath(licFilePath)

Sets path to the Sgp4 Open License file if the license file is not in the current working folder
"""
function sgp4SetLicFilePath(licFilePath)
    ccall((:Sgp4SetLicFilePath, sgp4prop),
        Nothing,
        (Cstring,),
        licFilePath)
end

"""
    sgp4InitSat(satkey)

Initializes an SGP4 satellite from an SGP or SGP4 TLE
"""
function sgp4InitSat(satkey)
    errCode = ccall((:Sgp4InitSat, sgp4prop),
        Clong,
        (Clonglong,),
        satkey)
    if errCode != 0
        error(getLastErrMsg())
    end
end

"""
    sgp4RemoveSat(satkey)

Removes a satellite, represented by the satKey, from the set of satellites
"""
function sgp4RemoveSat(satkey)
    retval = ccall((:Sgp4RemoveSat, sgp4prop),
        Cint,
        (Clonglong,),
        satkey)
    if retval != 0
        error(getLastErrMsg())
    end
end

"""
    sgp4RemoveAllSats()

Removes all currently loaded satellites from memory
"""
function sgp4RemoveAllSats()
    retval = ccall((:Sgp4RemoveAllSats, sgp4prop),
        Cint,
        ())
    if retval != 0
        error(getLastErrMsg())
    end
end

"""
    sgp4PropDs50UTC(satkey, ds50UTC)

Propagates a satellite, represented by the satKey, to the time expressed in days since 1950, UTC
"""
function sgp4PropDs50UTC(satkey, ds50UTC)
    pos = Array{Cdouble, 1}(undef, 3)
    vel = Array{Cdouble, 1}(undef, 3)
    llh = Array{Cdouble, 1}(undef, 3)
    mse = Ref{Cdouble}(0.0)

    retval = ccall((:Sgp4PropDs50UTC, sgp4prop),
        Cint,
        (Clonglong, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
        satkey,     ds50UTC, mse,          pos,          vel,          llh)
    if retval != 0
        error(getLastErrMsg())
    end

    pos, vel, llh, mse[]
end

"""
    sgp4PropDs50UTC!(pos, vel, llh, mse, satkey, ds50UTC)

Propagates a satellite, represented by the satKey, to the time expressed in days since 1950, UTC
"""
function sgp4PropDs50UTC!(pos, vel, llh, mse, satkey, ds50UTC)
    retval = ccall((:Sgp4PropDs50UTC, sgp4prop),
        Cint,
        (Clonglong, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
        satkey,     ds50UTC, mse,          pos,          vel,          llh)
    if retval != 0
        error(getLastErrMsg())
    end
end

"""
    sgp4PropMse(satkey, mse)

Propagates a satellite, represented by the satKey, to the time expressed in minutes since the satellite's epoch time
"""
function sgp4PropMse(satkey, mse)
    pos = Array{Cdouble, 1}(undef, 3)
    vel = Array{Cdouble, 1}(undef, 3)
    llh = Array{Cdouble, 1}(undef, 3)
    ds50UTC = Ref{Cdouble}(0.0)

    retval = ccall((:Sgp4PropMse, sgp4prop),
        Cint,
        (Clonglong, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
        satkey,     mse,     ds50UTC,      pos,          vel,          llh)
    if retval != 0
        error(getLastErrMsg())
    end

    pos, vel, llh, ds50UTC[]
end

"""
    sgp4PropMse!(pos, vel, llh, ds50UTC, satkey, mse)

Propagates a satellite, represented by the satKey, to the time expressed in minutes since the satellite's epoch time.
Stores the results in the provided variables.

# Examples
```julia-repl
julia> pos = Array{Cdouble, 1}(undef, 3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> vel = Array{Cdouble, 1}(undef, 3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> llh = Array{Cdouble, 1}(undef, 3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> ds50UTC = Ref{Cdouble}(0.0)
Base.RefValue{Float64}(0.0)

julia> SAASGP4.sgp4PropMse!(pos, vel, llh, ds50UTC, satkey, mse)

julia> pos
3-element Array{Float64,1}:
 42155.0257401623
   980.5342672211367
    -9.307135121360039

julia> vel
3-element Array{Float64,1}:
 -0.07133558166329691
  3.07373221015029
 -0.0018402927526264539

julia> llh
3-element Array{Float64,1}:
    -0.012659361941212956
    34.46169610621692
 35788.29393717118

julia> ds50UTC[]
18314.489312984446
```
"""
function sgp4PropMse!(pos, vel, llh, ds50UTC, satkey, mse)
    retval = ccall((:Sgp4PropMse, sgp4prop),
        Cint,
        (Clonglong, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
        satkey,     mse,     ds50UTC,            pos,          vel,          llh)
    if retval != 0
        error(getLastErrMsg())
    end
end
"""
    sgp4GetPropOut(satkey, index)

Retrieves propagator's precomputed results. This function can be used to obtain results from a propagation
which are not made available through calls to the propagation functions themselves.
"""
function sgp4GetPropOut(satkey, index)
    if index ∉ 1:4
        error("Invalid index, must be ∈ 1..4")
    end

    datalength = index === 1 ? 1 : index === 2 ? 3 : 6
    data = Array{Cdouble, 1}(undef, datalength)

    retval = ccall((:Sgp4GetPropOut, sgp4prop),
        Cint,
        (Clonglong, Cint,  Ref{Cdouble}),
        satkey,     index, data)
    if retval != 0
        error(getLastErrMsg())
    end

    data
end

end # module


text:
# Trajectory Parsing

"""
    parseTrajectory(sList, vList, uList)

Stores a parsed trajectory by separating the positions (sList), the velocities
(vList), and the controls (uList).

See [`getParseTrajectory`](@ref) for more info.
"""
struct parseTrajectory
    sList
    vList
    uList
end

"""
    getParseTrajectory(traj, nDim)

Converts a trajectory of the form `[x0; u0; x1; u1; ...]` where `xk = [sk; vk]`
to a [`parseTrajectory`](@ref) struct.
"""
function getParseTrajectory(traj, nDim)
    sList = []
    vList = []
    uList = []

    kStart = 1
    maxK = length(traj)
    # println("Max = $maxK")

    while kStart < (maxK - 2 * nDim)
        # println("Accessing $kStart to $(kStart + 3 * nDim - 1)")
        push!(sList, traj[kStart:kStart + nDim - 1])
        push!(vList, traj[kStart + nDim:kStart + 2 * nDim - 1])
        push!(uList, traj[kStart + 2 * nDim:kStart + 3 * nDim - 1])

        kStart = kStart + 3 * nDim
    end

    # Add the final state (has no control step)
    # println("Last Access is $(kStart + 2 * nDim - 1)")
    push!(sList, traj[kStart:kStart + nDim - 1])
    push!(vList, traj[kStart + nDim:kStart + 2 * nDim - 1])

    return parseTrajectory(sList, vList, uList)
end

"""
    splitDimensions(stateList)

Converts a list of `n` dimensional states into `n` lists of 1 dimension.
"""
function splitDimensions(stateList)
    nDim = size(stateList[1], 1) # How many dimensions?

    splitStateList = []

    for d in 1:nDim
        push!(splitStateList, [x[d] for x in stateList])
    end

    return splitStateList
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
# Functions for representing graphs in GraphViz's dot format
# http://www.graphviz.org/
# http://www.graphviz.org/Documentation/dotguide.pdf
# http://www.graphviz.org/pub/scm/graphviz2/doc/info/lang.html

# orientated and based more or less on :
# https://github.com/JuliaAttic/OldGraphs.jl/blob/master/src/dot.jl
# https://github.com/tkf/ShowGraphviz.jl
# and only simply modified (Roettgermann, 12/21)

#typealias
const AttributeDict = Dict{Symbol,Vector{String}}

"""
    default_attributes(g)

Return `attrs::AttributeDict default` as plotting paramter for `GraphViz` Plot.

See attributes on http://www.graphviz.org/. The AttributeDict is of `Dict{Symbol,Vector{String}}` with Symbols `P` (Plot), `G` (Graph), `N` (Node) and `E` (Edges).
For example the following dict-Entry set the layout engine to `dot`.

    attr[:layout] = ["G", "dot"]

#### Arguments
- `g::AbstractSimpleWeightedGraph`: a graph representation to export
"""
function default_attributes(graph::AbstractSimpleWeightedGraph; node_label::Bool = true, edge_label::Bool = false)::AttributeDict
    directed = Graphs.is_directed(graph)
    n = nv(graph)
    #size_graph=minimum(5.0 +*/sqrt, 15.0)  # ToDo: automated scaler
    attr = AttributeDict(
        :weights => ["P", (edge_label) ? "true" : "false"],
        :largenet => ["P", "200"],
        :arrowsize => ["E", "0.5"],
        :arrowtype => ["E", "normal"],
        :center => ["G", "1"],
        :overlap => ["G", "scale"],
        :color => ["N", "Turquoise"],
        :concentrate => ["G", "true"],
        :fontsize => ["N", (node_label) ? ((n < 100) ? "7.0" : "5") : "1.0"],
        :width => ["N", (node_label) ? "0.25" : "0.20"],
        :height => ["N", (node_label) ? "0.25" : "0.20"],
        :fixedsize => ["N", "true"],
        :fontsize => ["E", (edge_label) ? "8.0" : "1.0"],
        :layout => ["G", (directed) ? "dot" : "neato"], # dot or neato
        :size => ["G", (n < 20) ? "3.0" : ((n < 100) ? "7.0" : "10")],
        :shape => ["N", (node_label) ? "circle" : "point"]
    )

    # modify attr_ if it is a large network:
    if n > parse(Int64, attr[:largenet][2])
        attr = _mod_attr_large_network!(attr)
    end

    return attr
end

# internal function to modify `GraphViz` plot-paramter for large graphs
function _mod_attr_large_network!(attrs::AttributeDict)
    attrs[:shape] = ["N", "point"]
    attrs[:color] = ["N", "black"]
    attrs[:fontsize] = ["G", "1"]
    attrs[:concetrate] = ["G", "true"]
    attrs[:layout] = ["G", "sfdp"]
    attrs[:weights] = ["P", "false"]
    return attrs
end

# internal function to get `G`raph, `E`dge and `N`ode relateted attributes:
function get_GNE_attributes(attrs::AttributeDict, gne::String)
    if !isempty(attrs)
        GNE_attrs = Dict()
        for key in keys(attrs)
            if contains(attrs[key][1], gne)
                GNE_attrs[key] = attrs[key]
            end
        end
        return GNE_attrs
    else
        return ""
    end

end

# internal function to get a suitable string out of the attribute dictionary:
function _parse_attributes(mat::AbstractSimpleWeightedGraph, attrs::AttributeDict, gne::String)

    gne_attrs = get_GNE_attributes(attrs, gne)
    str_attr::String = ""

    if gne == "N" # node attributes
        str_attr = string("[", join(map(a -> to_dot(a[1], a[2][2]), collect(gne_attrs)), ","))
    elseif gne == "G" # graph attributes
        for key in keys(attrs)
            if contains(attrs[key][1], "G")
                str_attr = str_attr * string(to_dot(key, attrs[key][2]), ";\n ")
            end
        end
    elseif gne == "E" # edge attributes
        str_attr = string("[", join(map(a -> to_dot(a[1], a[2][2]), collect(gne_attrs)), ","))
    end
    return str_attr
end

# internal functions to identify strings:
to_dot(sym::Symbol, value::String) = "$sym=$value"
graph_type_string(graph::AbstractSimpleWeightedGraph) = WeightedNetwork.is_directed(graph) ? "digraph" : "graph"
edge_op(graph::AbstractSimpleWeightedGraph) = WeightedNetwork.is_directed(graph) ? "->" : "--"

# internal function to get the dot representation of a graph as a string.
function to_dot(graph::AbstractSimpleWeightedGraph; attributes::AttributeDict = default_attributes(graph), path = [])
    str = IOBuffer()
    to_dot(graph, str, attributes; path = path)
    String(take!(str)) #takebuf_string(str)
end

# internal function DOT-Language representation:
function to_dot(mat::AbstractSimpleWeightedGraph, stream::IO, attrs::AttributeDict; path = [])

    # check if `weighted` and labeled:
    edge_label = false
    if haskey(attrs, :weights)
        attr_ = attrs[:weights][2]
        (attr_ == "true") ? edge_label = true : edge_label = false
    end
    # write DOT:
    write(stream, "$(graph_type_string(mat)) graphname {\n")
    G = "G"
    write(stream, " $(_parse_attributes(mat,attrs, G))\n")
    n_vertices = nv(mat)
    N = "N"
    for node = 1:n_vertices
        if !isempty(path) && !Base.isnothing(findfirst(isequal(node), path))
            color_node = ",color=red]" # hard-coded> ToDo: could be part of attribute_list
        else
            color_node = "]"
        end
        write(stream, " $node $(_parse_attributes(mat,attrs, N)) $color_node;\n")
    end
    for node = 1:n_vertices
        childs = WeightedNetwork.children(mat, node)
        E = "E"
        for kid in childs
            # if n_vertices > kid # seems to be wrong / to think about that.
            if !isempty(path) && !Base.isnothing(findfirst(isequal(node), path)) && !Base.isnothing(findfirst(isequal(kid), path))
                edge_node = ",color=red]" # hard-coded> ToDo: could be part of attribute_list
            else
                edge_node = "]"
            end

            if edge_label
                w = mat.weights[node, kid]
                write(stream, " $node $(edge_op(mat)) $kid $(_parse_attributes(mat,attrs, E)), xlabel=$w $edge_node;\n")
            else
                write(stream, " $node $(edge_op(mat)) $kid $(_parse_attributes(mat,attrs, E)) $edge_node;\n")
            end
            # end
        end

    end
    write(stream, "}\n")
    return stream
end

"""
    plot_graphviz(g, node_label= true, edge_label=false; path = [])

Render graph `g` in iJulia using `Graphviz` engines.

#### Arguments
- `g::AbstractSimpleWeightedGraph`: a graph representation to export
- `node_label::Bool`: if true all nodes are numberd fom 1:N (default = true)
- `edge_label::Bool`: if true all edges are labeled with their weights (default = false)
- (optional) `path = []`: Int-Array of nodes. Nodes and their edges are drawn in red color (i.e. shortest path)
- (optional) `scale = 1.0`: Scale your plot
- (optional) `landscape = false`: render landscape
"""
function plot_graphviz(g::AbstractSimpleWeightedGraph; node_label::Bool = true, edge_label::Bool = false,
    path = [],
    scale = 3.0,
    landscape = false)

    attributes = default_attributes(g; node_label = node_label, edge_label = edge_label)
    attributes[:size] = ["G", string(scale)]
    (edge_label) ? attributes[:forcelabels] = ["G", "true"] : nothing
    (landscape) ? attributes[:orientation] = ["G", "LR"] : nothing

    gv_dot = to_dot(g; attributes = attributes, path = path)
    plot_graphviz(gv_dot)
end

"""
    plot_graphviz(g, attributes; path = [])

Render graph `g` in **iJulia** using `Graphviz` engines.

#### Arguments
- `g::AbstractSimpleWeightedGraph`: a graph representation to export
- `attributes::AttributeDict`: Render with own set of plotting attributes (see http://www.graphviz.org/ for details)
- (optional) `path = []`: Int-Array of nodes. Nodes and their edges are drawn in red color (i.e. shortest path)
"""
function plot_graphviz(g::AbstractSimpleWeightedGraph, attributes::AttributeDict; path = [])
    gv_dot = to_dot(g; attributes = attributes, path = path)
    plot_graphviz(gv_dot)
end

function plot_graphviz(str::AbstractString)
    ShowGraphviz.CONFIG.dot_option = `-q` # do not warn in iJulia!
    ShowGraphviz.DOT(str)
end

"""
    write_dot_file(g, file; attributes, path)

Export graph `g` to DOT-Format and store it in file `file`.

#### Arguments
- `g::AbstractSimpleWeightedGraph`: a graph representation to export
- `filename::AbstractString`: the filename to store (i.e. "graph.dot")
- (optional) `attributes::AttributeDict`: Dot-Attributes like node color stored in a dictionary
- (optional) `path = []`: Int-Array of nodes. Nodes and their edges are drawn in red color (i.e. shortest path)
"""
function write_dot_file(graph::AbstractSimpleWeightedGraph, filename::AbstractString; attributes::AttributeDict = default_attributes(graph), path = [])
    open(filename, "w") do f
        to_dot(graph, f, attributes; path = path)
    end
end

"""
    read_dot_file(file)

Import graph from DOT-Format and store it in file `SimpleWeightedGraph` or `SimpleWeightedDiGraph``.
ToDo: ERROR-Handling if not a suitable DOT-File is not implemented

#### Arguments
- `file::AbstractString`: the filename of dot-file (i.e. "graph.dot")
"""
function read_dot_file(filename::AbstractString)
    # to count total lines in the file
    node_count = 0

    directed = false

    # get size of graph
    f = open(filename, "r")
    for line in readlines(f)
        (line_type, nodes, weight) = _read_dotline(line)
        if line_type == "digraph"
            directed = true
        elseif line_type == "node"
            node_count += 1
        end
    end
    close(f)

    adj = zeros(node_count - 1, node_count - 1)

    # get edges
    f = open(filename, "r")
    for line in readlines(f)
        (line_type, nodes, weight) = _read_dotline(line)
        if line_type == "edge"
            adj[nodes[1], nodes[2]] = weight
            if directed == false
                adj[nodes[2], nodes[1]] = weight
            end
        end
    end
    # close file
    close(f)

    if directed
        adj = adj'
        return SimpleWeightedDiGraph(adj)
    else
        return SimpleWeightedGraph(adj)
    end

end

function _read_dotline(str::String)

    tokens = collect(tokenize(str))
    start_options = false
    weight_identifier = false
    weight = 1.0
    nodes = (Int)[]
    line_type = "node"
    idx_LSQARE = 0

    # handle input and identify edge, node or graph_line
    for token in tokens

        if token.val == "digraph"
            line_type = "digraph"
            break
        elseif token.val == "graph"
            line_type = "graph"
            break
        elseif token.kind == Tokenize.Tokens.INTEGER
            if weight_identifier == true
                weight = parse(Float64, token.val)
                weight_identifier = false
            else
                push!(nodes, parse(Int64, token.val))
            end
        elseif (Tokenize.Tokens.exactkind(token) == Tokenize.Tokens.ANON_FUNC) || Tokenize.Tokens.exactkind(token) == Tokenize.Tokens.ERROR
            line_type = "edge"
        elseif token.kind == Tokenize.Tokens.IDENTIFIER
            (token.val == "xlabel") ? weight_identifier = true : weight_identifier = false
        elseif token.kind == Tokenize.Tokens.FLOAT && (weight_identifier == true)
            weight = parse(Float64, token.val)
            weight_identifier = false
        elseif token.kind == Tokenize.Tokens.LSQUARE
            start_options = true
            idx_LSQARE = token.startpos[2]
        elseif token.kind == Tokenize.Tokens.EQ && (start_options == false)
            line_type = "graph"
        end
    end

    return (line_type, nodes, weight)

end

max_stars_repo_name: tragisch/NetworkOptimization

max_issues_repo_name: tragisch/NetworkOptimization

max_forks_repo_name: tragisch/NetworkOptimization


text:
include_implicit_dynamics()

"""
    Cartpole w/ internal friction
"""

struct Cartpole{T}
    dim::Dimensions

    mc::T     # mass of the cart in kg
    mp::T     # mass of the pole (point mass at the end) in kg
    l::T      # length of the pole in m
    g::T      # gravity m/s^2
end

function kinematics(model::Cartpole, q)
    [q[1] + model.l * sin(q[2]); -model.l * cos(q[2])]
end

lagrangian(model::Cartpole, q, q̇) = 0.0

function M_func(model::Cartpole, x)
    H = @SMatrix [model.mc + model.mp model.mp * model.l * cos(x[2]);
				  model.mp * model.l * cos(x[2]) model.mp * model.l^2.0]
    return H
end

function B_func(model::Cartpole, x)
    @SMatrix [1.0;
              0.0]
end

function P_func(model::Cartpole, x)
    @SMatrix [1.0 0.0;
              0.0 1.0]
end

function C_func(model::Cartpole, q, q̇)
    # c_func(model, q, q̇) * q̇ - τ_func(model, q)
    C = @SMatrix [0.0 -1.0 * model.mp * q̇[2] * model.l * sin(q[2]);
	 			  0.0 0.0]
    G = @SVector [0.0,
				  model.mp * model.g * model.l * sin(q[2])]

    return -C * q̇ + G
end

model = Cartpole(Dimensions(2, 1, 0, 4),
    1.0, 0.2, 0.5, 9.81)

model_no_friction = Cartpole(Dimensions(2, 1, 0, 4),
    1.0, 0.2, 0.5, 9.81)

function lagrangian_derivatives(model::Cartpole, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

function dynamics(model::Cartpole, h, q0, q1, u1, λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2
		+ B_func(model, qm2) * u1
        + transpose(P_func(model, q2)) * λ1)
end

function dynamics_no_friction(model::Cartpole, h, q0, q1, u1, λ1, q2)
	# evalutate at midpoint
	qm1 = 0.5 * (q0 + q1)
    vm1 = (q1 - q0) / h[1]
    qm2 = 0.5 * (q1 + q2)
    vm2 = (q2 - q1) / h[1]

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	return (0.5 * h[1] * D1L1 + D2L1 + 0.5 * h[1] * D1L2 - D2L2
		+ B_func(model, qm2) * u1)
        # + transpose(P_func(model, q2)) * λ1)
end

function residual(model::Cartpole, z, θ, κ)

    nq = model.dim.q
    nu = model.dim.u
    nc = model.dim.c

    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    h = θ[2nq + nu .+ (1:1)]
    μ_slider = θ[2nq + nu + 1 .+ (1:1)]
    μ_angle = θ[2nq + nu + 1 + 1 .+ (1:1)]

    q2 = z[1:nq]
    c1 = z[nq .+ (1:nc)]
    c2 = z[nq + nc .+ (1:nc)]

    vT1 = (q2[1] - q1[1]) / h
    vT2 = (q2[2] - q1[2]) / h

    β1 = c1[1:2]
    β2 = c2[1:2]

    η1 = c1[2 .+ (1:2)]
    η2 = c2[2 .+ (1:2)]

    λ1 = [β1[2]; β2[2]]

    [
     dynamics(model, h, q0, q1, u1, λ1, q2);
     η1[2] - vT1[1];
     # β1[1] - (μ_slider[1] * 1.0); # static normal
     β1[1] .- μ_slider[1] * (model.mp + model.mc) * model.g * h[1];
     second_order_cone_product(η1, β1) - κ .* [1.0; 0.0];
     η2[2] - vT2[1];
     # β2[1] - (μ_angle[1] * 1.0); # static "normal"
     β2[1] .- μ_angle[1] * (model.mp * model.g * model.l) * h[1];
     second_order_cone_product(η2, β2) - κ .* [1.0; 0.0];
    ]
end

function residual_no_friction(model::Cartpole, z, θ, κ)

    nq = model.dim.q
    nu = model.dim.u
    nc = model.dim.c

    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    h = θ[2nq + nu .+ (1:1)]

    q2 = z[1:nq]

    return dynamics_no_friction(model, h, q0, q1, u1, zeros(nc), q2);
end

nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nz = nq + nc + nc
nz_no_friction = nq
nθ = nq + nq + nu + 1 + 2

# Declare variables
@variables z[1:nz]
@variables z_no_friction[1:nz_no_friction]
@variables θ[1:nθ]
@variables κ[1:1]

# Residual
r = vec(residual(model, z, θ, κ))
r = Symbolics.simplify.(r)
rz = Symbolics.jacobian(r, z, simplify = true)
rθ = Symbolics.jacobian(r, θ, simplify = true)

# Build function
r_func = eval(build_function(r, z, θ, κ)[2])
rz_func = eval(build_function(rz, z, θ)[2])
rθ_func = eval(build_function(rθ, z, θ)[2])

rz_array = similar(rz, Float64)
rθ_array = similar(rθ, Float64)

@save joinpath(@__DIR__, "dynamics/residual.jl") r_func rz_func rθ_func rz_array rθ_array
@load joinpath(@__DIR__, "dynamics/residual.jl") r_func rz_func rθ_func rz_array rθ_array

# Residual
r_no_friction = residual_no_friction(model, z_no_friction, θ, κ)
r_no_friction = Symbolics.simplify.(r_no_friction)
rz_no_friction = Symbolics.jacobian(r_no_friction, z_no_friction, simplify = true)
rθ_no_friction = Symbolics.jacobian(r_no_friction, θ, simplify = true)

# Build function
r_no_friction_func = eval(build_function(r_no_friction, z_no_friction, θ, κ)[2])
rz_no_friction_func = eval(build_function(rz_no_friction, z_no_friction, θ)[2])
rθ_no_friction_func = eval(build_function(rθ_no_friction, z_no_friction, θ)[2])

rz_no_friction_array = similar(rz_no_friction, Float64)
rθ_no_friction_array = similar(rθ_no_friction, Float64)

@save joinpath(@__DIR__, "dynamics/residual_no_friction.jl") r_no_friction_func rz_no_friction_func rθ_no_friction_func rz_no_friction_array rθ_no_friction_array
@load joinpath(@__DIR__, "dynamics/residual_no_friction.jl") r_no_friction_func rz_no_friction_func rθ_no_friction_func rz_no_friction_array rθ_no_friction_array

idx_ineq = collect(1:0)
idx_soc = [collect([nq .+ (1:2)]...), collect([nq + 2 .+ (1:2)]...), collect([nq + nc .+ (1:2)]...), collect([nq + nc + 2 .+ (1:2)]...)]

z_subset_init = [1.0; 0.1; 1.0; 0.1; 1.0; 0.1; 1.0; 0.1]

q0 = [0.0; 0.0]
q1 = [0.0; 0.0]
u1 = [0.0]
h = 0.1

z0 = copy([q1; ones(nc + nc)])
θ0 = [q0; q1; u1; h; 0.1; 0.1]

# options
opts = InteriorPointOptions(
   κ_init = 0.1,
   κ_tol = 1.0e-4,
   r_tol = 1.0e-8,
   diff_sol = true)

# solver
ip = interior_point(z0, θ0,
   r! = r_func, rz! = rz_func,
   rz = similar(rz, Float64),
   rθ! = rθ_func,
   rθ = similar(rθ, Float64),
   # idx_ineq = idx_ineq,
   idx_soc = idx_soc,
   opts = opts)

# simulate
T = 100
q_hist = [q0, q1]

for t = 1:T-1
    u1 = [t == 1 ? 1.0 : t == 2 ? 0.0 : 0.0]
    ip.z .= copy([q_hist[end]; z_subset_init])
    ip.θ .= copy([q_hist[end-1]; q_hist[end]; u1; h; 0.1; 0.05])
    status = interior_point_solve!(ip)

    if status
        push!(q_hist, ip.z[1:nq])
    else
        println("dynamics failure t = $t")
        println("res norm = $(norm(ip.r, Inf))")
        # break
    end
end

include(joinpath(@__DIR__, "visuals.jl"))
include(joinpath(pwd(), "models/visualize.jl"))

vis = Visualizer()
# open(vis)
render(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", 1)

visualize!(vis, model, q_hist, Δt = h)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
#=
If the numbers 1 to 5 are written out in words: one, two, three, four, five, then there are 3 + 3 + 5 + 4 + 4 = 19 letters used in total.

If all the numbers from 1 to 1000 (one thousand) inclusive were written out in words, how many letters would be used?

NOTE: Do not count spaces or hyphens. For example, 342 (three hundred and forty-two) contains 23 letters and 115 (one hundred and fifteen) contains 20 letters. The use of "and" when writing out numbers is in compliance with British usage.
=#

# This solution is not as clean as I'd like, but it is fast
dict = Dict(
  1 => "one",
  2 => "two",
  3 => "three",
  4 => "four",
  5 => "five",
  6 => "six",
  7 => "seven",
  8 => "eight",
  9 => "nine",
  10 => "ten",
  11 => "eleven",
  12 => "twelve",
  13 => "thirteen",
  14 => "fourteen",
  15 => "fifteen",
  16 => "sixteen",
  17 => "seventeen",
  18 => "eighteen",
  19 => "nineteen",
)

tens_digits = Dict(
  2 => "twenty",
	3 => "thirty",
	4 => "forty",
	5 => "fifty",
	6 => "sixty",
	7 => "seventy",
	8 => "eighty",
	9 => "ninety"
)

function number_name(value, result, tens = false)
  value_digits = reverse(digits(value))

  if value == 0
    return result
  end

  if length(value_digits) == 4
    result = result * dict[first(value_digits)]
    result = result * "thousand"
    number_name(rem(value,1000), result)
  elseif length(value_digits) == 3
    result = result * dict[first(value_digits)]
    result = result * "hundred"
    number_name(rem(value,100), result)
  elseif length(value_digits) <= 2
    if value >= 20
      if length(result) > 0
        result = result * "and"
      end
      result = result * tens_digits[first(value_digits)]
      number_name(rem(value,10),result, true)
    else
      if length(result) > 0 && tens == false
        result = result * "and"
      end
      result = result * dict[value]
      number_name(0, result)
    end
  end
end

function calc()
  total = 0
  for i in 1:1000
    total += length(number_name(i, ""))
  end
  total
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
"""
    create_geometry(setup; recreate = true)

Create cells, surfaces and finite element mesh. If `recreate = false`, previous geometry
will be reused.

This function does the following:
- Check geometry setup consistency
- Create or load cell configuration
- Create or load surface triangulation
- Call TetGen
- Deform domain
- Split mesh into compartments
For custom geometries with more than one compartment, call `split_mesh`
directly instead. This requires facet and element labels.
"""
function create_geometry(setup; recreate = true)
    (; refinement) = setup

    # File name for saving or loading geometry
    filename = joinpath(setup.meshdir, setup.name)

    # Make sure that folder exists
    dir = dirname(filename)
    isdir(dir) || mkpath(dir)

    if hasfield(typeof(setup), :ecs_shape)
        setup.ecs_shape ∈ [:no_ecs, :box, :convex_hull, :tight_wrap] ||
            error("Invalid ECS shape")
    end

    if setup isa Union{CylinderSetup,SphereSetup}
        # Check if cell description file is already available
        cellfilename = filename * "_cells"
        if isfile(cellfilename) && !recreate
            cells = read_cells(cellfilename)
        else
            cells = create_cells(setup)
            save_cells(cells, cellfilename)
        end
    else
        cells = nothing
    end

    # Make directory for storing finite elements mesh
    is_stl = endswith(filename, ".stl")
    if is_stl
        setup.ecs_shape == :no_ecs || error("ECS is only available for surface meshes")
        filename = filename[1:end-4]
    end
    save_meshdir_path = filename * "_dir"
    isdir(save_meshdir_path) || mkpath(save_meshdir_path)

    # Use an existing finite elements mesh or create a new finite elements mesh. The name of
    # the finite elements mesh is stored in the string `fname_tetgen_femesh`
    refinement_str = isinf(refinement) ? "" : "_refinement$refinement"
    fname_tetgen = "$save_meshdir_path/$(split(filename, "/")[end])$(refinement_str)_mesh"

    # Read or create surface triangulation
    if is_stl
        surfaces = nothing
    elseif isfile(fname_tetgen * ".node") && isfile(fname_tetgen * ".poly") && !recreate
        surfaces = read_surfaces(fname_tetgen)
    else
        if isa(setup, NeuronSetup)
            surfaces = create_surfaces(setup, filename)
        else
            surfaces = create_surfaces(setup, cells)
        end
        save_surfaces(fname_tetgen, surfaces)
    end

    # Add ".1" suffix to output file name, since this is what Tetgen does
    fname_tetgen_femesh = fname_tetgen * ".1"

    if isfile(fname_tetgen_femesh * ".node") && !recreate
        # Read global mesh from Tetgen output
        mesh_all = read_tetgen(fname_tetgen_femesh)
    elseif is_stl
        error("Not implemented")
    else
        mesh_all = call_tetgen(surfaces, refinement)
        save_tetgen(mesh_all, fname_tetgen_femesh)
    end

    # Deform domain
    if isa(setup, CylinderSetup) && (setup.bend > 1e-16 || setup.twist > 1e-16)
        @debug "Deforming domain with bend $(setup.bend) and twist $(setup.twist)"
        deform_domain!(mesh_all.points, setup.bend, setup.twist)
    end

    # Split mesh into compartments
    mesh = split_mesh(mesh_all)

    mesh, surfaces, cells
end


text:
# All code is attached to its underlying database source
struct SourceCode{Source}
    source::Source
    code::Expr
end

# Every time `SourceCode` objects are combined, check to see whether they all come from the same source
function pop_source!(sources, something)
    something
end
function pop_source!(sources, source_code::SourceCode)
    push!(sources, source_code.source)
    source_code.code
end
function key_pop_source!(sources, (key, source_code))
    code = pop_source!(souces, source_code)
    Expr(:kw, key, code)
end

function combine_sources(a_function, source_codes...; key_source_codes...)
    sources = Set(Any[])
    codes = partial_map(pop_source!, sources, source_codes)
    key_codes = (
        key_pop_source!(sources, key_source_code)
        for key_source_code in key_source_codes
    )
    if length(sources) != 1
        error("Expected exactly one source; got ($(sources...))")
    else
        SourceCode(
            first(sources),
            Expr(:call, a_function, Expr(:parameters, key_codes...), codes...)
        )
    end
end

# `@code_instead` hijacks call to create Julia expressions instead of evaluating functions
function numbered_argument(number)
    Symbol(string("argument", number))
end

function assert_type(argument, type)
    Expr(:(::), argument, type)
end

# Splat Varargs
function maybe_splat(argument, a_type)
    if @capture a_type Vararg{AType_}
        Expr(:(...), argument)
    else
        argument
    end
end

function code_instead(location, a_function, types...)
    arguments = ntuple(numbered_argument, length(types))
    keywords = Expr(:parameters, Expr(:..., :keywords))
    Expr(:function,
        Expr(:call, a_function, keywords, map_unrolled(assert_type, arguments, types)...),
        Expr(:block, location, Expr(:call,
            combine_sources,
            keywords,
            a_function,
            map_unrolled(maybe_splat, arguments, types)...
        ))
    )
end

macro code_instead(a_function, types...)
    code_instead(__source__, a_function, types...) |> esc
end

# Currently, query doesn't do anything
function QueryOperators.query(source_code::SourceCode)
    source_code
end

@code_instead (==) SourceCode Any
@code_instead (==) Any SourceCode
@code_instead (==) SourceCode SourceCode

@code_instead (!=) SourceCode Any
@code_instead (!=) Any SourceCode
@code_instead (!=) SourceCode SourceCode

@code_instead (!) SourceCode

@code_instead (&) SourceCode Any
@code_instead (&) Any SourceCode
@code_instead (&) SourceCode SourceCode

@code_instead (|) SourceCode Any
@code_instead (|) Any SourceCode
@code_instead (|) SourceCode SourceCode

@code_instead (*) SourceCode Any
@code_instead (*) Any SourceCode
@code_instead (*) SourceCode SourceCode

@code_instead (+) SourceCode Any
@code_instead (+) Any SourceCode
@code_instead (+) SourceCode SourceCode

@code_instead (%) SourceCode Any
@code_instead (%) Any SourceCode
@code_instead (%) SourceCode SourceCode

@code_instead abs SourceCode

@code_instead chop SourceCode

@code_instead chop SourceCode

@code_instead coalesce SourceCode Vararg{Any}

@code_instead QueryOperators.drop SourceCode Integer

@code_instead QueryOperators.filter SourceCode Any Expr

@code_instead QueryOperators.groupby SourceCode Any Expr Any Expr

@code_instead if_else SourceCode Any Any
@code_instead if_else Any SourceCode Any
@code_instead if_else Any Any SourceCode
@code_instead if_else Any SourceCode SourceCode
@code_instead if_else SourceCode Any SourceCode
@code_instead if_else SourceCode SourceCode Any
@code_instead if_else SourceCode SourceCode SourceCode

@code_instead in SourceCode Any
@code_instead in Any SourceCode
@code_instead in SourceCode SourceCode

@code_instead isequal SourceCode Any
@code_instead isequal Any SourceCode
@code_instead isequal SourceCode SourceCode

@code_instead isless SourceCode Any
@code_instead isless Any SourceCode
@code_instead isless SourceCode SourceCode

@code_instead ismissing SourceCode

@code_instead QueryOperators.join SourceCode SourceCode Any Expr Any Expr Any Expr

@code_instead max SourceCode Vararg{Any}
@code_instead min SourceCode Vararg{Any}

@code_instead length SourceCode

@code_instead QueryOperators.map SourceCode Any Expr

@code_instead occursin AbstractString SourceCode
@code_instead occursin Regex SourceCode

@code_instead QueryOperators.orderby SourceCode Any Expr

@code_instead QueryOperators.orderby_descending SourceCode Any Expr

@code_instead secondary SourceCode

@code_instead startswith SourceCode Any
@code_instead startswith Any SourceCode
@code_instead startswith SourceCode SourceCode

@code_instead string SourceCode Any
@code_instead string Any SourceCode
@code_instead string SourceCode SourceCode

@code_instead QueryOperators.take SourceCode Any

@code_instead QueryOperators.thenby SourceCode Any Expr

@code_instead QueryOperators.thenby_descending SourceCode Any Expr

@code_instead QueryOperators.unique SourceCode Any Expr

@code_instead uppercase SourceCode


text:
using CategoricalArrays

# CategoricalArrays
for ord in (true, false)
    v = CategoricalArray(repeat(["b", "a"], inner = 5), ordered=ord)
    v2 = rcopy(CategoricalArray, RObject(v))
    @test isequal(v2, v)
    @test v2 isa CategoricalVector{String}
    @test levels(v2) == levels(v)
    @test isordered(v2) === ord
end
@test levels(rcopy(CategoricalArray, R"factor(c('c','a','a'))")) == ["a","c"]
@test !isordered(rcopy(CategoricalArray, R"factor(c('c','a','a'))"))
@test isordered(rcopy(CategoricalArray, R"ordered(c('c','a','a'))"))

a = Array{Union{String, Missing}}(repeat(["b", "a"], inner = 5))
a[repeat([true, false], outer = 5)] .= missing
for ord in (true, false)
    v = CategoricalArray(a, ordered=ord)
    v2 = rcopy(CategoricalArray, RObject(v))
    @test isequal(v2, v)
    @test v2 isa CategoricalVector{Union{String,Missing}}
    @test levels(v2) == levels(v)
    @test isordered(v2) === ord
end
@test levels(rcopy(CategoricalArray, R"factor(c('c',NA,'a'))")) == ["a","c"]
@test !isordered(rcopy(CategoricalArray, R"factor(c('c',NA,'a'))"))
@test isordered(rcopy(CategoricalArray, R"ordered(c('c',NA,'a'))"))


text:
module TreeKnit

# External modules
using Comonicon
using Logging
using LoggingExtras
using Parameters
using Random
using Setfield
# Personal modules
using TreeTools

include("mcc_base.jl")
export naive_mccs

include("mcc_splits.jl")

include("mcc_tools.jl")

include("mcc_IO.jl")
export read_mccs, write_mccs

include("resolving.jl")
export resolve!

include("SplitGraph/SplitGraph.jl")
using TreeKnit.SplitGraph

include("objects.jl")
export OptArgs

include("main.jl")
export computeMCCs, inferARG

include("SimpleReassortmentGraph/SimpleReassortmentGraph.jl")
import TreeKnit.SimpleReassortmentGraph: SRG
export SRG

include("cli.jl")

# TreeTools re-exports for docs
import TreeTools: node2tree, parse_newick
export node2tree
export parse_newick

end # module

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
using CUDA
using Adapt
#using Random
include("../util/Matrix_utils.jl")
include("../util/Enum_Collection.jl")
include("Button.jl")
include("Drawables.jl")

function get_window(name::Enum,windows::All_Windows)#::Window
    if name == main_window
        return windows.main_window
    else
        return windows.main_window
    end
end

function draw_self(window::Window,windows::All_Windows,buttons::All_Buttons)
    #kernel_blit_image_inplace(threadIdx().x,window.current_matrix,window.matrix_self,window.x_Coord,window.y_Coord)
    kernel_blit_image_inplace(threadIdx().x,blockIdx().x,window.current_matrix,window.matrix_self,window.x_Coord,window.y_Coord)

    for button in window.buttons
        current_button = get_button(button,buttons)
        draw_self(current_button,window.current_matrix)
    end
end

function click(x,y,window::Window,all_buttons::All_Buttons)
    click_on_window = false

    if includes_point(x,y,window.x_Coord,window.y_Coord,window.width,window.height)
        click_on_window = true
        #for button in window.buttons
        for button in window.buttons

            current_button = get_button(button,all_buttons)
            reward,click_on_child,coords_x,coords_y = click(x,y,current_button,window.x_Coord,window.y_Coord)
            sync_threads()
            if click_on_child
                if current_button.status[1,blockIdx().x]
                    #kernel_blit_image_inplace(threadIdx().x,blockIdx().x,window.current_matrix,current_button.matrix_clicked,coords_x,coords_y)
                else
                    #kernel_blit_image_inplace(threadIdx().x,blockIdx().x,window.current_matrix,current_button.matrix_unclicked,coords_x,coords_y)
                end
                return reward,click_on_window,window.x_Coord,window.y_Coord
            end
        end
    end
    return 0,false,0,0
end


text:
using Test
using LinearAlgebra

points2D = [SVector(0.1, 0.1), #1
            SVector(1.0, 1.0), #2
            SVector(0.2, 0.2), #3
            SVector(1.1, 1.1)] #4

root = create_tree(points2D, KMeansTreeOptions())

@test length(root.children) == 2
@test root.children[1].data[1] == 1
@test root.children[1].data[2] == 3
@test root.children[2].data[1] == 2
@test root.children[2].data[2] == 4
@test root.children[1].radius ≈ norm([0.05 0.05]) atol=1e-15

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
using LightGraphs.SimpleGraphs: SimpleEdge, AbstractSimpleGraph
using LightGraphs.SimpleGraphs: SimpleGraphs

export color, colors

mutable struct ColoredGraph{SGT<:AbstractSimpleGraph,T} <: AbstractColoredGraph{T}
    uncolored_graph::SGT
    colmap::Dict{SimpleEdge{T}, Vector{T}}
end

function ColoredGraph(g::AbstractSimpleGraph)
    T=eltype(g)
    cmap = Dict{SimpleEdge{T}, Vector{T}}()
    for e=edges(g)
        cmap[e] = [1]
        if !is_directed(g)
            cmap[reverse(e)] = [1]
        end
    end
    return ColoredGraph(g, cmap)
end

ColoredGraph(n::T=0) where T = ColoredGraph(SimpleGraph(n), Dict{SimpleEdge{T}, Vector{T}}())
ColoredDiGraph(n::T=0) where T = ColoredGraph(SimpleDiGraph(n), Dict{SimpleEdge{T}, Vector{T}}())

Base.eltype(x::ColoredGraph{SGT,T}) where {SGT,T} = T
LightGraphs.edgetype(x::ColoredGraph{SGT,T}) where {SGT,T}  = ColoredEdge{T}

LightGraphs.nv(g::AbstractColoredGraph) = nv(g.uncolored_graph)
LightGraphs.ne(g::AbstractColoredGraph) = begin
    tot = sum(length.(values(colmap(g))))
    is_directed(g) && return tot
    return Int(tot/2)
end
LightGraphs.vertices(g::AbstractColoredGraph) = vertices(g.uncolored_graph)
LightGraphs.edges(g::AbstractColoredGraph) = ColoredEdgeIter(g)

colmap(g::AbstractColoredGraph) = g.colmap
uncolored(g::AbstractColoredGraph) = g.uncolored_graph

set_colors!(g::AbstractColoredGraph, src, dst, cols::AbstractVector) = begin
    colmap(g)[SimpleEdge(src,dst)] = cols
    !is_directed(g) ? colmap(g)[SimpleEdge(dst, src)] = cols : nothing;
    return nothing
end

add_color!(g::AbstractColoredGraph, src, dst, col) = begin
    edge = SimpleEdge(src,dst)

    if edge ∉ keys(colmap(g))
        colmap(g)[edge] = Int[]
    end
    edge_cols = colmap(g)[edge]
    col ∈ edge_cols && return false
    push!(edge_cols, col)

    # If is not directed add the reverse connection
    if !is_directed(g)
        edge = SimpleEdge(dst,src)
        if edge ∉ keys(colmap(g))
            colmap(g)[edge] = Int[]
        end
        edge_cols = colmap(g)[edge]
        push!(edge_cols, col)
    end
    return true
end

colors(g::AbstractColoredGraph, src, dst) = colors(g, SimpleEdge(src,dst))
colors(g::AbstractColoredGraph, e) = colmap(g)[e]
colored(g::AbstractColoredGraph, e::SimpleEdge, col_id) = ColoredEdge(e, colors(g, e)[col_id])
colored(g::AbstractColoredGraph, src, dst, col)         = ColoredEdge(src, dst, colors(g, src, dst)[col])

LightGraphs.has_vertex(g::AbstractColoredGraph, v::Integer) = v in vertices(g)

# vertex properties
LightGraphs.inneighbors(g::AbstractColoredGraph, v::Integer) = inneighbors(uncolored(g),v)
LightGraphs.outneighbors(g::AbstractColoredGraph, v::Integer) = outneighbors(uncolored(g),v)

# specific implementations
#LightGraphs.is_directed(g::ColoredGraph) = is_directed(uncolored(g))
LightGraphs.is_directed(::Type{<:ColoredGraph{SGT,T}}) where {SGT, T}= is_directed(SGT)
LightGraphs.has_edge(g::ColoredGraph, s, d) = has_edge(uncolored(g), s, d)

function LightGraphs.has_edge(g::ColoredGraph, s, d, c)
    has_edge(uncolored(g), s, d) || return false
    return c ∈ colors(g, s, d)
end
LightGraphs.has_edge(g::ColoredGraph, e::ColoredEdge) = has_edge(g, src(e), dst(e), color(e))

# modifications
SimpleGraphs.add_edge!(g::AbstractColoredGraph, e::AbstractColoredEdge) =
    add_edge!(g, src(e), dst(e), color(e))
SimpleGraphs.add_edge!(g::AbstractColoredGraph, src, dst, col) = begin
    add_edge!(uncolored(g), src, dst)
    return add_color!(g, src, dst, col)
end

SimpleGraphs.rem_edge!(g::AbstractColoredGraph, e) = rem_edge!(g, src(e), dst(e))
SimpleGraphs.rem_edge!(g::AbstractColoredGraph, src, dst) = begin
    !rem_edge!(uncolored(g))  && return false
    delete!(colmap(g), SimpleEdge(src,dst))
    !is_directed(g) && delete!(colmap(g), SimpleEdge(dst, src))
    return true
end

SimpleGraphs.add_vertex!(g::AbstractColoredGraph) = add_vertex!(uncolored(g))
SimpleGraphs.add_vertices!(g::AbstractColoredGraph, n) = add_vertices!(uncolored(g), n)

SimpleGraphs.rem_vertex!(g::AbstractColoredGraph, v) = begin
    dsts = outneighbors(g, v)
    !rem_vertex!(uncolored(g)) && return false
    for dst=dsts
        delete!(colmap(g), SimpleEdge(v,dst))
    end
end

# I/O
show(io::IO, g::ColoredGraph) = print(io, "Colored Graph with $(nv(g)) vertices, $(ne(g)) edges")