text:
## Code for CategoricalArray

import Base: Array, convert, collect, copy, getindex, setindex!, similar, size,
             unique, vcat, in, summary, float, complex, copyto!

# Used for keyword argument default value
_isordered(x::AbstractCategoricalArray) = isordered(x)
_isordered(x::Any) = false

function reftype(sz::Int)
    if sz <= typemax(UInt8)
        return UInt8
    elseif sz <= typemax(UInt16)
        return UInt16
    elseif sz <= typemax(UInt32)
        return UInt32
    else
        return UInt64
    end
end

"""
    CategoricalArray{T}(undef, dims::Dims; ordered::Bool=false)
    CategoricalArray{T}(undef, dims::Int...; ordered::Bool=false)

Construct an uninitialized `CategoricalArray` with levels of type `T` and dimensions `dim`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

    CategoricalArray{T, N, R}(undef, dims::Dims; ordered::Bool=false)
    CategoricalArray{T, N, R}(undef, dims::Int...; ordered::Bool=false)

Similar to definition above, but uses reference type `R` instead of the default type
(`$DefaultRefType`).

    CategoricalArray(A::AbstractArray; ordered::Bool=false)

Construct a new `CategoricalArray` with the values from `A` and the same element type.

If the element type supports it, levels are sorted in ascending order;
else, they are kept in their order of appearance in `A`. The `ordered` keyword
argument determines whether the array values can be compared according to the
ordering of levels or not (see [`isordered`](@ref)).

    CategoricalArray(A::CategoricalArray; ordered::Bool=false)

If `A` is already a `CategoricalArray`, its levels are preserved;
the same applies to the ordered property and the reference type unless
explicitly overriden.
"""
function CategoricalArray end

"""
    CategoricalVector{T}(undef, m::Int; ordered::Bool=false)

Construct an uninitialized `CategoricalVector` with levels of type `T` and dimensions `dim`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

    CategoricalVector{T, R}(undef, m::Int; ordered::Bool=false)

Similar to definition above, but uses reference type `R` instead of the default type
(`$DefaultRefType`).

    CategoricalVector(A::AbstractVector; ordered::Bool=false)

Construct a `CategoricalVector` with the values from `A` and the same element type.

If the element type supports it, levels are sorted in ascending order;
else, they are kept in their order of appearance in `A`. The `ordered` keyword
argument determines whether the array values can be compared according to the
ordering of levels or not (see [`isordered`](@ref)).

    CategoricalVector(A::CategoricalVector; ordered::Bool=false)

If `A` is already a `CategoricalVector`, its levels are preserved;
the same applies to the ordered property and the reference type unless
explicitly overriden.
"""
function CategoricalVector end

"""
    CategoricalMatrix{T}(undef, m::Int, n::Int; ordered::Bool=false)

Construct an uninitialized `CategoricalMatrix` with levels of type `T` and dimensions `dim`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

    CategoricalMatrix{T, R}(undef, m::Int, n::Int; ordered::Bool=false)

Similar to definition above, but uses reference type `R` instead of the default type
(`$DefaultRefType`).

    CategoricalMatrix(A::AbstractVector; ordered::Bool=false)

Construct a `CategoricalMatrix` with the values from `A` and the same element type.

If the element type supports it, levels are sorted in ascending order;
else, they are kept in their order of appearance in `A`. The `ordered` keyword
argument determines whether the array values can be compared according to the
ordering of levels or not (see [`isordered`](@ref)).

    CategoricalMatrix(A::CategoricalMatrix; ordered::Bool=isordered(A))

If `A` is already a `CategoricalMatrix`, its levels are preserved;
the same applies to the ordered property and the reference type unless
explicitly overriden.
"""
function CategoricalMatrix end

# UndefInitializer array constructors

CategoricalArray(::UndefInitializer, dims::Int...; ordered=false) =
    CategoricalArray{String}(undef, dims, ordered=ordered)

function CategoricalArray{T, N, R}(::UndefInitializer, dims::NTuple{N,Int};
                                   ordered=false) where {T, N, R}
    U = leveltype(nonmissingtype(T))
    S = T >: Missing ? Union{U, Missing} : U
    V = CategoricalValue{U, R}
    CategoricalArray{S, N}(zeros(R, dims), CategoricalPool{U, R, V}(ordered))
end

CategoricalArray{T, N}(::UndefInitializer, dims::NTuple{N,Int}; ordered=false) where {T, N} =
    CategoricalArray{T, N, DefaultRefType}(undef, dims, ordered=ordered)
CategoricalArray{T, N}(::UndefInitializer, dims::NTuple{N,Int}; ordered=false) where
    {R, T <: Union{Missing, CategoricalValue{<:Any, R}}, N} =
    CategoricalArray{T, N, R}(undef, dims, ordered=ordered)
CategoricalArray{T}(::UndefInitializer, dims::NTuple{N,Int}; ordered=false) where {T, N} =
    CategoricalArray{T, N}(undef, dims, ordered=ordered)
CategoricalArray{T, 1}(::UndefInitializer, m::Int; ordered=false) where {T} =
    CategoricalArray{T, 1}(undef, (m,), ordered=ordered)
CategoricalArray{T, 2}(::UndefInitializer, m::Int, n::Int; ordered=false) where {T} =
    CategoricalArray{T, 2}(undef, (m, n), ordered=ordered)
CategoricalArray{T, 1, R}(::UndefInitializer, m::Int; ordered=false) where {T, R} =
    CategoricalArray{T, 1, R}(undef, (m,), ordered=ordered)
# R <: Integer is required to prevent default constructor from being called instead
CategoricalArray{T, 2, R}(::UndefInitializer, m::Int, n::Int; ordered=false) where {T, R <: Integer} =
    CategoricalArray{T, 2, R}(undef, (m, n), ordered=ordered)
CategoricalArray{T, 3, R}(::UndefInitializer, m::Int, n::Int, o::Int; ordered=false) where {T, R} =
    CategoricalArray{T, 3, R}(undef, (m, n, o), ordered=ordered)
CategoricalArray{T}(::UndefInitializer, m::Int; ordered=false) where {T} =
    CategoricalArray{T}(undef, (m,), ordered=ordered)
CategoricalArray{T}(::UndefInitializer, m::Int, n::Int; ordered=false) where {T} =
    CategoricalArray{T}(undef, (m, n), ordered=ordered)
CategoricalArray{T}(::UndefInitializer, m::Int, n::Int, o::Int; ordered=false) where {T} =
    CategoricalArray{T}(undef, (m, n, o), ordered=ordered)

CategoricalVector(::UndefInitializer, m::Integer; ordered=false) =
    CategoricalArray(undef, m, ordered=ordered)
CategoricalVector{T}(::UndefInitializer, m::Int; ordered=false) where {T} =
    CategoricalArray{T}(undef, (m,), ordered=ordered)

CategoricalMatrix(::UndefInitializer, m::Int, n::Int; ordered=false) =
    CategoricalArray(undef, m, n, ordered=ordered)
CategoricalMatrix{T}(::UndefInitializer, m::Int, n::Int; ordered=false) where {T} =
    CategoricalArray{T}(undef, (m, n), ordered=ordered)

## Constructors from arrays

function CategoricalArray{T, N, R}(A::CategoricalArray{S, N, Q};
                                   ordered=_isordered(A)) where {S, T, N, Q, R}
    V = unwrap_catvaluetype(T)
    res = convert(CategoricalArray{V, N, R}, A)
    refs = res.refs === A.refs ? copy(res.refs) : res.refs
    pool = res.pool === A.pool ? copy(res.pool) : res.pool
    ordered!(CategoricalArray{V, N}(refs, pool), ordered)
end

function CategoricalArray{T, N, R}(A::AbstractArray; ordered=_isordered(A)) where {T, N, R}
    V = unwrap_catvaluetype(T)
    ordered!(convert(CategoricalArray{V, N, R}, A), ordered)
end

# From AbstractArray
CategoricalArray{T, N}(A::AbstractArray{S, N}; ordered=_isordered(A)) where {S, T, N} =
    CategoricalArray{T, N, DefaultRefType}(A, ordered=ordered)
CategoricalArray{T}(A::AbstractArray{S, N}; ordered=_isordered(A)) where {S, T, N} =
    CategoricalArray{T, N}(A, ordered=ordered)
CategoricalArray(A::AbstractArray{T, N}; ordered=_isordered(A)) where {T, N} =
    CategoricalArray{T, N}(A, ordered=ordered)

CategoricalVector{T}(A::AbstractVector{S}; ordered=_isordered(A)) where {S, T} =
    CategoricalArray{T, 1}(A, ordered=ordered)
CategoricalVector(A::AbstractVector{T}; ordered=_isordered(A)) where {T} =
    CategoricalArray{T, 1}(A, ordered=ordered)

CategoricalMatrix{T}(A::AbstractMatrix{S}; ordered=_isordered(A)) where {S, T} =
    CategoricalArray{T, 2}(A, ordered=ordered)
CategoricalMatrix(A::AbstractMatrix{T}; ordered=_isordered(A)) where {T} =
    CategoricalArray{T, 2}(A, ordered=ordered)

# From CategoricalArray (preserve R)
CategoricalArray{T, N}(A::CategoricalArray{S, N, R};
                       ordered=_isordered(A)) where {S, T, N, R} =
    CategoricalArray{T, N, R}(A, ordered=ordered)
CategoricalArray{T}(A::CategoricalArray{S, N, R};
                    ordered=_isordered(A)) where {S, T, N, R} =
    CategoricalArray{T, N, R}(A, ordered=ordered)
CategoricalArray(A::CategoricalArray{T, N, R};
                 ordered=_isordered(A)) where {T, N, R} =
    CategoricalArray{T, N, R}(A, ordered=ordered)

CategoricalVector{T}(A::CategoricalArray{S, 1, R};
                     ordered=_isordered(A)) where {S, T, R} =
    CategoricalArray{T, 1, R}(A, ordered=ordered)
CategoricalVector(A::CategoricalArray{T, 1, R};
                  ordered=_isordered(A)) where {T, R} =
    CategoricalArray{T, 1, R}(A, ordered=ordered)

CategoricalMatrix{T}(A::CategoricalArray{S, 2, R};
                     ordered=_isordered(A)) where {S, T, R} =
    CategoricalArray{T, 2, R}(A, ordered=ordered)
CategoricalMatrix(A::CategoricalArray{T, 2, R};
                  ordered=_isordered(A)) where {T, R} =
    CategoricalArray{T, 2, R}(A, ordered=ordered)

## Conversion methods

# From AbstractArray
convert(::Type{CategoricalArray{T, N}}, A::AbstractArray{S, N}) where {S, T, N} =
    convert(CategoricalArray{T, N, DefaultRefType}, A)
convert(::Type{CategoricalArray{T}}, A::AbstractArray{S, N}) where {S, T, N} =
    convert(CategoricalArray{T, N}, A)
convert(::Type{CategoricalArray}, A::AbstractArray{T, N}) where {T, N} =
    convert(CategoricalArray{T, N}, A)

convert(::Type{CategoricalVector{T}}, A::AbstractVector) where {T} =
    convert(CategoricalVector{T, DefaultRefType}, A)
convert(::Type{CategoricalVector}, A::AbstractVector{T}) where {T} =
    convert(CategoricalVector{T}, A)
convert(::Type{CategoricalVector{T}}, A::CategoricalVector{T}) where {T} = A
convert(::Type{CategoricalVector}, A::CategoricalVector) = A

convert(::Type{CategoricalMatrix{T}}, A::AbstractMatrix) where {T} =
    convert(CategoricalMatrix{T, DefaultRefType}, A)
convert(::Type{CategoricalMatrix}, A::AbstractMatrix{T}) where {T} =
    convert(CategoricalMatrix{T}, A)
convert(::Type{CategoricalMatrix{T}}, A::CategoricalMatrix{T}) where {T} = A
convert(::Type{CategoricalMatrix}, A::CategoricalMatrix) = A

function convert(::Type{CategoricalArray{T, N, R}}, A::AbstractArray{S, N}) where {S, T, N, R}
    res = CategoricalArray{T, N, R}(undef, size(A))
    copyto!(res, A)

    # if order is defined for level type, automatically apply it
    L = leveltype(res)
    if hasmethod(isless, Tuple{L, L})
        levels!(res.pool, sort(levels(res.pool)))
    end

    res
end

# From CategoricalArray (preserve levels, ordering and R)
function convert(::Type{CategoricalArray{T, N, R}}, A::CategoricalArray{S, N}) where {S, T, N, R}
    if length(A.pool) > typemax(R)
        throw(LevelsException{T, R}(levels(A)[typemax(R)+1:end]))
    end

    if !(T >: Missing) && S >: Missing && any(iszero, A.refs)
        throw(MissingException("cannot convert CategoricalArray with missing values to a CategoricalArray{$T}"))
    end

    pool = convert(CategoricalPool{unwrap_catvaluetype(nonmissingtype(T)), R}, A.pool)
    refs = convert(Array{R, N}, A.refs)
    CategoricalArray{unwrap_catvaluetype(T), N}(refs, pool)
end
convert(::Type{CategoricalArray{T, N}}, A::CategoricalArray{S, N, R}) where {S, T, N, R} =
    convert(CategoricalArray{T, N, R}, A)
convert(::Type{CategoricalArray{T}}, A::CategoricalArray{S, N, R}) where {S, T, N, R} =
    convert(CategoricalArray{T, N, R}, A)
convert(::Type{CategoricalArray}, A::CategoricalArray{T, N, R}) where {T, N, R} =
    convert(CategoricalArray{T, N, R}, A)

# R<:Integer is needed for this method to be considered more specific
# than the generic one above (JuliaLang/julia#18443)
convert(::Type{CategoricalArray{T, N, R}}, A::CategoricalArray{T, N, R}) where {T, N, R<:Integer} = A
convert(::Type{CategoricalArray{T, N}}, A::CategoricalArray{T, N}) where {T, N} = A
convert(::Type{CategoricalArray{T}}, A::CategoricalArray{T}) where {T} = A
convert(::Type{CategoricalArray}, A::CategoricalArray) = A

function Base.:(==)(A::CategoricalArray{S}, B::CategoricalArray{T}) where {S, T}
    if size(A) != size(B)
        return false
    end
    anymissing = false
    if A.pool === B.pool
        @inbounds for (a, b) in zip(A.refs, B.refs)
            if a == 0 || b == 0
                (S >: Missing || T >: Missing) && (anymissing = true)
            elseif a != b
                return false
            end
        end
    else
        @inbounds for (a, b) in zip(A, B)
            eq = (a == b)
            if eq === false
                return false
            elseif S >: Missing || T >: Missing
                anymissing |= ismissing(eq)
            end
        end
    end
    return anymissing ? missing : true
end

function Base.isequal(A::CategoricalArray, B::CategoricalArray)
    if size(A) != size(B)
        return false
    end
    if A.pool === B.pool
        @inbounds for (a, b) in zip(A.refs, B.refs)
            if a != b
                return false
            end
        end
    else
        @inbounds for (a, b) in zip(A, B)
            if !isequal(a, b)
                return false
            end
        end
    end
    return true
end

size(A::CategoricalArray) = size(A.refs)
Base.IndexStyle(::Type{<:CategoricalArray}) = IndexLinear()

@inline function setindex!(A::CategoricalArray, v::Any, I::Real...)
    @boundscheck checkbounds(A, I...)
    @inbounds A.refs[I...] = get!(A.pool, v)
end

Base.fill!(A::CategoricalArray, v::Any) =
    (fill!(A.refs, get!(A.pool, v)); A)

# Methods preserving levels and more efficient than AbstractArray fallbacks
copy(A::CategoricalArray{T, N}) where {T, N} =
    CategoricalArray{T, N}(copy(A.refs), copy(A.pool))

CatArrOrSub{T, N} = Union{CategoricalArray{T, N},
                          SubArray{<:Any, N, <:CategoricalArray{T}}} where {T, N}

function copyto!(dest::CatArrOrSub{T, N}, dstart::Integer,
                 src::CatArrOrSub{<:Any, N}, sstart::Integer,
                 n::Integer) where {T, N}
    n == 0 && return dest
    n < 0 && throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    (dstart ∈ destinds && dstart+n-1 ∈ destinds) || throw(BoundsError(dest, dstart:dstart+n-1))
    (sstart ∈ srcinds  && sstart+n-1 ∈ srcinds)  || throw(BoundsError(src,  sstart:sstart+n-1))

    drefs = refs(dest)
    srefs = refs(src)
    dpool = pool(dest)
    spool = pool(src)

    # try converting src to dest type to avoid partial copy corruption of dest
    # in the event that the src cannot be copied into dest
    slevs = convert(Vector{T}, levels(src))
    dlevs = levels(dest)
    if eltype(src) >: Missing && !(eltype(dest) >: Missing) && !all(x -> x > 0, srefs)
        throw(MissingException("cannot copy array with missing values to an array with element type $T"))
    end

    newlevels, ordered = mergelevels(isordered(dest), dlevs, slevs)
    if isordered(dest) && (length(newlevels) != length(dlevs))
        # Uncomment this when removing deprecation
        # throw(OrderedLevelsException(newlevels[findfirst(!in(Set(dlevs)), newlevels)],
        #                              dlevs))
        Base.depwarn("adding new levels to ordered CategoricalArray destination " *
                     "will throw an error in the future", :copyto!)
        ordered &= isordered(src) | (length(newlevels) == length(dlevs))
    end
    # Exception: empty pool marked as ordered if new value is ordered
    if isempty(dlevs) && isordered(src)
        ordered = true
    end
    if ordered != isordered(dest)
        isa(dest, SubArray) && throw(ArgumentError("cannot set ordered=$ordered on dest SubArray as it would affect the parent. Found when trying to set levels to $newlevels."))
        ordered!(dest, ordered)
    end

    # Simple case: replace all values
    if !isa(dest, SubArray) && dstart == dstart == 1 && n == length(dest) == length(src)
        # Set index to reflect refs
        levels!(dpool, T[]) # Needed in case src and dest share some levels
        levels!(dpool, index(spool))

        # Set final levels in their visible order
        levels!(dpool, newlevels)

        copyto!(drefs, srefs)
    else # More work to do: preserve some values (and therefore index)
        levels!(dpool, newlevels)

        indexmap = indexin(index(spool), index(dpool))

        @inbounds for i = 0:(n-1)
            x = srefs[sstart+i]
            drefs[dstart+i] = x > 0 ? indexmap[x] : 0
        end

    end

    dest
end

copyto!(dest::CatArrOrSub, src::CatArrOrSub) =
    copyto!(dest, 1, src, 1, length(src))

copyto!(dest::CatArrOrSub, dstart::Integer, src::CatArrOrSub) =
    copyto!(dest, dstart, src, 1, length(src))

if VERSION >= v"1.1"
    import Base: copy!
else
    import Future: copy!
end
copy!(dest::CatArrOrSub, src::CatArrOrSub) = copyto!(dest, 1, src, 1, length(src))

similar(A::CategoricalArray{S, M, R}, ::Type{T},
        dims::NTuple{N, Int}) where {T, N, S, M, R} =
    Array{T, N}(undef, dims)
similar(A::CategoricalArray{S, M, R}, ::Type{Missing},
        dims::NTuple{N, Int}) where {N, S, M, R} =
    Array{Missing, N}(missing, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{CategoricalValue{T, R}},
        dims::NTuple{N, Int}) where {R, T, N, S, M, Q} =
    CategoricalArray{T, N, R}(undef, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{CategoricalValue{T}},
        dims::NTuple{N, Int}) where {T, N, S, M, Q} =
    CategoricalArray{T, N, Q}(undef, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{Union{CategoricalValue{T, R}, Missing}},
        dims::NTuple{N, Int}) where {R, T, N, S, M, Q} =
    CategoricalArray{Union{T, Missing}, N, R}(undef, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{Union{CategoricalValue{T}, Missing}},
        dims::NTuple{N, Int}) where {T, N, S, M, Q} =
    CategoricalArray{Union{T, Missing}, N, Q}(undef, dims)

for A in (:Array, :Vector, :Matrix) # to fix ambiguities
    @eval begin
        similar(A::$A, ::Type{CategoricalValue{T, R}},
                dims::NTuple{N, Int}=size(A)) where {T, R, N} =
            CategoricalArray{T, N, R}(undef, dims)
        similar(A::$A, ::Type{CategoricalValue{T}},
                dims::NTuple{N, Int}=size(A)) where {T, N} =
            CategoricalArray{T, N}(undef, dims)
        similar(A::$A, ::Type{Union{CategoricalValue{T, R}, Missing}},
                dims::NTuple{N, Int}=size(A)) where {T, R, N} =
            CategoricalArray{Union{T, Missing}, N, R}(undef, dims)
        similar(A::$A, ::Type{Union{CategoricalValue{T}, Missing}},
                dims::NTuple{N, Int}=size(A)) where {T, N} =
            CategoricalArray{Union{T, Missing}, N}(undef, dims)
    end
end

similar(::Type{T}, dims::Dims) where {U, R, T<:Array{CategoricalValue{U, R}}} =
    CategoricalArray{eltype(T)}(undef, dims)
similar(::Type{T}, dims::Dims) where {U, T<:Array{CategoricalValue{U}}} =
    CategoricalArray{eltype(T)}(undef, dims)
similar(::Type{T}, dims::Dims) where {U, R, T<:Array{Union{CategoricalValue{U, R}, Missing}}} =
    CategoricalArray{eltype(T)}(undef, dims)
similar(::Type{T}, dims::Dims) where {U, T<:Array{Union{CategoricalValue{U}, Missing}}} =
    CategoricalArray{eltype(T)}(undef, dims)

"""
    compress(A::CategoricalArray)

Return a copy of categorical array `A` using the smallest reference type able to hold the
number of [`levels`](@ref) of `A`.

While this will reduce memory use, this function is type-unstable, which can affect
performance inside the function where the call is made. Therefore, use it with caution.
"""
function compress(A::CategoricalArray{T, N}) where {T, N}
    R = reftype(length(index(A.pool)))
    convert(CategoricalArray{T, N, R}, A)
end

"""
    decompress(A::CategoricalArray)

Return a copy of categorical array `A` using the default reference type ($DefaultRefType).
If `A` is using a small reference type (such as `UInt8` or `UInt16`) the decompressed array
will have room for more levels.

To avoid the need to call decompress, ensure [`compress`](@ref) is not called when creating
the categorical array.
"""
decompress(A::CategoricalArray{T, N}) where {T, N} =
    convert(CategoricalArray{T, N, DefaultRefType}, A)

function vcat(A::CategoricalArray...)
    ordered = any(isordered, A) && all(a->isordered(a) || isempty(levels(a)), A)
    newlevels, ordered = mergelevels(ordered, map(levels, A)...)

    refsvec = map(A) do a
        ii = convert(Vector{Int}, indexin(index(a.pool), newlevels))
        [x==0 ? 0 : ii[x] for x in a.refs]::Array{Int,ndims(a)}
    end

    T = Base.promote_eltype(A...) >: Missing ?
        Union{eltype(newlevels), Missing} : eltype(newlevels)
    refs = DefaultRefType[refsvec...;]
    pool = CategoricalPool(newlevels, ordered)
    CategoricalArray{T, ndims(refs)}(refs, pool)
end

@inline function getindex(A::CategoricalArray{T}, I...) where {T}
    @boundscheck checkbounds(A, I...)
    # Let Array indexing code handle everything
    @inbounds r = A.refs[I...]

    if isa(r, Array)
        res = CategoricalArray{T, ndims(r)}(r, copy(A.pool))
        return ordered!(res, isordered(A))
    else
        r > 0 || throw(UndefRefError())
        @inbounds res = A.pool[r]
        return res
    end
end

leveltype(::Type{T}) where {T <: CategoricalArray} = leveltype(nonmissingtype(eltype(T)))

"""
    levels(x::CategoricalArray)
    levels(x::CategoricalValue)

Return the levels of categorical array or value `x`.
This may include levels which do not actually appear in the data
(see [`droplevels!`](@ref)).
"""
DataAPI.levels(A::CategoricalArray) = levels(A.pool)

"""
    levels!(A::CategoricalArray, newlevels::Vector; allow_missing::Bool=false)

Set the levels categorical array `A`. The order of appearance of levels will be respected
by [`levels`](@ref), which may affect display of results in some operations; if `A` is
ordered (see [`isordered`](@ref)), it will also be used for order comparisons
using `<`, `>` and similar operators. Reordering levels will never affect the values
of entries in the array.

If `A` accepts missing values (i.e. `eltype(A) >: Missing`) and `allow_missing=true`,
entries corresponding to omitted levels will be set to `missing`.
Else, `newlevels` must include all levels which appear in the data.
"""
function levels!(A::CategoricalArray{T}, newlevels::Vector; allow_missing=false) where {T}
    if !allunique(newlevels)
        throw(ArgumentError(string("duplicated levels found: ",
                                   join(unique(filter(x->sum(newlevels.==x)>1, newlevels)), ", "))))
    end

    # first pass to check whether, if some levels are removed, changes can be applied without error
    # TODO: save original levels and undo changes in case of error to skip this step
    # equivalent to issubset but faster due to JuliaLang/julia#24624
    if !isempty(setdiff(index(A.pool), newlevels))
        deleted = [!(l in newlevels) for l in index(A.pool)]
        @inbounds for (i, x) in enumerate(A.refs)
            if T >: Missing
                !allow_missing && x > 0 && deleted[x] &&
                    throw(ArgumentError("cannot remove level $(repr(index(A.pool)[x])) as it is used at position $i and allow_missing=false."))
            else
                deleted[x] &&
                    throw(ArgumentError("cannot remove level $(repr(index(A.pool)[x])) as it is used at position $i. " *
                                        "Change the array element type to Union{$T, Missing} using convert if you want to transform some levels to missing values."))
            end
        end
    end

    # actually apply changes
    oldindex = copy(index(A.pool))
    levels!(A.pool, newlevels)

    if index(A.pool) != oldindex
        levelsmap = similar(A.refs, length(oldindex)+1)
        # 0 maps to a missing value
        levelsmap[1] = 0
        levelsmap[2:end] .= something.(indexin(oldindex, index(A.pool)), 0)

        @inbounds for (i, x) in enumerate(A.refs)
            A.refs[i] = levelsmap[x+1]
        end
    end

    A
end

function _unique(::Type{S},
                 refs::AbstractArray{T},
                 pool::CategoricalPool) where {S, T<:Integer}
    nlevels = length(index(pool)) + 1
    order = fill(0, nlevels) # 0 indicates not seen
    # If we don't track missings, short-circuit even if none has been seen
    count = S >: Missing ? 0 : 1
    @inbounds for i in refs
        if order[i + 1] == 0
            count += 1
            order[i + 1] = count
            count == nlevels && break
        end
    end
    S[i == 1 ? missing : index(pool)[i - 1] for i in sortperm(order) if order[i] != 0]
end

"""
    unique(A::CategoricalArray)

Return levels which appear in `A` in their order of appearance.
This function is significantly slower than [`levels`](@ref)
since it needs to check whether levels are used or not.
"""
unique(A::CategoricalArray{T}) where {T} = _unique(T, A.refs, A.pool)

if VERSION >= v"0.7.0-DEV.4882"
    """
        droplevels!(A::CategoricalArray)

    Drop levels which do not appear in categorical array `A` (so that they will no longer be
    returned by [`levels`](@ref)).
    """
    droplevels!(A::CategoricalArray) = levels!(A, intersect!(levels(A), unique(A)))
else # intersect! method missing on Julia 0.6
    """
        droplevels!(A::CategoricalArray)

    Drop levels which do not appear in categorical array `A` (so that they will no longer be
    returned by [`levels`](@ref)).
    """
    droplevels!(A::CategoricalArray) = levels!(A, intersect(levels(A), filter!(!ismissing, unique(A))))
end

"""
    isordered(A::CategoricalArray)

Test whether entries in `A` can be compared using `<`, `>` and similar operators,
using the ordering of levels.
"""
isordered(A::CategoricalArray) = isordered(A.pool)

"""
    ordered!(A::CategoricalArray, ordered::Bool)

Set whether entries in `A` can be compared using `<`, `>` and similar operators,
using the ordering of levels. Return the modified `A`.
"""
ordered!(A::CategoricalArray, ordered) = (ordered!(A.pool, ordered); return A)

function Base.resize!(A::CategoricalVector, n::Integer)
    n_orig = length(A)
    resize!(A.refs, n)
    if n > n_orig
        A.refs[n_orig+1:end] .= 0
    end
    A
end

function Base.push!(A::CategoricalVector, item)
    r = get!(A.pool, item)
    push!(A.refs, r)
    A
end

function Base.append!(A::CategoricalVector, B::CatArrOrSub)
    levels!(A, union(levels(A), levels(B)))
    len = length(A)
    len2 = length(B)
    resize!(A.refs, len + len2)
    for i = 1:len2
        A[len + i] = B[i]
    end
    return A
end

Base.empty!(A::CategoricalArray) = (empty!(A.refs); return A)

function Base.reshape(A::CategoricalArray{T, N}, dims::Dims) where {T, N}
    x = reshape(A.refs, dims)
    res = CategoricalArray{T, ndims(x)}(x, A.pool)
    ordered!(res, isordered(res))
end

"""
    categorical{T}(A::AbstractArray{T}[, compress::Bool]; ordered::Bool=false)

Construct a categorical array with the values from `A`.

If the element type supports it, levels are sorted in ascending order;
else, they are kept in their order of appearance in `A`. The `ordered` keyword
argument determines whether the array values can be compared according to the
ordering of levels or not (see [`isordered`](@ref)).

If `compress` is provided and set to `true`, the smallest reference type able to hold the
number of unique values in `A` will be used. While this will reduce memory use, passing
this parameter will also introduce a type instability which can affect performance inside
the function where the call is made. Therefore, use this option with caution (the
one-argument version does not suffer from this problem).

    categorical{T}(A::CategoricalArray{T}[, compress::Bool]; ordered::Bool=isordered(A))

If `A` is already a `CategoricalArray`, its levels are preserved;
the same applies to the ordered property, and to the reference type
unless `compress` is passed.
"""
function categorical end

categorical(A::AbstractArray; ordered=_isordered(A)) = CategoricalArray(A, ordered=ordered)

# Type-unstable methods
function categorical(A::AbstractArray{T, N}, compress; ordered=_isordered(A)) where {T, N}
    RefType = compress ? reftype(length(unique(A))) : DefaultRefType
    CategoricalArray{T, N, RefType}(A, ordered=ordered)
end
function categorical(A::CategoricalArray{T, N, R}, compress; ordered=_isordered(A)) where {T, N, R}
    RefType = compress ? reftype(length(levels(A))) : R
    CategoricalArray{T, N, RefType}(A, ordered=ordered)
end

function in(x::Any, y::CategoricalArray{T, N, R}) where {T, N, R}
    ref = get(y.pool, x, zero(R))
    ref != 0 ? ref in y.refs : false
end

function in(x::CategoricalValue, y::CategoricalArray{T, N, R}) where {T, N, R}
    if x.pool === y.pool
        return x.level in y.refs
    else
        ref = get(y.pool, index(x.pool)[x.level], zero(R))
        return ref != 0 ? ref in y.refs : false
    end
end

Array(A::CategoricalArray{T}) where {T} = Array{T}(A)
collect(A::CategoricalArray{T}) where {T} = Array{T}(A)

function float(A::CategoricalArray{T}) where T
    if !isconcretetype(T)
        error("`float` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(float(zero(T)))}, A)
end
function complex(A::CategoricalArray{T}) where T
    if !isconcretetype(T)
        error("`complex` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(complex(zero(T)))}, A)
end

# Override AbstractArray method to avoid printing useless type parameters
if VERSION >= v"0.7.0-DEV.2657"
    summary(io::IO, A::CategoricalArray{T, N, R}) where {T, N, R} =
        print(io, Base.dims2string(size(A)), " $CategoricalArray{$T,$N,$R}")
else
    summary(A::CategoricalArray{T, N, R}) where {T, N, R} =
        string(Base.dims2string(size(A)), " $CategoricalArray{$T,$N,$R}")
end

refs(A::CategoricalArray) = A.refs
pool(A::CategoricalArray) = A.pool

Base.deleteat!(A::CategoricalArray, inds) = (deleteat!(A.refs, inds); A)

Base.Broadcast.broadcasted(::typeof(ismissing), A::CategoricalArray{T}) where {T} =
    T >: Missing ? Base.Broadcast.broadcasted(==, A.refs, 0) :
                   Base.Broadcast.broadcasted(_ -> false, A.refs)

Base.Broadcast.broadcasted(::typeof(!ismissing), A::CategoricalArray{T}) where {T} =
    T >: Missing ? Base.Broadcast.broadcasted(>, A.refs, 0) :
                   Base.Broadcast.broadcasted(_ -> true, A.refs)

function Base.Broadcast.broadcasted(::typeof(levelcode), A::CategoricalArray{T}) where {T}
    ord = order(A.pool)
    if T >: Missing
        Base.Broadcast.broadcasted(i -> i > 0 ? Signed(widen(ord[i])) : missing, A.refs)
    else
        Base.Broadcast.broadcasted(i -> Signed(widen(ord[i])), A.refs)
    end
end

function Base.sort!(v::CategoricalVector;
                    # alg is ignored since counting sort is more efficient
                    alg::Base.Algorithm=Base.Sort.defalg(v),
                    lt=isless,
                    by=identity,
                    rev::Bool=false,
                    order::Base.Ordering=Base.Forward)
    counts = zeros(UInt, length(v.pool) + (eltype(v) >: Missing))

    # do a count/histogram of the references
    @inbounds for ref in v.refs
        counts[ref + (eltype(v) >: Missing)] += 1
    end

    # compute the order in which to read from counts
    ord = Base.Sort.ord(lt, by, rev, order)
    index = eltype(v) >: Missing ? [missing; v.pool.valindex] : v.pool.valindex
    seen = counts .> 0
    anymissing = eltype(v) >: Missing && seen[1]
    perm = sortperm(view(index, seen), order=ord)
    nzcounts = counts[seen]
    j = 0
    @inbounds for ref in perm
        tmpj = j + nzcounts[ref]
        v.refs[(j+1):tmpj] .= ref - anymissing
        j = tmpj
    end

    return v
end


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

mutable struct AuthenticationMethod <: SwaggerModel

    function AuthenticationMethod(;)
        o = new()
        o
    end
end # type AuthenticationMethod

const _property_map_AuthenticationMethod = Dict{Symbol,Symbol}()
const _property_types_AuthenticationMethod = Dict{Symbol,String}()
Base.propertynames(::Type{ AuthenticationMethod }) = collect(keys(_property_map_AuthenticationMethod))
Swagger.property_type(::Type{ AuthenticationMethod }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AuthenticationMethod[name]))}
Swagger.field_name(::Type{ AuthenticationMethod }, property_name::Symbol) =  _property_map_AuthenticationMethod[property_name]

function check_required(o::AuthenticationMethod)
    true
end

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

ext: jl


text:
struct DistanceForest
   forests::Vector{PyCall.PyObject}
   intervals::Vector{Tuple{Int,Int}}
   midpoints::Vector{Int}
   midvalues::Vector{Vector{Float64}}
   heuristic::Vector{Function}

   function DistanceForest( outratio::Float64=0.1, heuristic=[>, <, >, >, >] )
      inter, mids = distances()
      forests = [IsolationForest(contamination=outratio) for i in 1:length(mids)]
      new(forests, inter, mids, Vector{Vector{Float64}}(length(mids)), heuristic)
   end
end

function train!( df::DistanceForest, array::Array{Float64,2}, dist_col::Int=1 )
   data = sortrows( array )
   training = data[:, 2:end]
   dist = data[:, 1]
   for i in 1:length(df.intervals)
      curint = df.intervals[i]
      start = searchsortedfirst(dist, curint[1])
      stop  = searchsortedlast(dist, curint[2])
      df.midvalues[i] = Float64[]
      for j in 1:size(training, 2)
         push!( df.midvalues[i], mean(training[start:stop, j]) )
      end
      if length(start:stop) > 1
         fit!( df.forests[i], training[start:stop, :] )
      end
   end
end

predict( df::DistanceForest, array::Vector{Float64} ) = predict( df, array' )

function predict( df::DistanceForest, array::Array{Float64,2} )
   dist = array[:, 1]
   assigned = map(x->searchsortedfirst( df.midpoints, x ), dist)
   outliers = Int[]

   function validate( row::Int, meanvalues::Vector{Float64} )
      for k in 2:size(array, 2)
         if !(df.heuristic[k-1]( array[row, k], meanvalues[k-1] ))
            return false
         end
      end
      return true
   end

   for i in 1:length(df.midpoints)
      range = searchsorted( assigned, i )
      if length(range) >= 1
         res = ScikitLearn.predict( df.forests[i], array[range, 2:end] )
         for j in 1:length(res)
            if res[j] == -1 && validate( j + first(range) - 1, df.midvalues[i] )
               push!( outliers, j + first(range) - 1 )
            end
         end
      end
   end
   outliers
end

#function distances( b::Int=100, range::UnitRange=1:1 )
function distances( b::Int=2, range::UnitRange=2:12 )
   first = [(b^i, b^(i+1)) for i in range]
   last  = [(b^i+b^(i-1), b^(i+1)+b^i) for i in range]
   first[1] = (1, first[1][2])
   last[1]  = (1, last[1][2])
   ret = Vector{Tuple{Int,Int}}()
   retmean = Int[]
   midpoint( left, right ) = ((right - left) >> 1) + left
   for i in 1:length(first)
      push!(ret, first[i])
      push!(retmean, midpoint(first[i][1], first[i][2]))
      push!(ret, last[i])
      push!(retmean, midpoint(last[i][1], last[i][2]))
   end
   ret, retmean
end


text:
"""
```
    Matching <: AbstractHistogramAdjustmentAlgorithm
    Matching(targetimg; nbins = 256, edges = nothing)

    adjust_histogram([T,] img, f::Matching)
    adjust_histogram!([out,] img, f::Matching)
```

Returns a histogram matched image with a granularity of `nbins` number of bins.
The first argument `img` is the image to be matched, whereas the argument
`targetimg` in `Matching()` is the image having the desired histogram to be
matched to.

# Details
The purpose of histogram matching is to transform the intensities in a source
image so that the intensities distribute according to the histogram of a
specified target image. If one interprets histograms as piecewise-constant
models of probability density functions (see [`build_histogram`](@ref
build_histogram(::AbstractArray, ::Integer, ::Union{Real,AbstractGray},
::Union{Real,AbstractGray}))), then the histogram matching task can be modelled
as the problem of transforming one probability distribution into another [1].
It turns out that the solution to this transformation problem involves the
cumulative and inverse cumulative distribution functions of the source and
target probability density functions.

In particular, let the random variables ``x \\thicksim p_{x} `` and ``z
\\thicksim p_{z}``  represent an intensity in the source and target image
respectively, and let

```math
 S(x) = \\int_0^{x}p_{x}(w)\\mathrm{d} w \\quad \\text{and} \\quad
 T(z) = \\int_0^{z}p_{z}(w)\\mathrm{d} w
```
represent their concomitant cumulative distribution functions. Then the
sought-after mapping ``Q(\\cdot)`` such that ``Q(x) \\thicksim p_{z} `` is given
by

```math
Q(x) =  T^{-1}\\left( S(x) \\right),
```

where ``T^{-1}(y) = \\operatorname{min} \\{ x \\in \\mathbb{R} : y \\leq T(x)
\\}`` is the inverse cumulative distribution function of ``T(x)``.

The mapping suggests that one can conceptualize histogram matching as performing
histogram equalization on the source and target image and relating the two
equalized histograms. Refer to [`adjust_histogram`](@ref
adjust_histogram(::Equalization, ::AbstractArray, ::Integer,
::Union{Real,AbstractGray}, ::Union{Real,AbstractGray})) for more details on
histogram equalization.

# Options

Various options for the parameters of the `adjust_histogram` function and `Matching`
type are described in more detail below.

## Choices for `img` and `targetimg`

The `adjust_histogram(img, Matching())` function can handle a variety of input
types. The type of the returned image matches the input type.

For colored images, the inputs are converted to
[YIQ](https://en.wikipedia.org/wiki/YIQ)  type and the distributions of the Y
channels are matched. The modified Y channel is then combined with the I and Q
channels and the resulting image converted to the same type as the input.

## Choices for `nbins`

You can specify the total number of bins in the histogram. If you do not
specify the number of bins then a default value of 256 bins is utilized.

## Choices for `edges`

If you do not designate the number of bins, then you have the option to directly
stipulate how the intervals will be divided by specifying a `AbstractRange`
type.

# Example

```julia
using Images, TestImages, ImageView

img_source = testimage("mandril_gray")
img_target = adjust_gamma(img_source, 1/2)
img_transformed = adjust_histogram(img_source, Matching(targetimg = img_target))
#=
    A visual inspection confirms that img_transformed resembles img_target
    much more closely than img_source.
=#
imshow(img_source)
imshow(img_target)
imshow(img_transformed)
```

# References
1. W. Burger and M. J. Burge. *Digital Image Processing*. Texts in Computer Science, 2016. [doi:10.1007/978-1-4471-6684-9](https://doi.org/10.1007/978-1-4471-6684-9)
"""
Base.@kwdef struct Matching{T₁ <: AbstractArray,
                            T₂ <: Union{Integer, Nothing},
                            T₃ <: Union{AbstractRange, Nothing}} <: AbstractHistogramAdjustmentAlgorithm
    targetimg::T₁
    nbins::T₂ = 256
    edges::T₃ = nothing
end

function (f::Matching)(out::GenericGrayImage, img::GenericGrayImage)
    #TODO Throw error/warning if user specifies both edges and nbins simultaneously.
    out .= img
    edges, pdf, target_pdf = isnothing(f.edges) ? construct_pdfs(out, f.targetimg, f.nbins) : construct_pdfs(out, f.targetimg, f.edges)
    match_pdf!(out, edges, pdf, target_pdf)
    return out
end

function (f::Matching)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    T = eltype(img)
    yiq = convert.(YIQ, img)
    yiq_view = channelview(yiq)

    targetimg = f.targetimg
    yiq_targetimg = convert.(YIQ, targetimg)
    yiq_targetimg_view = channelview(yiq_targetimg)

    adjust_histogram!(view(yiq_view,1,:,:), Matching(targetimg = view(yiq_targetimg_view,1,:,:),
                                                         edges = f.edges,
                                                         nbins = f.nbins))
    out .= convert.(T, yiq)
end

(f::Matching)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))

function match_pdf!(img::AbstractArray, edges::AbstractArray, pdf::AbstractArray, target_pdf::AbstractArray)
    cdf = parent(cumsum(pdf))
    target_cdf = parent(cumsum(target_pdf))
    # Precompute the inverse cummulative distribution function of target_cdf.
    lookup_table = lookup_icdf(cdf, target_cdf)
    # Transform the intensities in img so that they are distributed according
    # to the distribution of the target_histogram.
    T = eltype(img)
    step_size = step(edges)
    first_edge = first(edges)
    last_edge = last(edges)
    map!(img, img) do val
        if isnan(val)
            return val
        else
            if val >= last_edge
                newval = edges[last(lookup_table)-1] + step_size
            elseif val < first_edge
                newval = edges[first(lookup_table)]
            else
                index = Int(Base.div(val-first_edge,step_size)) + 1
                newval = edges[lookup_table[index]]
            end
            return T <: Integer ? ceil(newval) : newval
        end
    end
end


text:
###############################################
# Solve laplace equation

# by: Jose Miguel Muñoz
# please contact me if this doesn't works
###############################################

using Plots, Statistics
###############################################
#       BOUNDARY
###############################################
N = 100
a , b = 0, 2

Rx = a:(b-a)/(N-1):b
Ry = a:(b-a)/(N-1):b

boundary_X_Up = map(rx -> 0,Rx)

boundary_X_Down = map(rx -> 0, Rx)

boundary_Y_R = map(rx -> 10,Rx)

boundary_Y_L = map(rx -> 0,Rx)

V = zeros(N, N)

V[1,:] = boundary_X_Up
V[:,1] = boundary_Y_R
V[N,:] = boundary_X_Down
V[:,N] = boundary_Y_L

###############################################
#       SOLVER
###############################################
function solver!(V,n_inter,N=N, tol=1)
    for rep = 1:n_inter
        nV =  (V)
        for row = 2:N-1
            for column = 2:N-1
                V[row, column] = mean([nV[row-1, column],nV[row+1, column],nV[row, column-1],nV[row, column+1]])
                #if (row+10)^2+(column+10)^2<1000
                #    V[row, column]=1
                #end
    end end
    V =  (nV)
 end
    return V
end

###############################################
#       Ask for the boundaries
###############################################
V = solver!(V,100000)
@benchmark solver!(V,1000)
#contour(V,levels=100)
plot(Ry, Ry, V, linetype=:contourf,title = "Solucion por Jacobi")

#gr()
#anim = @animate for n = 1:100:1000
#    heatmap(vss[:,:,n])
#end
#gif(anim)
##############################
#   Errores
##############################
using Statistics
plot(Ry, Ry, (V-Space*1.01), linetype=:contourf,title = "Solucion por Jacobi")

std(V.^2-Space.^2)

#######################################################
#       Convergence Analysis
#######################################################
listconvJacobi=[]
for n_reps=500:100:5000
    V = rand(N, N)

    V[1,:] = boundary_X_Up
    V[:,1] = boundary_Y_R
    V[N,:] = boundary_X_Down
    V[:,N] = boundary_Y_L

    V = solver!(V,n_reps)
    append!(listconvJacobi,abs(mean(V-Space)))
end
listconvJacobi

max_stars_repo_name: munozariasjm/Laplace-and-Julia

max_issues_repo_name: munozariasjm/Laplace-and-Julia

max_forks_repo_name: munozariasjm/Laplace-and-Julia


text:
# constraint Jacobian
function constraint_jacobian_configuration(mechanism, contact::ContactConstraint, body::Body)
    relative = (body.id == contact.parent_id ? :parent : :child)

    timestep = mechanism.timestep
    pbody = get_body(mechanism, contact.parent_id)
    cbody = get_body(mechanism, contact.child_id)

    return constraint_jacobian_configuration(relative,
        contact.model,
        next_configuration_velocity(pbody.state, timestep)...,
        next_configuration_velocity(cbody.state, timestep)...,
        mechanism.timestep)
end

function constraint_jacobian_velocity(mechanism, contact::ContactConstraint, body::Body)
    relative = (body.id == contact.parent_id ? :parent : :child)

    timestep = mechanism.timestep
    pbody = get_body(mechanism, contact.parent_id)
    cbody = get_body(mechanism, contact.child_id)

    return constraint_jacobian_velocity(relative,
        contact.model,
        next_configuration_velocity(pbody.state, timestep)...,
        next_configuration_velocity(cbody.state, timestep)...,
        mechanism.timestep)
end

# impulses
function impulses!(mechanism, body::Body, contact::ContactConstraint)
    body.state.d -= impulse_map(mechanism, contact, body) * contact.impulses[2]
    return
end

function impulse_map(mechanism, contact::ContactConstraint, body::Body)
    relative = (body.id == contact.parent_id ? :parent : :child)
    pbody = get_body(mechanism, contact.parent_id)
    cbody = get_body(mechanism, contact.child_id)
    return impulse_map(relative, contact.model, pbody, cbody, mechanism.timestep)
end

function impulse_map_jacobian_configuration(mechanism, body::Body{T}, contact::ContactConstraint{T}) where T
    relative = (body.id == contact.parent_id ? :parent : :child)

    return impulse_map_jacobian(relative, relative, contact.model,
        get_body(mechanism, contact.parent_id),
        get_body(mechanism, contact.child_id),
        contact.impulses[2],
        mechanism.timestep)
end

function impulses_jacobian_velocity!(mechanism, body::Body, contact::ContactConstraint)
    timestep = mechanism.timestep
    body.state.D -= impulse_map_jacobian_configuration(mechanism, body, contact) * integrator_jacobian_velocity(body, timestep)
    return
end

# off-diagonal terms for linear system
function off_diagonal_jacobians(mechanism, body::Body, contact::ContactConstraint{T,N,Nc,Cs,N½}) where {T,N,Nc,Cs,N½}
    Z = szeros(T,N½,6)
    return [Z' -impulse_map(mechanism, contact, body)], [Z; constraint_jacobian_velocity(mechanism, contact, body)]
end

# off-diagonal terms for linear system
function off_diagonal_jacobians(mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½}, body::Body) where {T,N,Nc,Cs,N½}
    Z = szeros(T,N½,6)
    return [Z; constraint_jacobian_velocity(mechanism, contact, body)], [Z' -impulse_map(mechanism, contact, body)]
end

# linear system entries
function set_matrix_vector_entries!(mechanism::Mechanism, matrix_entry::Entry, vector_entry::Entry, contact::ContactConstraint)
    matrix_entry.value = constraint_jacobian(contact)
    vector_entry.value = vcat(-complementarityμ(mechanism, contact), -constraint(mechanism, contact))
    return
end

# reset variables using cone-specific neutral vector
function reset!(contact::ContactConstraint;
    scale=1.0)
    contact.impulses[1] = scale * neutral_vector(contact.model)
    contact.impulses[2] = scale * neutral_vector(contact.model)
    contact.impulses_dual[1] = scale * neutral_vector(contact.model)
    contact.impulses_dual[2] = scale * neutral_vector(contact.model)
    return
end

# dimension of cone
cone_degree(contact::ContactConstraint) = sum(cone_degree(contact.model))


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

name = "libinchi"
version = v"1.5.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mojaie/libinchi.git", "3679104c28bd404ecb9faf93cf2be3ac883c13c0")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libinchi
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
if [[ "${target}" == *-mingw* ]]; then
    mkdir -p "${libdir}"
    cp "libinchi.${dlext}" "${libdir}"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "windows"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "linux"; libc="glibc")
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libinchi", :libinchi)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module FieldsTests

using Test

@testset "Fields" begin include("FieldInterfacesTests.jl") end

#@testset "FieldOperations" begin include("FieldOperationsTests.jl") end

@testset "FieldArrays" begin include("FieldArraysTests.jl") end

@testset "DiffOperators" begin include("DiffOperatorsTests.jl") end

@testset "AffineMaps" begin include("AffineMapsTests.jl") end

@testset "FieldArraysOperations" begin include("FieldArraysOperationsTests.jl") end

#@testset "LazyFieldArrays" begin include("LazyFieldArraysTests.jl") end

@testset "BlockFieldArrays" begin include("BlockFieldArraysTests.jl") end

end


text:
using Pkg
cd(@__DIR__)
Pkg.activate("../"); Pkg.instantiate()

using Plots; gr()
using LaTeXStrings
using DataFrames
using CSV
using OrdinaryDiffEq
using LinearAlgebra
using DataInterpolations
using BenchmarkTools
using DiffEqFlux
using Optim
using Flux
using Zygote
using ModelingToolkit
using DataDrivenDiffEq

model_univ = FastChain(FastDense(6, 16, tanh),
                       FastDense(16, 2))
include("../optimization_results/neural_augmentation_results_u.jl")
include("../optimization_results/globaloptimizationresults.jl")

# Plot path
figpath = "../figs/";
# Linewidths and styles
LW1 = 2.5
LW2 = 1.5
LS1 = :solid
LS2 = :dot
LS3 = :dash
LC1 = :blue
LC2 = :red
LC3 = :green
LC4 = :black
LA1 = 1.0
LA2 = 0.7
LA3 = 0.3
;

df=CSV.read("../data/Lab-Helicopter_Experimental-data.csv")
ENV["COLUMNS"]=100 # change the number of IJulia columns from default 80 to 100
df[1:5,:]

u_th = df[:,1]
u_psi = df[:,2]
theta = df[:,3]
psi = df[:,4]
dt = df[1,5]
tm = 0:dt:(length(psi)-1)*dt
;

timeseriesdata = [theta';psi']
Nd = length(tm);

U_th = ConstantInterpolation(u_th,tm)
U_psi = ConstantInterpolation(u_psi,tm)
Theta = ConstantInterpolation(theta,tm)
Psi = ConstantInterpolation(psi,tm);
u_t = (th=U_th, psi=U_psi)
rad2deg2(t,u) = (t,u*180/pi)

#### Sparsify

# Create a Basis
@variables u[1:2]

# Array of potential basis functions
h = [cos.(u)...; sin.(u)...; u[1]^0; u[1]^2; u[2]^2; u[1]*u[2]]
basis_u = Basis(h, u)
println(basis_u)
# Create an optimizer for the SINDY problem
opt = SR3()
# Create the thresholds which should be used in the search process
λ = exp10.(-10:0.1:10)
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
f_target(x, w) = iszero(x[1]) ? Inf : norm(w.*x, 2)

# Test on original data and without further knowledge
X = reduce(hcat,utimeseries)
DX = reduce(hcat,Ktimeseries)
Ψ = SInDy(X, DX, basis_u, λ, opt = opt, maxiter = 10000, f_target = f_target)
println(Ψ)
print_equations(Ψ)

###########
## Check the fit
###########

function helicopter_uode(dx,x,Ps,P0,u,p,t)
# P0 = (rc=rc,h=h,m=m,Jbc=Jbc,Js=Js,g=g,K_th_th=K_th_th,K_th_psi=K_th_psi,K_psi_th=K_psi_th,K_psi_psi=K_psi_psi,D_th=D_th,D_psi=D_psi);
    th = x[1]
    psi = x[2]
    w_th = x[3]
    w_psi = x[4]

    U = [u.th(t),u.psi(t)]
    K = Ψ.equations.f_(U,Ψ.parameters,t)

    #
    M11 = Ps[4]*P0.Jbc + Ps[3]*P0.m*((Ps[1]*P0.rc)^2+(Ps[2]*P0.h)^2)
    M22 = Ps[4]*P0.Jbc*cos(th)^2 + Ps[3]*P0.m*(Ps[1]*P0.rc*cos(th)-Ps[2]*P0.h*sin(th))^2 + Ps[5]*P0.Js
    #
    dLdth = -w_psi^2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th)-(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2)) - (Ps[3]*P0.m)*(Ps[6]*P0.g)*((Ps[1]*P0.rc)*cos(th) - (Ps[2]*P0.h)*sin(th))
    #
    dp_psi_dth = -2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th) + (Ps[3]*P0.m)*(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2))*w_psi
    #
    F_th = Ps[7]*P0.K_th_th*U[1] - Ps[8]*P0.K_th_psi*U[2] - P0.D_th*w_th + K[1]
    F_psi = Ps[9]*P0.K_psi_th*U[1] - Ps[10]*P0.K_psi_psi*U[2] - P0.D_psi*w_psi + K[2]

    #
    dx[1] = w_th
    dx[2] = w_psi
    dx[3] = (dLdth + F_th)/M11
    dx[4] = (-dp_psi_dth*w_th + F_psi)/M22
end

uode_f! = (dx,x,p,t) -> helicopter_uode(dx,x,Ps_best,P0,u_t,p,t)
prob = ODEProblem(uode_f!,x0_best,tspan)
sol_best = solve(prob,AutoTsit5(TRBDF2(autodiff=false)),reltol=1e-6)

model_univ_u = FastChain(FastDense(2, 16, tanh),
                       FastDense(16, 2))

function helicopter_uode(dx,x,Ps,P0,u,p,t)
# P0 = (rc=rc,h=h,m=m,Jbc=Jbc,Js=Js,g=g,K_th_th=K_th_th,K_th_psi=K_th_psi,K_psi_th=K_psi_th,K_psi_psi=K_psi_psi,D_th=D_th,D_psi=D_psi);
    th = x[1]
    psi = x[2]
    w_th = x[3]
    w_psi = x[4]

    U = [u.th(t),u.psi(t)]
    K = model_univ_u(U,p)
    #
    M11 = Ps[4]*P0.Jbc + Ps[3]*P0.m*((Ps[1]*P0.rc)^2+(Ps[2]*P0.h)^2)
    M22 = Ps[4]*P0.Jbc*cos(th)^2 + Ps[3]*P0.m*(Ps[1]*P0.rc*cos(th)-Ps[2]*P0.h*sin(th))^2 + Ps[5]*P0.Js
    #
    dLdth = -w_psi^2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th)-(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2)) - (Ps[3]*P0.m)*(Ps[6]*P0.g)*((Ps[1]*P0.rc)*cos(th) - (Ps[2]*P0.h)*sin(th))
    #
    dp_psi_dth = -2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th) + (Ps[3]*P0.m)*(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2))*w_psi
    #
    F_th = Ps[7]*P0.K_th_th*U[1] - Ps[8]*P0.K_th_psi*U[2] - P0.D_th*w_th + K[1]
    F_psi = Ps[9]*P0.K_psi_th*U[1] - Ps[10]*P0.K_psi_psi*U[2] - P0.D_psi*w_psi + K[2]

    #
    dx[1] = w_th
    dx[2] = w_psi
    dx[3] = (dLdth + F_th)/M11
    dx[4] = (-dp_psi_dth*w_th + F_psi)/M22
end

uode_f! = (dx,x,p,t) -> helicopter_uode(dx,x,Ps_best,P0,u_t,p,t)

tspan = (tm[1],tm[end])

#
prob = ODEProblem(uode_f!,x0_best,tspan,pnn)
#
sol_nn = solve(prob,AutoTsit5(TRBDF2(autodiff=false)),reltol=1e-6)

#
plot(sol_best,vars=(rad2deg2,0,1),lw=LW1,lc=LC1,label=L"$\theta$")
plot!(sol_nn,vars=(rad2deg2,0,1),lw=LW1,lc=LC3,label=L"$\theta$")
plot!(tm,theta*180/pi,lw=LW1,ls=LS2,lc=LC2,label=L"\theta^\mathrm{d}")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel=L"$\theta$ [$^\circ$]")
plot!(title="Pitch angle: discovered (blue) vs. data (red) vs NN (green)",box=true)
figname="Helicopter_pitch-angle_model_fit_FNN_u_eq.svg"
savefig(figpath*figname)

plot(sol_best,vars=(rad2deg2,0,2),lw=LW1,lc=LC1,label=L"$\psi$")
plot!(sol_nn,vars=(rad2deg2,0,2),lw=LW1,lc=LC3,label=L"$\psi$")
plot!(tm,psi*180/pi,lw=LW1,ls=LS2,lc=LC2,label=L"\psi^\mathrm{d}")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel=L"$\psi$ [$^\circ$]")
plot!(title="Yaw angle: discovered(blue) vs. data (red) vs NN (green)",box=true)
figname="Helicopter_yaw-angle_model_fit_FNN_u_eq.svg"
savefig(figpath*figname)

open("../optimization_results/equation_discovery_u.txt", "w") do f
    println(f,Ψ)
    println(f)
    print_equations(f,Ψ)
    println(f)
    print_equations(f,Ψ,show_parameter=true)
    println(f)
    write(f,"parameters = $(Ψ.parameters) \n")
end


text:
export policy_evaluation!, policy_improvement!, policy_iteration!, value_iteration!

"""
    policy_evaluation!(V::AbstractApproximator, π, model::AbstractEnvironmentModel; γ::Float64=0.9, θ::Float64=1e-4)

See more details at Section (4.1) on Page 75 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
function policy_evaluation!(
    ;
    V::AbstractApproximator,
    π::AbstractPolicy,
    model::AbstractEnvironmentModel,
    γ::Float64 = 0.9,
    θ::Float64 = 1e-4,
)
    states, actions = 1:length(get_observation_space(model)), 1:length(get_action_space(model))
    while true
        Δ = 0.0
        for s in states
            v = sum(
                a -> get_prob(π, s, a) *
                     sum(p * (r + γ * V(s′)) for (s′, r, p) in model(s, a)),
                actions,
            )
            error = v - V(s)
            update!(V, s => error)
            Δ = max(Δ, abs(error))
        end
        Δ < θ && break
    end
    V
end

"""
    policy_improvement!(;V::AbstractApproximator, π::AbstractPolicy, model::AbstractEnvironmentModel, γ::Float64 = 0.9)

See more details at Section (4.2) on Page 76 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
function policy_improvement!(
    ;
    V::AbstractApproximator,
    π::AbstractPolicy,
    model::AbstractEnvironmentModel,
    γ::Float64 = 0.9,
)
    states, actions = 1:length(get_observation_space(model)), 1:length(get_action_space(model))
    is_policy_stable = true
    for s in states
        old_a = π(s)
        best_action_inds = find_all_max([sum(p * (r + γ * V(s′)) for (s′, r, p) in model(
            s,
            a,
        )) for a in actions])[2]
        new_a = actions[sample(best_action_inds)]
        if new_a != old_a
            update!(π, s => new_a)
            is_policy_stable = false
        end
    end
    is_policy_stable
end

"""
    policy_iteration!(V::AbstractApproximator, π, model::AbstractEnvironmentModel; γ::Float64=0.9, θ::Float64=1e-4, max_iter=typemax(Int))

See more details at Section (4.3) on Page 80 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
function policy_iteration!(
    ;
    V::AbstractApproximator,
    π::AbstractPolicy,
    model::AbstractEnvironmentModel,
    γ::Float64 = 0.9,
    θ::Float64 = 1e-4,
    max_iter = typemax(Int),
)
    for i = 1:max_iter
        @debug "iteration: $i"
        policy_evaluation!(; V = V, π = π, model = model, γ = γ, θ = θ)
        policy_improvement!(; V = V, π = π, model = model, γ = γ) && break
    end
end

"""
    value_iteration!(V::AbstractApproximator, model::AbstractEnvironmentModel; γ::Float64=0.9, θ::Float64=1e-4, max_iter=typemax(Int))

See more details at Section (4.4) on Page 83 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
function value_iteration!(
    ;
    V::AbstractApproximator,
    model::AbstractEnvironmentModel,
    γ::Float64 = 0.9,
    θ::Float64 = 1e-4,
    max_iter = typemax(Int),
)
    states, actions = 1:length(get_observation_space(model)), 1:length(get_action_space(model))
    for i = 1:max_iter
        Δ = 0.0
        for s in states
            v = maximum(sum(p * (r + γ * V(s′)) for (s′, r, p) in model(s, a)) for a in actions)
            error = v - V(s)
            update!(V, s => error)
            Δ = max(Δ, abs(error))
        end
        Δ < θ && break
    end
end


text:
using Oceananigans.Solvers: solve_poisson_equation!

function ∇²!(grid, f, ∇²f)
    @loop_xyz i j k grid begin
        @inbounds ∇²f[i, j, k] = ∇²(i, j, k, grid, f)
    end
end

function pressure_solver_instantiates(FT, Nx, Ny, Nz, planner_flag)
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), length=(100, 100, 100))
    solver = PressureSolver(CPU(), grid, PressureBoundaryConditions(grid), planner_flag)
    return true  # Just making sure the PressureSolver does not error/crash.
end

function poisson_ppn_planned_div_free_cpu(FT, Nx, Ny, Nz, planner_flag)
    arch = CPU()
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6))
    fbcs = TracerBoundaryConditions(grid)
    pbcs = PressureBoundaryConditions(grid)
    solver = PressureSolver(arch, grid, fbcs)

    RHS = CellField(FT, arch, grid, fbcs)
    interior(RHS) .= rand(Nx, Ny, Nz)
    interior(RHS) .= interior(RHS) .- mean(interior(RHS))

    RHS_orig = deepcopy(RHS)
    solver.storage .= interior(RHS)
    solve_poisson_equation!(solver, grid)

    ϕ   = CellField(FT, arch, grid, pbcs)
    ∇²ϕ = CellField(FT, arch, grid, pbcs)

    interior(ϕ) .= real.(solver.storage)

    fill_halo_regions!(ϕ, arch)
    ∇²!(grid, ϕ, ∇²ϕ)

    fill_halo_regions!(∇²ϕ, arch)

    interior(∇²ϕ) ≈ interior(RHS_orig)
end

function poisson_pnn_planned_div_free_cpu(FT, Nx, Ny, Nz, planner_flag)
    arch = CPU()
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6), topology=(Periodic, Bounded, Bounded))
    fbcs = TracerBoundaryConditions(grid)
    pbcs = PressureBoundaryConditions(grid)
    solver = PressureSolver(arch, grid, fbcs)

    RHS = CellField(FT, arch, grid, fbcs)
    interior(RHS) .= rand(Nx, Ny, Nz)
    interior(RHS) .= interior(RHS) .- mean(interior(RHS))

    RHS_orig = deepcopy(RHS)

    solver.storage .= interior(RHS)

    solve_poisson_equation!(solver, grid)

    ϕ   = CellField(FT, arch, grid, pbcs)
    ∇²ϕ = CellField(FT, arch, grid, pbcs)

    interior(ϕ) .= real.(solver.storage)

    fill_halo_regions!(ϕ, arch)
    ∇²!(grid, ϕ, ∇²ϕ)

    fill_halo_regions!(∇²ϕ, arch)

    interior(∇²ϕ) ≈ interior(RHS_orig)
end

function poisson_ppn_planned_div_free_gpu(FT, Nx, Ny, Nz)
    arch = GPU()
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6))
    pbcs = PressureBoundaryConditions(grid)
    solver = PressureSolver(arch, grid, pbcs)

    RHS = rand(Nx, Ny, Nz)
    RHS .= RHS .- mean(RHS)
    RHS = CuArray(RHS)

    RHS_orig = copy(RHS)

    solver.storage .= RHS

    # Performing the permutation [a, b, c, d, e, f] -> [a, c, e, f, d, b]
    # in the z-direction in preparation to calculate the DCT in the Poisson
    # solver.
    solver.storage .= cat(solver.storage[:, :, 1:2:Nz], solver.storage[:, :, Nz:-2:2]; dims=3)

    solve_poisson_equation!(solver, grid)

    # Undoing the permutation made above to complete the IDCT.
    solver.storage .= CuArray(reshape(permutedims(cat(solver.storage[:, :, 1:Int(Nz/2)],
                                                      solver.storage[:, :, end:-1:Int(Nz/2)+1]; dims=4), (1, 2, 4, 3)), Nx, Ny, Nz))

    ϕ   = CellField(FT, arch, grid, pbcs)
    ∇²ϕ = CellField(FT, arch, grid, pbcs)

    interior(ϕ) .= real.(solver.storage)

    fill_halo_regions!(ϕ, arch)
    ∇²!(grid, ϕ.data, ∇²ϕ.data)

    fill_halo_regions!(∇²ϕ, arch)
    interior(∇²ϕ) ≈ RHS_orig
end

function poisson_pnn_planned_div_free_gpu(FT, Nx, Ny, Nz)
    arch = GPU()
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), length=(1.0, 2.5, 3.6), topology=(Periodic, Bounded, Bounded))
    pbcs = PressureBoundaryConditions(grid)
    solver = PressureSolver(arch, grid, pbcs)

    RHS = rand(Nx, Ny, Nz)
    RHS .= RHS .- mean(RHS)
    RHS = CuArray(RHS)

    RHS_orig = copy(RHS)

    storage = solver.storage.storage1
    storage .= RHS

    storage .= cat(storage[:, :, 1:2:Nz], storage[:, :, Nz:-2:2]; dims=3)
    storage .= cat(storage[:, 1:2:Ny, :], storage[:, Ny:-2:2, :]; dims=2)

    solve_poisson_equation!(solver, grid)

    ϕ   = CellField(FT, arch, grid)
    ∇²ϕ = CellField(FT, arch, grid)

    # Indices used when we need views to permuted arrays where the odd indices
    # are iterated over first followed by the even indices.
    p_y_inds = [1:2:Ny..., Ny:-2:2...] |> CuArray
    p_z_inds = [1:2:Nz..., Nz:-2:2...] |> CuArray

    ϕ_p = view(interior(ϕ), 1:Nx, p_y_inds, p_z_inds)

    @. ϕ_p = real(storage)

    fill_halo_regions!(ϕ, arch)
    ∇²!(grid, ϕ.data, ∇²ϕ.data)

    fill_halo_regions!(∇²ϕ, arch)
    interior(∇²ϕ) ≈ RHS_orig
end

"""
    poisson_ppn_recover_sine_cosine_solution(FT, Nx, Ny, Nz, Lx, Ly, Lz, mx, my, mz)

Test that the Poisson solver can recover an analytic solution. In this test, we
are trying to see if the solver can recover the solution

    ``\\Psi(x, y, z) = cos(\\pi m_z z / L_z) sin(2\\pi m_y y / L_y) sin(2\\pi m_x x / L_x)``

by giving it the source term or right hand side (RHS), which is

    ``f(x, y, z) = \\nabla^2 \\Psi(x, y, z) =
    -((\\pi m_z / L_z)^2 + (2\\pi m_y / L_y)^2 + (2\\pi m_x/L_x)^2) \\Psi(x, y, z)``.
"""
function poisson_ppn_recover_sine_cosine_solution(FT, Nx, Ny, Nz, Lx, Ly, Lz, mx, my, mz)
    grid = RegularCartesianGrid(FT, size=(Nx, Ny, Nz), length=(Lx, Ly, Lz))
    solver = PressureSolver(CPU(), grid, TracerBoundaryConditions(grid))

    xC, yC, zC = grid.xC, grid.yC, grid.zC
    xC = reshape(xC, (Nx, 1, 1))
    yC = reshape(yC, (1, Ny, 1))
    zC = reshape(zC, (1, 1, Nz))

    Ψ(x, y, z) = cos(π*mz*z/Lz) * sin(2π*my*y/Ly) * sin(2π*mx*x/Lx)
    f(x, y, z) = -((mz*π/Lz)^2 + (2π*my/Ly)^2 + (2π*mx/Lx)^2) * Ψ(x, y, z)

    @. solver.storage = f(xC, yC, zC)
    solve_poisson_equation!(solver, grid)
    ϕ = real.(solver.storage)

    error = norm(ϕ - Ψ.(xC, yC, zC)) / √(Nx*Ny*Nz)

    @info "Error (ℓ²-norm), $FT, N=($Nx, $Ny, $Nz), m=($mx, $my, $mz): $error"

    isapprox(ϕ, Ψ.(xC, yC, zC), rtol=5e-2)
end

@testset "Pressure solvers" begin
    @info "Testing pressure solvers..."

    @testset "Pressure solver instantiation" begin
        @info "  Testing pressure solver instantiation..."

        for FT in float_types
            @test pressure_solver_instantiates(FT, 32, 32, 32, FFTW.ESTIMATE)
            @test pressure_solver_instantiates(FT, 1,  32, 32, FFTW.ESTIMATE)
            @test pressure_solver_instantiates(FT, 32,  1, 32, FFTW.ESTIMATE)
            @test pressure_solver_instantiates(FT,  1,  1, 32, FFTW.ESTIMATE)
        end
    end

    @testset "Divergence-free solution [CPU]" begin
        @info "  Testing divergence-free solution [CPU]..."

        for N in [7, 10, 16, 20]
            for FT in float_types
                for planner_flag in (FFTW.ESTIMATE, FFTW.MEASURE)
                    @test poisson_ppn_planned_div_free_cpu(FT, N, N, N, planner_flag)
                    @test poisson_ppn_planned_div_free_cpu(FT, 1, N, N, planner_flag)
                    @test poisson_ppn_planned_div_free_cpu(FT, N, 1, N, planner_flag)
                    @test poisson_ppn_planned_div_free_cpu(FT, 1, 1, N, planner_flag)
                end
            end
        end

        Ns = [5, 11, 20, 32]
        for Nx in Ns, Ny in Ns, Nz in Ns, FT in float_types
            @test poisson_ppn_planned_div_free_cpu(FT, Nx, Ny, Nz, FFTW.ESTIMATE)
            @test poisson_pnn_planned_div_free_cpu(FT, Nx, Ny, Nz, FFTW.ESTIMATE)
        end
    end

    @testset "Divergence-free solution [GPU]" begin
        @info "  Testing divergence-free solution [GPU]..."
        @hascuda begin
            for FT in [Float64]
                @test poisson_ppn_planned_div_free_gpu(FT, 16, 16, 16)
                @test poisson_ppn_planned_div_free_gpu(FT, 32, 32, 32)
                @test poisson_ppn_planned_div_free_gpu(FT, 32, 32, 16)
                @test poisson_ppn_planned_div_free_gpu(FT, 16, 32, 24)

                @test poisson_pnn_planned_div_free_gpu(FT, 16, 16, 16)
                @test poisson_pnn_planned_div_free_gpu(FT, 32, 32, 32)
                @test poisson_pnn_planned_div_free_gpu(FT, 32, 32, 16)
                @test poisson_pnn_planned_div_free_gpu(FT, 16, 32, 24)
            end
        end
    end

    @testset "Analytic solution reconstruction" begin
        @info "  Testing analytic solution reconstruction..."
        for N in [32, 48, 64], m in [1, 2, 3]
            @test poisson_ppn_recover_sine_cosine_solution(Float64, N, N, N, 100, 100, 100, m, m, m)
        end
    end
end


text:
type CtrDrbg  <: AbstractRNG
    data::Ptr{Void}
    entropy::Nullable{Entropy}

    function CtrDrbg()
        ctx = new()
        ctx.data = Libc.malloc(1000)  # 344
        ccall((:mbedtls_ctr_drbg_init, MBED_CRYPTO), Void, (Ptr{Void},), ctx.data)
        finalizer(ctx, ctx->begin
            ccall((:mbedtls_ctr_drbg_free, MBED_CRYPTO), Void, (Ptr{Void},), ctx.data)
            Libc.free(ctx.data)
        end
        )
        ctx
    end

    CrtDrbg(data) = new(data)
end

function f_rng(c_ctx, c_buf, sz)
    jl_ctx = unsafe_pointer_to_objref(c_ctx)
    jl_buf = unsafe_wrap(Array, c_buf, sz, false)
    rand!(jl_ctx, jl_buf)
    MBED_SUCCESS
end

function seed!(rng::CtrDrbg, entropy, pdata)
    rng.entropy = Nullable(entropy)
    entropy_func = cglobal((:mbedtls_entropy_func, MBED_CRYPTO))
    @err_check ccall((:mbedtls_ctr_drbg_seed, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Csize_t),
        rng.data, entropy_func, entropy.data, pdata, sizeof(pdata))
    rng
end

seed!(rng::CtrDrbg, entropy) = seed!(rng, entropy, UInt8[])

function Base.rand!(rng::CtrDrbg, buf::Array)
    @err_check ccall((:mbedtls_ctr_drbg_random, MBED_CRYPTO), Cint,
            (Ptr{Void}, Ptr{Void}, Csize_t),
        rng.data, buf, sizeof(buf))
    buf
end

function Base.rand(rng::CtrDrbg, size::Integer)
    buf = Vector{UInt8}(size)
    rand!(rng, buf)
end

function __ctr_drbg__init__()
    const global c_rng=cfunction(f_rng, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t))
end


text:
module RandomTestMatrices

export randsym, randantisym, randposdef, randstable, randnormal, randorth, randunitary

#=

Random test matrices. These are for testing and debugging purposes, for random matrix ensembles with known limiting distributions, see the package RandomMatrices.jl .

=#

"""
randsym(d)

	Gaussian random symmetric matrix of dimension ``d``.
"""
function randsym(d)
    y = randn(d, d)
	triu(y) + triu(y,1)'
end

"""
randantisym(d)

    Anti symmetric matrix
"""
function randantisym(d)
	y = randn(d, d)
	triu(y,1) - triu(y,1)'
end

"""
randposdef(d, n = d)

	Random positive definite matrix of dimension ``d`` from the Laguerre ensemble ``n``.
"""
function randposdef(d,n = d)
	x = randn(d, n)
	x'* x/n
end

"""
randstable(d)

	Random stable matrix (matrix with eigenvalues with negative real part) with dimension ``d``.
"""
function randstable(d)
	# positive definite matrix
	x = randn(d, d)
	a = x'*x/d

	# antisymmetric matrix
	b = randantisym(d)

	# return stable matrix
	b - a
end

"""
randunitary(d)

	Random unitary matrix of dimension ``d``.
"""
randunitary(d) = expm(im* randposdef(d))

"""
randorth(d)

	Orthogonal matrix drawn according to the Haar measure on the group of orthogonal matrices.
"""
randorth(d) = qr(randn(d,d))[1]

function tridiag(di, u, l)
	assert(length(di) -1 == length(u) == length(l))
	M = diagm(di)
	for i in 1:length(u)
	  M[i,i+1] = l[i]
  	  M[i+1,i] = u[i]
	end
	M

end

"""
randnormal(d)

	Random normal matrix of dimension ``d``.
"""
function randnormal(d)
	Q = randorth(d)
	B = schur(randn(d,d))[1] #same eigenvalues as a randn(d,d)
	m = d
	alpha = zeros(d)
	beta = zeros(d-1)
	while m > 1
	   s = abs(B[m-1,m-1]) + abs(B[m,m])
	   if s + abs(B[m,m-1]) > s # if significant offdiagonal value: evaluate submatrix m-1:m;
		spur = B[m,m] + B[m-1,m-1]
		dis2 = abs2(B[m,m]- B[m-1,m-1]) + 4.0*B[m,m-1]*B[m-1,m] # = spur^2 - 4det
		alpha[m] = alpha[m-1] = 0.5*spur

		beta[m-1] = 0.5*sqrt(-dis2)
		m -= 2
	   else
	        alpha[m] = B[m,m]
	        beta[m-1] = 0
	        m -= 1
	   end

	end
	if (m == 1)
		alpha[1] = B[1,1]
	end

	Q' * tridiag(alpha, beta, -beta) * Q

end

end # module


text:
using ProjAssistant
using Test
using Plots
ENV["GKSwstype"]="nul" # avoid external display

@testset "ProjAssistant.jl" begin
    include("cache_tests.jl")
    include("proj_gen_tests.jl")
end


text:
function do_nothing()
end

function show_node(node, search_helper)
    # search_helper.depth
    # show(stdio, node)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using PowerSystems
using TimeSeries

base_dir = PowerSystems.download(PowerSystems.TestData; branch = "master");

sys = System(joinpath(base_dir, "matpower", "case5_re.m"))
sys

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
@inline Base.@propagate_inbounds Base.getindex(a::AbstractVec,I::Integer) = getfield(a,I)
Base.IndexStyle(a::Type{<:AbstractVec}) = Base.IndexLinear()

@inline xpos(a::AbstractVec) = a.x
@inline ypos(a::AbstractVec) = a.y
@inline zpos(a::AbstractVec) = a.z

#pos of a node returns a tuple with the positions
pos(a::AbstractVec) = (xpos(a),ypos(a),zpos(a))

Base.norm(a::AbstractVec) = sqrt(muladd(xpos(a), xpos(a), muladd(ypos(a), ypos(a), zpos(a)^2)))
distance(a::AbstractVec,b::AbstractVec) = sqrt((xpos(b)-xpos(a))^2 + (ypos(b)-ypos(a))^2 + (zpos(b)-zpos(a))^2)

Base.:+(a::AbstractVec,b::AbstractVec) = Vec(xpos(a)+xpos(b), ypos(a)+ypos(b), zpos(a)+zpos(b))
Base.:-(a::AbstractVec,b::AbstractVec) = Vec(xpos(a)-xpos(b), ypos(a)-ypos(b), zpos(a)-zpos(b))
Base.:-(a::AbstractVec) = Vec(-xpos(a), -ypos(a), -zpos(a))
Base.:*(a::Number,b::AbstractVec) = Vec(a*xpos(b), a*ypos(b), a*zpos(b))
Base.:*(b::AbstractVec,a::Number) = Vec(a*xpos(b), a*ypos(b), a*zpos(b))
Base.:/(b::AbstractVec,a::Number) = Vec(xpos(b)/a, ypos(b)/a, zpos(b)/a)

Base.dot(a::AbstractVec,b::AbstractVec) = muladd(xpos(a), xpos(b), muladd(ypos(a), ypos(b), zpos(a)*zpos(b)))

Base.cross(a::AbstractVec,b::AbstractVec) = Vec(ypos(a)*zpos(b) - zpos(a)*ypos(b), zpos(a)*xpos(b) - xpos(a)*zpos(b), xpos(a)*ypos(b) - ypos(a)*xpos(b))

struct Vec{T<:Number} <: AbstractVec{T}
    x::T
    y::T
    z::T
end

Base.size(a::Vec) = (3,)
Base.length(a::Vec) = 3
Base.linearindices(a::Vec) = Base.OneTo(3)
Base.indices(a::Vec) = (Base.OneTo(3),)
Base.zero(a::Type{Vec{T}}) where {T} = Vec{T}(zero(T),zero(T),zero(T))
Vec(x,y,z) = Vec(promote(x,y,z)...)

struct Vec2D{T<:Number} <: AbstractVec{T}
    x::T
    y::T
end

Base.size(a::Vec2D) = (2,)
Base.length(a::Vec2D) = 2
Base.linearindices(a::Vec2D) = Base.OneTo(2)
Base.indices(a::Vec2D) = (Base.OneTo(2),)
Base.zero(a::Type{Vec2D{T}}) where {T} = Vec2D{T}(zero(T),zero(T))

Vec2D(x,y) = Vec2D(promote(x,y)...)

@inline zpos(a::Vec2D{T}) where {T} = zero(T)

Base.:+(a::Vec2D{T},b::Vec2D{T2}) where {T,T2} = Vec2D{promote_type(T,T2)}(xpos(a)+xpos(b), ypos(a)+ypos(b))
Base.:-(a::Vec2D{T},b::Vec2D{T2}) where {T,T2} = Vec2D{promote_type(T,T2)}(xpos(a)-xpos(b), ypos(a)-ypos(b))
Base.:-(a::Vec2D{T}) where {T} = Vec2D{T}(-xpos(a), -ypos(a))
Base.:*(a::T,b::Vec2D{T2}) where {T<:Number,T2} = Vec2D{promote_type(T,T2)}(a*xpos(b), a*ypos(b))
Base.:*(b::Vec2D{T},a::T2) where {T,T2<:Number} = Vec2D{promote_type(T,T2)}(a*xpos(b), a*ypos(b))
Base.:/(b::Vec2D{T},a::T2) where {T,T2<:Number} = Vec2D{promote_type(T,T2)}(xpos(b)/a, ypos(b)/a)

Base.dot(a::Vec2D,b::Vec2D) = muladd(xpos(a), xpos(b), ypos(a)*ypos(b))
Base.norm(a::Vec2D) = sqrt(muladd(xpos(a), xpos(a), ypos(a)^2))
Base.cross(a::Vec2D,b::Vec2D) = xpos(a)*ypos(b) - ypos(a)*xpos(b)

struct Vec1D{T<:Number} <: AbstractVec{T}
    x::T
end

Base.size(a::Vec1D) = (1,)
Base.length(a::Vec1D) = 1
Base.linearindices(a::Vec1D) = Base.OneTo(1)
Base.indices(a::Vec1D) = (Base.OneTo(1),)
Base.zero(a::Type{Vec1D{T}}) where {T} = Vec1D{T}(zero(T))

@inline ypos(a::Vec1D{T}) where {T} = zero(T)
@inline zpos(a::Vec1D{T}) where {T} = zero(T)

Base.:+(a::Vec1D{T},b::Vec1D{T2}) where {T,T2} = Vec1D{promote_type(T,T2)}(xpos(a)+xpos(b))
Base.:-(a::Vec1D{T},b::Vec1D{T2}) where {T,T2} = Vec1D{promote_type(T,T2)}(xpos(a)-xpos(b))
Base.:-(a::Vec1D{T}) where {T} = Vec1D{T}(-xpos(a))
Base.:*(a::T,b::Vec1D{T2}) where {T<:Number,T2} = Vec1D{promote_type(T,T2)}(a*xpos(b))
Base.:*(b::Vec1D{T},a::T2) where {T,T2<:Number} = Vec1D{promote_type(T,T2)}(a*xpos(b))
Base.:/(b::Vec1D{T},a::T2) where {T,T2<:Number} = Vec1D{promote_type(T,T2)}(xpos(b)/a)

Base.dot(a::Vec1D,b::Vec1D) = xpos(a)*xpos(b)
Base.norm(a::Vec1D) = abs(xpos(a))
Base.cross(a::Vec1D,b::Vec1D) = 0


text:
using Documenter, ACEtb

makedocs(sitename="ACEtb.jl Documentation",
         pages = [
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "Getting Started" => "gettingstarted.md",
        "Installation" => "installation.md",
        "ACEtb Model" => [
            "What is ACEtb Model" => "acetb_intro.md",
            "How to fit ACEtb Model" => "fitting.md",
            "Basic Usage" => "acetb_usage.md",
        ],
        "ACEtb API: `lib_acetb`" => [
            "DFTB+ Integration" => "dftbplus_acetb.md",
        ],
        "Developer Docs" => "devel.md"
        ])

# deploydocs(
#     repo = "github.com/ACEsuit/ACEtb.jl.git",
# )


text:
include("test_images/synth_images.jl")

@testset "synthetic images" begin
    for (case, image) in synth_imgs
        @testset "$(case)" begin
            expected = collect(PNGFiles._prepare_buffer(image))
            fpath = joinpath(PNG_TEST_PATH, "test_img_$(case).png")
            @testset "write" begin
                open(io->PNGFiles.save(io, image), fpath, "w") #test IO method
                @test PNGFiles.save(fpath, image) === nothing
            end
            @testset "read" begin
                global read_in_pngf = PNGFiles.load(fpath)
                @test typeof(read_in_pngf) <: AbstractMatrix
            end
            @testset "compare" begin
                @test eltype(eltype(expected)) <: N0f8 ? expected == read_in_pngf : all(expected .≈ read_in_pngf)
            end
            global read_in_immag = _standardize_grayness(ImageMagick.load(fpath))
            @testset "$(case): ImageMagick read type equality" begin
                # The lena image is Grayscale saved as RGB...
                @test eltype(_standardize_grayness(read_in_pngf)) == eltype(read_in_immag)
            end
            @testset "$(case): ImageMagick read values equality" begin
                imdiff_val = imdiff(read_in_pngf, read_in_immag)
                onfail(@test imdiff_val < 0.01) do
                    PNGFiles._inspect_png_read(fpath)
                    _add_debugging_entry(fpath, case, imdiff_val)
                end
            end
            path, ext = splitext(fpath)
            newpath = path * "_new" * ext
            PNGFiles.save(newpath, read_in_pngf)
            @testset "$(case): Fpath save IO is idempotent" begin
                @test all(read_in_pngf .≈ PNGFiles.load(newpath))
                @test all(read_in_pngf .≈ open(io->PNGFiles.load(io), newpath))
            end
            open(io->PNGFiles.save(io, read_in_pngf), newpath, "w") #test IO method
            @testset "$(case): Stream save IO is idempotent" begin
                @test all(read_in_pngf .≈ PNGFiles.load(newpath))
                @test all(read_in_pngf .≈ open(io->PNGFiles.load(io), newpath))
            end
        end
    end

    for (case, func_in, image) in edge_case_imgs
        @testset "$(case)" begin
            fpath = joinpath(PNG_TEST_PATH, "test_img_$(case).png")
            @testset "write" begin
                @test PNGFiles.save(fpath, image) === nothing
            end
            @testset "read" begin
                global read_in_pngf = PNGFiles.load(fpath)
                @test read_in_pngf isa Matrix
            end
            @testset "compare" begin
                @test all(read_in_pngf .== func_in(image))
            end
            global read_in_immag = _standardize_grayness(ImageMagick.load(fpath))
            @testset "$(case): ImageMagick read type equality" begin
                @test eltype(read_in_pngf) == eltype(read_in_immag)
            end
            @testset "$(case): ImageMagick read values equality" begin
                imdiff_val = imdiff(read_in_pngf, read_in_immag)
                onfail(@test imdiff_val < 0.01) do
                    PNGFiles._inspect_png_read(fpath)
                    _add_debugging_entry(fpath, case, imdiff_val)
                end
            end
            path, ext = splitext(fpath)
            newpath = path * "_new" * ext
            PNGFiles.save(newpath, read_in_pngf)
            @testset "$(case): IO is idempotent" begin
                @test imdiff(read_in_pngf, PNGFiles.load(newpath)) < 0.01
                @test imdiff(read_in_pngf, open(io->PNGFiles.load(io), newpath)) < 0.01
            end
        end
    end
end


text:
export ARCSpectral

function ARCSpectral(nlp 		:: AbstractNLPModel,
               		 nlpstop 	:: NLPStopping;
					 kwargs...
               		 )

	return TRARC2(nlp,
				  nlpstop;
				  TR = TrustRegion(10.0),
				  c = Combi(ARCTR.hessian_dense,ARCTR.PDataSpectral{eltype(nlp.meta.x0)},ARCTR.solve_modelARCDiag,ARCTR.preprocessSpectral,ARCTR.decreaseFact,ARCTR.Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
using LinearAlgebra

function loadTitlePhrases(fileNameIn::String)::Tuple{Vector{Vector{String}}, Dict{String,Int}}
    fileIn = open(fileNameIn)
    phraseGroups = Vector{Vector{String}}()
    phraseFreqs = Dict{String,Int}()
    for line in readlines(fileIn)
        phrases = split(line)
        if length(phrases) > 0
            push!(phraseGroups, phrases)
            for phrase in phrases
                phraseFreqs[phrase] = get(phraseFreqs, phrase, 0) + 1
            end
        end
    end
    close(fileIn)
    (phraseGroups, phraseFreqs)
end

function getSubPhrasePosition(words1::Vector{SubString{String}},words2::Vector{SubString{String}})::Tuple{Float64, Float64}
    n1 = length(words1)
    n2 = length(words2)
    if n1 <= n2
        return (-1.0,-1.0)
    end
    for i in 1:n1-n2+1
        matched = true
        for j in 1:n2
            if words1[i+j-1] != words2[j]
                matched = false
                break
            end
        end
        if matched
            return (i/n1,(n1-i-n2+1)/n1)
        end
    end
    return (-1.0,-1.0)
end

SubPhrasePositions = Dict{String,Tuple{Float64,Float64}}
function findSubPhrases(phraseFreqs::Dict{String,Int})::Tuple{Dict{String,SubPhrasePositions}, Dict{String, Int}}
    phrases = [phrase for (phrase, freq) in phraseFreqs]
    numPhrases = length(phrases)
    phraseWords = [split(phrases[i], "-") for i in 1:numPhrases]

    numCPUs = length(Sys.cpu_info())
    jobs = Channel{Tuple{Int,Int}}(numCPUs)
    jobOutputs = Channel{Dict{String, SubPhrasePositions}}(numCPUs)
    function runJob(iCPU::Int)
        myOutputs = Dict{String, SubPhrasePositions}()
        for job in jobs
            for phraseID1 in job[1]:job[2]
                phrase1 = phrases[phraseID1]
                words1 = phraseWords[phraseID1]
                subPhrases1 = SubPhrasePositions()
                for phraseID2 in 1:numPhrases
                    words2 = phraseWords[phraseID2]
                    subPhrasePosition = getSubPhrasePosition(words1, words2)
                    if subPhrasePosition[1] >= 0.0
                        phrase2 = phrases[phraseID2]
                        subPhrases1[phrase2] = subPhrasePosition
                    end
                end
                myOutputs[phrase1] = subPhrases1
            end
        end
        put!(jobOutputs, myOutputs)
    end
    function makeJobs(batchSize::Int)
        for i in 1:batchSize:numPhrases
            put!(jobs, (i, min(i + batchSize - 1, numPhrases)))
        end
    end
    bind(jobs, @async makeJobs(100))
    for iCPU in 1:numCPUs
        Threads.@spawn runJob(iCPU)
        # runJob(iCPU)
    end

    subPhrases = Dict{String, SubPhrasePositions}()
    for iCPU in 1:numCPUs
        localSubPhrases = take!(jobOutputs)
        for (phrase,subPhrasePositions) in localSubPhrases
            subPhrases[phrase] = subPhrasePositions
        end
    end

    subPhraseFreqs = copy(phraseFreqs)
    for (phrase,subPhrasePositions) in subPhrases
        freq = phraseFreqs[phrase]
        for (subPhrase, position) in subPhrasePositions
            subPhraseFreqs[subPhrase] += freq
        end
    end

    (subPhrases, subPhraseFreqs)
end

function getScores(fileNameIn::String, minCount::Int)::Tuple{Array{Float64,2},Array{String,1},Dict{String,Int64}}
    phraseGroups, phraseFreqs = loadTitlePhrases(fileNameIn)
    subPhrases, subPhraseFreqs = findSubPhrases(phraseFreqs)

    phrases = Vector{String}()
    phraseIDs = Dict{String,Int}()
    for (phrase, freq) in subPhraseFreqs
        if freq >= minCount
            push!(phrases, phrase)
            phraseIDs[phrase] = length(phrases)
        end
    end

    numPhrases = length(phrases)
    scores = fill(0.0, (numPhrases, numPhrases))
    for titlePhrases in phraseGroups
        highFreqs = Vector{Bool}(undef, length(titlePhrases))
        subHighFreqs = Vector{Bool}(undef, length(titlePhrases))
        for (j, phraseJ) in enumerate(titlePhrases)
            highFreqs[j] = haskey(phraseIDs, phraseJ)
            subHighFreqs[j] = false
            for (subPhrase, position) in subPhrases[phraseJ]
                if haskey(phraseIDs, subPhrase)
                    subHighFreqs[j] = true
                end
            end
        end
        for (j, phraseJ) in enumerate(titlePhrases)
            if subHighFreqs[j]
                for (subPhrase1, position1) in subPhrases[phraseJ]
                    if haskey(phraseIDs, subPhrase1)
                        id1 = phraseIDs[subPhrase1]
                        for (subPhrase2, position2) in subPhrases[phraseJ]
                            if subPhrase1 == subPhrase2
                                continue
                            end
                            if haskey(phraseIDs, subPhrase2)
                                id2 = phraseIDs[subPhrase2]
                                ldist = abs(position1[1] - position2[1])
                                rdist = abs(position1[2] - position2[2])
                                dist = max(ldist, rdist)
                                @assert dist > 0.0
                                score = 1.0 / dist
                                # make sure the following update is one sided
                                scores[id1, id2] += score
                            end
                        end
                        if highFreqs[j]
                            id2 = phraseIDs[phraseJ]
                            dist = max(position1[1], position1[2])
                            @assert dist > 0.0
                            score = 1.0 / dist
                            scores[id1, id2] += score
                            scores[id2, id1] += score
                        end
                    end
                end
            end
            for k in j+1:length(titlePhrases)
                phraseK = titlePhrases[k]
                if highFreqs[j]
                    idJ = phraseIDs[phraseJ]
                    if highFreqs[k]
                        idK = phraseIDs[phraseK]
                        score = 1.0 / (k - j)
                        scores[idJ, idK] += score
                        scores[idK, idJ] += score
                    end
                    if subHighFreqs[k]
                        for (subPhraseK, positionK) in subPhrases[phraseK]
                            if haskey(phraseIDs, subPhraseK)
                                idK = phraseIDs[subPhraseK]
                                score = 1.0 / (k - j + positionK[1])
                                scores[idJ, idK] += score
                                scores[idK, idJ] += score
                            end
                        end
                    end
                end
                if subHighFreqs[j]
                    for (subPhraseJ, positionJ) in subPhrases[phraseJ]
                        if haskey(phraseIDs, subPhraseJ)
                            idJ = phraseIDs[subPhraseJ]
                            if highFreqs[k]
                                idK = phraseIDs[phraseK]
                                score = 1.0 / (k - j + positionJ[2])
                                scores[idJ, idK] += score
                                scores[idK, idJ] += score
                            end
                            if subHighFreqs[k]
                                for (subPhraseK, positionK) in subPhrases[phraseK]
                                    if haskey(phraseIDs, subPhraseK)
                                        idK = phraseIDs[subPhraseK]
                                        score = 1.0 / (k - j + positionJ[2] + positionK[1])
                                        scores[idJ, idK] += score
                                        scores[idK, idJ] += score
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    A = fill(0.0, (numPhrases, numPhrases))
    y = fill(0.0, numPhrases)
    for i in 1:numPhrases
        for j in 1:numPhrases
            if j == i
                continue
            end
            if scores[i, j] > 0.0
                A[i,i] += 1.0
                A[i,j] += 1.0
                y[i] += log(scores[i,j])
            end
        end
    end
    bias = A \ y
    for i in 1:numPhrases
        scores[i,i] = 1.0
        for j in 1:numPhrases
            if j == i
                continue
            end
            if scores[i, j] > 0.0
                unbiasedScore = log(scores[i,j]) - bias[i] - bias[j]
                scores[i,i] += abs(unbiasedScore)
                scores[i,j] = unbiasedScore
            end
        end
    end

    # this enforce symmetry to eliminate numerical corruption of symmetry
    scores = 0.5(scores + scores')

    (scores, phrases, phraseIDs)
end

function trainWordEmbedding(fileNameIn::String, fileNameOut::String, numDims::Int, minCount::Int)
    scores, phrases, phraseIDs = getScores(fileNameIn, minCount)

    eigScores = eigen(Symmetric(scores))
    eigOrder = sortperm(eigScores.values, rev=true)

    # note: eigScores.vectors * diagm(eigScores.values) * eigScores.vectors' == scores
    numPhrases = length(phrases)
    w = fill(0.0, min(numDims, numPhrases))
    for k in 1:min(numDims, numPhrases)
        w[k] = sqrt(eigScores.values[eigOrder[k]])
    end
    fileOut = open(fileNameOut, "w")
    println(fileOut, "$numPhrases $numDims")
    for i in 1:numPhrases
        phrase = phrases[i]
        print(fileOut, phrase)
        for k in 1:min(numDims, numPhrases)
            value = eigScores.vectors[i,eigOrder[k]] * w[k]
            print(fileOut, " $value")
        end
        println(fileOut)
    end
    close(fileOut)
end

max_stars_repo_path: src/x3 aspect_aware_word_embedding.jl

max_issues_repo_path: src/x3 aspect_aware_word_embedding.jl

max_forks_repo_path: src/x3 aspect_aware_word_embedding.jl


text:
module AvalonPackets

#  This module has functions to create and manipulate Avalon Packets.

export do_transaction, tx_packet, get_readback_data
export SEQUENTIAL_WRITE, analyze_response

using Printf

#  Special Packet characters
SOP = 0x7a
EOP = 0x7b
CHANNEL = 0x7c
ESC = 0x7d

#  There is a 3-step process to creating a packet.
#  do_transaction:  Simply glues an Avalon header to the data.
#  tx_packet:  Inserts special characters where required.
#

#  Transaction op codes (for reference):
SEQUENTIAL_WRITE = 0x04
SEQUENTIAL_READ  = 0x14
NON_SEQUENTIAL_WRITE = 0x00
NON_SEQUENTIAL_READ = 0x10

HEADER_LENGTH = 8
RESPONSE_LENGTH = 4

#  This function simply attaches an Avalon packet header to the data.
function do_transaction(type, size, address, data)

    if ((type == SEQUENTIAL_WRITE) | (type == NON_SEQUENTIAL_WRITE) | (type == SEQUENTIAL_READ) | (type == NON_SEQUENTIAL_READ))
      transaction = zeros(UInt8, (size + HEADER_LENGTH))
#transaction = Array{UInt8,1}(UndefInitializer(), size + HEADER_LENGTH)
    end

    #  Make header:
    transaction[1] = type
    transaction[2] = 0x00
    transaction[3] = (size >> 8) & 0xff  #  Works exactly like C?
    transaction[4] = (size & 0xff)
    transaction[5] = (address >> 24) & 0xff
    transaction[6] = (address >> 16) & 0xff
    transaction[7] = (address >>  8) & 0xff
    transaction[8] = (address & 0xff)

    #  Load the data into the transaction:
    if((type == SEQUENTIAL_WRITE) | (type == NON_SEQUENTIAL_WRITE))
    for i in 9:(size + HEADER_LENGTH)
        transaction[i] = data[i - HEADER_LENGTH]
    end
end

    length = size + 8
    return transaction, length

end # do_transaction

#  Complete the Transmit packet.  Pass in the result from do_transaction.
#  The address is already included in the data's header.
#  This function will insert the special characters.
function tx_packet(data, data_length)
#  Response data is a fixed value of 2 bytes to form a 16 bit unsigned integer,
#  which is the number of bytes transmitted or received.  This is the data
#  returned by the Avalon response packet.
#  There are 4 bytes total in the response.  The first byte is the transaction
#  op code with the MSB bit inverted.  The 2nd byte is a "reserved for future use".
response_length = 4

# data_length = length(data)
#println("data_length = $data_length")
send_max_length = 2 * data_length + 4
response_max_length = 2 * response_length + 4

#  Buffers.  Make sure all buffers are initialized to zero!
#  Stray special characters in uninitialized memory can wreak havoc!
send_packet = zeros(UInt8, send_max_length + response_max_length)
#  response_packet = zeros(UInt8, response_max_length)

#  Set up header
send_packet[1] = SOP
send_packet[2] = CHANNEL
send_packet[3] = 0x00

#  Build the packet, adding in special Packet characters as required.
#  Have to also track position in send_packet.
#  This was done with tricky pointer math in original C code.
j = 4  #  This is the index of send_packet. First 3 indices already filled.
for i in 1:data_length
#    println("j = $j")
    current_byte = data[i]
    #  Add EOP before the last byte
    if i == (data_length) #  This is the last loop.
        send_packet[j] = EOP
        j += 1
        send_packet[j] = current_byte  #  Last data byte.
    end
    if current_byte == SOP
        send_packet[j] = ESC
        j += 1
        send_packet[j] = xor_20(current_byte)
        j += 1
    elseif current_byte == EOP
        send_packet[j] = ESC
        j += 1
        send_packet[j] = xor_20(current_byte)
        j += 1
    elseif current_byte == CHANNEL
        send_packet[j] = ESC
        j += 1
        send_packet[j] = xor_20(current_byte)
        j += 1
    elseif current_byte == ESC
        send_packet[j] = ESC
        j += 1
        send_packet[j] = xor_20(current_byte)
        j += 1
    else
        send_packet[j] = current_byte
        j += 1
    end
end
#println("j = $j")  #  j should be length of packet

#    for i in 1:j
#    Printf.@printf("send_packet[%d] = %x\n", i, send_packet[i])
#    end
    length = j - 1 + response_max_length  # Compute the length after inserting special characters.
    return send_packet[1:length], length  # TEMPORARY Need to use response length!!!
end

 BYTESIDLECHAR = 0x4a
 BYTESESCCHAR  = 0x4d

function byte_to_core(send_data, data_length)

#data_length = length(send_data)
send_max_length = 2 * data_length + 4

send_packet = zeros(UInt8, send_max_length)
#send_packet = Array{UInt8,1}(undef, send_max_length)

j = 1
for i in 1:data_length
    current_byte = send_data[i]
    if current_byte == BYTESIDLECHAR
        send_packet[j] = BYTESESCCHAR
        j += 1
        send_packet[j] = xor_20(current_byte)
        j += 1
    elseif current_byte == BYTESESCCHAR
        send_packet[j] = BYTESESCCHAR
        j += 1
        send_packet[j] = xor_20(current_byte)
        j += 1
    else
        send_packet[j] = current_byte
        j += 1
    end
end
length = j - 1  #  Compute the length after inserting special characters.
return send_packet[1:length], length
end

function xor_20(val::UInt8)
    return xor(val,0x20)
end

function clean_response(response_data)
    #  The response packet will always be ? bytes long. Use 12 for now.
    response_packet = zeros(UInt8, 12)
    j = 1
    i = 1
    #  This loop cleans out the idles, which are removed, and removes escapes plus XORing the following byte.
    while j < 13
        current_byte = response_data[j]
        #  Must XOR packets which are preceded by certain special characters.
        if response_data[j] == BYTESIDLECHAR#  Skip and don't write anything to response packet.
            j += 1  #  Skip over the idle.
            continue
        elseif response_data[j] == BYTESESCCHAR
            response_packet[i] = xor(response_data[j + 1]) #  XOR the byte AFTER the escape!
            j += 2 #  jump escape and following byte.
        else
            response_packet[i] = current_byte
            j += 1
        end
        i += 1
    #    Printf.@printf("current_byte = %x\n", current_byte)
    end
    #  Print the response packet:
#    for i in 1:12
#Printf.@printf("response_packet[%d] = %x\n", i, response_packet[i])
#    end
    return response_packet #  This should return the entire packet plus 12 zeros which account for the response, which must be clocked out.
end

#  This function extracts the response data payload.
#  Note on the response data.  The first byte is the op-code with most signiificant bit
#  reversed.  For example, SEQUENTIAL_WRITE = 0x04 with MSB reversed is 0x84.
#  The 2nd byte is reserved; the 3rd and 4th bytes is the 16 bit length of the transaction.
function get_response_data(response_packet)
    data = zeros(UInt16, 12)
i = 1
j = 1
#  First, find SOP and increment i to next byte after.
while i < 12
    current_byte = response_packet[i]
    if response_packet[i] == SOP
        i += 1
        break #  Does this break out of while loop?
    end
    i += 1
end
while i < 12
    current_byte = response_packet[i]
    if response_packet[i] == SOP
        i += 1
        data[j] = xor_20(response_packet[i])
        j += 1
    elseif response_packet[i] == CHANNEL
     i += 2  # Skip special byte and following channel number.
     elseif response_packet[i] == EOP
     # Check to see if byte after EOP is an escape:
       if response_packet[i + 1]  == BYTESESCCHAR
         data[j] = xor_20(response_packet[i+2])
       else
         data[j] = response_packet[i + 1]  #  Byte after EOP.
         i = 13  #  Force break from while loop.
       end
     else
     data[j] = current_byte
     i += 1
     j += 1
 end
 end
# for i in 1:4
#Printf.@printf("data[%d] = %x\n", i, data[i])
# end

 bytes_written::UInt16 = (data[3] << 8) | data[4]
# Printf.@printf("Number of bytes written = %d\n", bytes_written)

 return bytes_written
end

function get_readback_data(avalon_packet)
    read_data = 0x00
    i = 1
    for i in 1:12
        if avalon_packet[i] == EOP
            read_data = avalon_packet[i+1]
            break
        end
    end
    return read_data
end

end # module

max_stars_repo_name: Greg-R/spiavalonfpga

max_issues_repo_name: Greg-R/spiavalonfpga

max_forks_repo_name: Greg-R/spiavalonfpga