text:
export AbConstraintSet, L1Constraint, L2Constraint, BasisL2Constraint, LInftyConstraint
export project, projectedge, normcs, scaletonorm, scaleproject

abstract type AbConstraintSet end

normcs(cs::CS, u) where CS<:AbConstraintSet = error("Not implemented for $CS")
project(cs::CS, u) where CS<:AbConstraintSet = error("Not implemented for $CS")
projectedge(cs::CS, u) where CS<:AbConstraintSet = error("Not implemented for $CS")
scaletonorm(cs::CS, u::Vector) where CS<:AbConstraintSet = u*cs.radius/normcs(cs, u)

###########################################################################  ℓ1
struct L1Constraint <: AbConstraintSet
  radius :: Real
  fromindex ::Int
end
L1Constraint() = L1Constraint(1., 2)
L1Constraint(r::Real) = L1Constraint(r, 2)

normcs(cs::L1Constraint, u::Vector) = norm(u[cs.fromindex:end], 1)

function project(cs::L1Constraint, u::Vector)
    uproj = projectl1ball(u[cs.fromindex:end], cs.radius)
    vcat(u[1:(cs.fromindex-1)], uproj)
end

function projectedge(cs::L1Constraint, u::Vector)
    uproj = projectl1ballegde(u[cs.fromindex:end], cs.radius)
    vcat(u[1:(cs.fromindex-1)], uproj)
end

function projectsimplex(u0, radius=1)
    u = sort(u0, rev=true)
    cssv = cumsum(u)
    rho = findlast(u.*(1:length(u)) .> (cssv-radius))
    theta = (cssv[rho]-radius)/rho
    max.(u0-theta, 0.)
end

projectl1balledge(u, radius=1) = sign.(u) .* projectsimplex(abs.(u), radius)
projectl1ball(u, radius=1) = sum(abs, u) <= radius ?  u : projectl1balledge(u, radius)

###########################################################################  ℓ2
struct L2Constraint <: AbConstraintSet
  radius :: Real
  fromindex :: Int
end
L2Constraint() = L2Constraint(1., 2)
L2Constraint(r::Real) = L2Constraint(r, 2)

normcs(cs::L2Constraint, u::Vector) = norm(u[cs.fromindex:end], 2)

function project(cs::L2Constraint, u::Vector)
    uproj = u[cs.fromindex:end] * min(1, cs.radius/norm(u[cs.fromindex:end]))
    vcat(u[1:(cs.fromindex-1)], uproj)
end

function projectedge(cs::L2Constraint, u::Vector)
    uproj = u[cs.fromindex:end] * cs.radius/norm(u[cs.fromindex:end])
    vcat(u[1:(cs.fromindex-1)], uproj)
end

###########################################################################  ℓ∞
struct LInftyConstraint <: AbConstraintSet
  radius :: Real
  fromindex :: Int
end
LInftyConstraint() = LInftyConstraint(1., 2)
LInftyConstraint(r::Real) = LInftyConstraint(r, 2)

normcs(cs::LInftyConstraint, u::Vector) = maximum(abs, u[cs.fromindex:end])

function project(cs::LInftyConstraint, u::Vector)
    uproj = clamp.(u[cs.fromindex:end], -cs.radius, cs.radius)
    vcat(u[1:(cs.fromindex-1)], uproj)
end

function projectedge(cs::LInftyConstraint, u::Vector)
    uproj = cs.radius * sign.(u[cs.fromindex:end])
    vcat(u[1:(cs.fromindex-1)], uproj)
end

#####################################################################  Basis ℓ2
struct BasisL2Constraint{B} <: AbConstraintSet
  basis::Type{B}
  X
  radius::Real
end
BasisL2Constraint(::Type{B}, X) where B = BasisL2Constraint{B}(B, X, 1.)

normcs(cs::BasisL2Constraint{B}, u::Vector) where B = norm(calc(B, cs.X, u))

_scaler(::Type{R}, x) where {R<:Relu} = x
_scaler(::Type{R}, x) where {R<:Requ} = √x

function project(cs::BasisL2Constraint{B}, u::Vector) where {B<:AbUniDirBasis}
  u * _scaler(B, min(1, cs.radius/normcs(cs, u)))
end

function projectedge(cs::BasisL2Constraint{B}, u::Vector) where {B<:AbUniDirBasis}
  u * _scaler(B, cs.radius/normcs(cs, u))
end

function scaleproject(cs::BasisL2Constraint{B}, u::Vector) where {B<:AbUniDirBasis}
    scalefactor = cs.radius/normcs(cs, u) #min(1, cs.radius/normcs(cs, u))
    1/scalefactor, u*_scaler(B, scalefactor)
end

#########################################################################  ℓend


text:
#=== Simplex ===#
struct Simplex{N, P} <: AbstractSimplex
    vs::NTuple{N, P}
    hash::UInt64
    Simplex{N,P}(vals::NTuple{N, P}) where {N, P} = new(vals, hash(Set(vals)))
    Simplex{0,P}() where {P} = new((), zero(UInt))
end
Simplex(splx::P...) where {P} = Simplex{length(splx),P}(splx)
Simplex(splx::Vector{P}) where {P} = Simplex(splx...)
Simplex(::Type{P}) where {P} = Simplex{0,P}()

# Private methods

show(io::IO, splx::Simplex) = show(io, "σ$(splx.vs)")
eltype(splx::Simplex{N,P}) where {N,P} = P
eltype(::Type{Simplex{N,P}}) where {N,P} = P

# Public methods

dim(splx::Simplex{N,P}) where {N,P} = N-1
hash(splx::Simplex) = splx.hash
==(a::Simplex, b::Simplex) = a.hash == b.hash

function vertices(splx::Simplex)
    vtxs = typeof(hash(splx))[]
    for v in splx.vs
        push!(vtxs, hash(Simplex(v)))
    end
    return vtxs
end

union(u::Simplex, v::Simplex) = Simplex(collect(u.vs ∪ v.vs))

# TODO make generated
function faces(σ::Simplex)
    d = dim(σ)
    d == 0 && return Simplex[]
    (Simplex([ifelse(j < i, σ.vs[j], σ.vs[j+1]) for j in 1:d]) for i in 1:(d+1))
end

function boundary(::Type{R}, σ::Simplex) where {R}
    d = dim(σ)
    ch = Chain(d-1, UInt, R)
    d == 0 && return ch
    o = one(R)
    for (i,face) in enumerate(faces(σ))
        push!(ch, hash(face)=>( isodd(i) ? o : -o))
    end
    return ch
end

# Public methods for AbstractSimplex

indices(splx::Simplex) = collect(splx.vs)

# Misc. methods

function volume(S::AbstractMatrix)
    d, vc = size(S)
    @assert d == vc-1 "Number of vertices in simplex should be dim+1"
    v0 = S[:,1]
    return abs(det(S[:,2:end] .- v0))/prod(1:d)
end
volume(s::Simplex{N,Int}, X::AbstractMatrix) where {N}= volume(X[:, indices(s)])

parse(::Type{Simplex{N,P}}, str::AbstractString) where {N,P} =
    Simplex(map(e->parse(P, e), split(str, ' ')))

ext: jl


text:
include("generic_functions.jl")
#example
output_instance_path = "..\\data_ROPF"
output_decomposition_path = "..\\data_sdp"
list_generation = ["generation"]
generation_files_path = "..\\data_ROPF"
network_name = "118nodes"

networks = ["118nodes_disconnect"]

for network_name in networks
    f = open("results_test_allshunts1_$(network_name).csv", "w")
    write(f, "Numero ; UB ; LB ; Gap (%)\n")
    for i in 0:2499
        # network_name, numero = split(instance_name, '-')
        numero = "test_$i"
        instance_name = network_name*"-"*numero
        matpower_instance_path = "..\\mp_data\\$(network_name)\\$(numero).mat"
        for generation in list_generation
            construct_dat_file_ROPF(instance_name, matpower_instance_path, output_instance_path)
            # generate_clique_decomposition(instance_name, matpower_instance_path, output_decomposition_path)
            ROPF = ROPF_infos(instance_name,
            matpower_instance_path,
            output_instance_path,
            "cholesky",
            output_decomposition_path,
            generation,
            generation_files_path)
            UB = Inf
            UB, LB = solve1_refgen_allshunts1(ROPF)
            println("UB=$UB ; LB=$LB")
            if UB != Inf
                cp("solutions\\solve1_allshunts1_refgen_solution_$(instance_name).csv", "solutions_mp_data\\allshunts1_refgen_solution_$(instance_name).csv")
            end
            gap = (UB-LB)/UB
            write(f, "$numero ; $UB ; $LB ; $gap \n")
        end
    end
    close(f)
end


text:
module PooledArrays

import DataAPI
import Future.copy!

export PooledArray, PooledVector, PooledMatrix

##############################################################################
##
## PooledArray type definition
##
##############################################################################

const DEFAULT_POOLED_REF_TYPE = UInt32
const DEFAULT_SIGNED_REF_TYPE = Int32

# This is used as a wrapper during PooledArray construction only, to distinguish
# arrays of pool indices from normal arrays
mutable struct RefArray{R}
    a::R
end

function _invert(d::Dict{K,V}) where {K,V}
    d1 = Vector{K}(undef, length(d))
    for (k, v) in d
        d1[v] = k
    end
    return d1
end

mutable struct PooledArray{T, R<:Integer, N, RA} <: AbstractArray{T, N}
    refs::RA
    pool::Vector{T}
    invpool::Dict{T,R}
    # refcount[] is 1 if only one PooledArray holds a reference to pool and invpool
    refcount::Threads.Atomic{Int}

    function PooledArray{T,R,N,RA}(rs::RefArray{RA}, invpool::Dict{T, R},
                                   pool::Vector{T}=_invert(invpool),
                                   refcount::Threads.Atomic{Int}=Threads.Atomic{Int}(1)) where {T,R,N,RA<:AbstractArray{R, N}}
        # we currently support only 1-based indexing for refs
        # TODO: change to Base.require_one_based_indexing after we drop Julia 1.0 support
        for ax in axes(rs.a)
            if first(ax) != 1
                throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
            end
        end

        # this is a quick but incomplete consistency check
        if length(pool) != length(invpool)
            throw(ArgumentError("inconsistent pool and invpool"))
        end
        if length(rs.a) > 0
            # 0 indicates #undef
            # refs mustn't overflow pool
            minref, maxref = extrema(rs.a)
            if (minref < 0 || maxref > length(invpool))
                throw(ArgumentError("Reference array points beyond the end of the pool"))
            end
        end
        pa = new{T,R,N,RA}(rs.a, pool, invpool, refcount)
        finalizer(x -> Threads.atomic_sub!(x.refcount, 1), pa)
        return pa
    end
end
const PooledVector{T,R} = PooledArray{T,R,1}
const PooledMatrix{T,R} = PooledArray{T,R,2}

const PooledArrOrSub = Union{SubArray{T, N, <:PooledArray{T, R}},
                             PooledArray{T, R, N}} where {T, N, R}

##############################################################################
##
## PooledArray constructors
##
# Algorithm:
# * Start with:
#   * A null pool
#   * A pre-allocated refs
#   * A hash from T to Int
# * Iterate over d
#   * If value of d in pool already, set the refs accordingly
#   * If value is new, add it to the pool, then set refs
##############################################################################

# Echo inner constructor as an outer constructor
@inline PooledArray(refs::RefArray{RA}, invpool::Dict{T,R}, pool::Vector{T}=_invert(invpool),
            refcount::Threads.Atomic{Int}=Threads.Atomic{Int}(1)) where {T,R,RA<:AbstractArray{R}} =
    PooledArray{T,R,ndims(RA),RA}(refs, invpool, pool, refcount)

# workaround https://github.com/JuliaLang/julia/pull/39809
_our_copy(x) = copy(x)

function _our_copy(x::SubArray{<:Any, 0})
    y = similar(x)
    y[] = x[]
    return y
end

@inline function PooledArray(d::PooledArrOrSub)
    Threads.atomic_add!(refcount(d), 1)
    return PooledArray(RefArray(_our_copy(DataAPI.refarray(d))),
                       DataAPI.invrefpool(d), DataAPI.refpool(d), refcount(d))
end

function _label(xs::AbstractArray,
                ::Type{T}=eltype(xs),
                ::Type{I}=DEFAULT_POOLED_REF_TYPE,
                start = 1,
                labels = Array{I}(undef, size(xs)),
                invpool::Dict{T,I} = Dict{T, I}(),
                pool::Vector{T} = T[],
                nlabels = 0,
               ) where {T, I<:Integer}

    @inbounds for i in start:length(xs)
        idx = i + firstindex(xs) - 1
        if !isassigned(xs, idx)
            labels[i] = zero(I)
        else
            x = xs[idx]
            lbl = get(invpool, x, zero(I))
            if lbl !== zero(I)
                labels[i] = lbl
            else
                if nlabels == typemax(I)
                    I2 = _widen(I)
                    return _label(xs, T, I2, i, convert(Vector{I2}, labels),
                                convert(Dict{T, I2}, invpool), pool, nlabels)
                end
                nlabels += 1
                labels[i] = nlabels
                invpool[x] = nlabels
                push!(pool, x)
            end
        end
    end
    labels, invpool, pool
end

_widen(::Type{UInt8}) = UInt16
_widen(::Type{UInt16}) = UInt32
_widen(::Type{UInt32}) = UInt64
_widen(::Type{Int8}) = Int16
_widen(::Type{Int16}) = Int32
_widen(::Type{Int32}) = Int64

# Constructor from array, invpool, and ref type

"""
    PooledArray(array, [reftype]; signed=false, compress=false)

Freshly allocate `PooledArray` using the given array as a source where each
element will be referenced as an integer of the given type.

If `reftype` is not specified then `PooledArray` constructor is not type stable.
In this case Boolean keyword arguments `signed` and `compress`
determine the type of integer references. By default (`signed=false`), unsigned integers
are used, as they have a greater range.
However, the Arrow standard at https://arrow.apache.org/, as implemented in
the Arrow package, requires signed integer types, which are provided when `signed=true`.
When `compress=false`, `reftype` is a 32-bits type (`UInt32` for unsigned, `Int32` for signed);
when `compress=true`, `reftype` is chosen to be the smallest integer type that is
large enough to hold all unique values in `array`.

Note that if you hold mutable objects in `PooledArray` it is not allowed to modify them
after they are stored in it.

In order to improve performance of `getindex` and `copyto!` operations `PooledArray`s
may share pools. This sharing is automatically undone by copying a shared pool before
adding new values to it.

It is not safe to assign values that are not already present in a `PooledArray`'s pool
from one thread while either reading or writing to the same array from another thread
(even if pools are not shared). However, reading and writing from different threads is safe
if all values already exist in the pool.
"""
PooledArray

@inline function PooledArray{T}(d::AbstractArray, r::Type{R}) where {T,R<:Integer}
    refs, invpool, pool = _label(d, T, R)

    if length(invpool) > typemax(R)
        throw(ArgumentError("Cannot construct a PooledArray with type $R with a pool of size $(length(pool))"))
    end

    # Assertions are needed since _label is not type stable
    return PooledArray(RefArray(refs::Array{R, ndims(d)}), invpool::Dict{T,R}, pool)
end

@inline function PooledArray{T}(d::AbstractArray; signed::Bool=false, compress::Bool=false) where {T}
    R = signed ? (compress ? Int8 : DEFAULT_SIGNED_REF_TYPE) : (compress ? UInt8 : DEFAULT_POOLED_REF_TYPE)
    refs, invpool, pool = _label(d, T, R)
    return PooledArray(RefArray(refs), invpool, pool)
end

@inline PooledArray(d::AbstractArray{T}, r::Type) where {T} = PooledArray{T}(d, r)
@inline PooledArray(d::AbstractArray{T}; signed::Bool=false, compress::Bool=false) where {T} =
    PooledArray{T}(d, signed=signed, compress=compress)

# Construct an empty PooledVector of a specific type
@inline PooledArray(t::Type) = PooledArray(Array(t,0))
@inline PooledArray(t::Type, r::Type) = PooledArray(Array(t,0), r)

##############################################################################
##
## Basic interface functions
##
##############################################################################

DataAPI.refarray(pa::PooledArray) = pa.refs
DataAPI.refvalue(pa::PooledArray, i::Integer) = pa.pool[i]
DataAPI.refpool(pa::PooledArray) = pa.pool
DataAPI.invrefpool(pa::PooledArray) = pa.invpool
refcount(pa::PooledArray) = pa.refcount

DataAPI.refarray(pav::SubArray{<:Any, <:Any, <:PooledArray}) = view(parent(pav).refs, pav.indices...)
DataAPI.refvalue(pav::SubArray{<:Any, <:Any, <:PooledArray}, i::Integer) = parent(pav).pool[i]
DataAPI.refpool(pav::SubArray{<:Any, <:Any, <:PooledArray}) = parent(pav).pool
DataAPI.invrefpool(pav::SubArray{<:Any, <:Any, <:PooledArray}) = parent(pav).invpool
refcount(pav::SubArray{<:Any, <:Any, <:PooledArray}) = parent(pav).refcount

Base.size(pa::PooledArray) = size(pa.refs)
Base.length(pa::PooledArray) = length(pa.refs)
Base.lastindex(pa::PooledArray) = lastindex(pa.refs)

Base.copy(pa::PooledArrOrSub) = PooledArray(pa)

# here we do not allow dest to be SubArray as copy! is intended to replace whole arrays
# slow path will be used for SubArray
function copy!(dest::PooledArray{T, R, N},
               src::PooledArrOrSub{T, N, R}) where {T, N, R}
    copy!(dest.refs, DataAPI.refarray(src))
    src_refcount = refcount(src)

    if dest.pool !== DataAPI.refpool(src)
        Threads.atomic_sub!(dest.refcount, 1)
        Threads.atomic_add!(src_refcount, 1)
        dest.pool = DataAPI.refpool(src)
        dest.invpool = DataAPI.invrefpool(src)
        dest.refcount = src_refcount
    else
        @assert dest.invpool === DataAPI.invrefpool(src)
        @assert dest.refcount === src_refcount
    end
    return dest
end

# this is needed as Julia Base uses a special path for this case we want to avoid
Base.copyto!(dest::PooledArrOrSub{T, N, R}, src::PooledArrOrSub{T, N, R}) where {T, N, R} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::PooledArrOrSub{T, N, R}, doffs::Union{Signed, Unsigned},
                      src::PooledArrOrSub{T, N, R}, soffs::Union{Signed, Unsigned},
                      n::Union{Signed, Unsigned}) where {T, N, R}
    n == 0 && return dest
    n > 0 || Base._throw_argerror()
    if soffs < 1 || doffs < 1 || soffs + n - 1 > length(src) || doffs + n - 1 > length(dest)
        throw(BoundsError())
    end

    dest_pa = dest isa PooledArray ? dest : parent(dest)
    src_refcount = refcount(src)

    # if dest_pa.pool is empty we can safely replace it as we are sure it holds
    # no information; having this path is useful because then we can efficiently
    # `copyto!` into a fresh `PooledArray` created using the `similar` function
    if DataAPI.refpool(dest) === DataAPI.refpool(src)
        @assert DataAPI.invrefpool(dest) === DataAPI.invrefpool(src)
        @assert refcount(dest) === refcount(src)
        copyto!(DataAPI.refarray(dest), doffs, DataAPI.refarray(src), soffs, n)
    elseif length(dest_pa.pool) == 0
        @assert length(dest_pa.invpool) == 0
        Threads.atomic_add!(src_refcount, 1)
        dest_pa.pool = DataAPI.refpool(src)
        dest_pa.invpool = DataAPI.invrefpool(src)
        Threads.atomic_sub!(dest_pa.refcount, 1)
        dest_pa.refcount = src_refcount
        copyto!(DataAPI.refarray(dest), doffs, DataAPI.refarray(src), soffs, n)
    else
        @inbounds for i in 0:n-1
            dest[doffs+i] = src[soffs+i]
        end
    end
    return dest
end

function Base.resize!(pa::PooledArray{T,R,1}, n::Integer) where {T,R}
    oldn = length(pa.refs)
    resize!(pa.refs, n)
    pa.refs[oldn+1:n] .= zero(R)
    return pa
end

function Base.reverse(x::PooledArray)
    Threads.atomic_add!(x.refcount, 1)
    PooledArray(RefArray(reverse(x.refs)), x.invpool, x.pool, x.refcount)
end

function Base.permute!!(x::PooledArray, p::AbstractVector{T}) where T<:Integer
    Base.permute!!(x.refs, p)
    return x
end

function Base.invpermute!!(x::PooledArray, p::AbstractVector{T}) where T<:Integer
    Base.invpermute!!(x.refs, p)
    return x
end

Base.similar(pa::PooledArray{T,R}, S::Type, dims::Dims) where {T,R} =
    PooledArray(RefArray(zeros(R, dims)), Dict{S,R}())

Base.findall(pdv::PooledVector{Bool}) = findall(convert(Vector{Bool}, pdv))

##############################################################################
##
## map
## Calls `f` only once per pool entry.
##
##############################################################################

"""
    map(f, x::PooledArray; pure::Bool=false)

Transform `PooledArray` `x` by applying `f` to each element.

If `pure=true` then `f` is applied to each element of pool of `x`
exactly once (even if some elements in pool are not present it `x`).
This will typically be much faster when the proportion of unique values
in `x` is small.

If `pure=false`, the returned array will use the same reference type
as `x`, or `Int` if the number of unique values in the result is too large
to fit in that type.
"""
function Base.map(f, x::PooledArray{<:Any, R, N, RA}; pure::Bool=false)::Union{PooledArray{<:Any, R, N, RA},
                                                                               PooledArray{<:Any, Int, N,
                                                                                typeof(similar(x.refs, Int, ntuple(i -> 0, ndims(x.refs))))}} where {R, N, RA}
    pure && return _map_pure(f, x)
    length(x) == 0 && return PooledArray([f(v) for v in x])
    v1 = f(x[1])
    invpool = Dict(v1 => one(eltype(x.refs)))
    pool = [v1]
    labels = similar(x.refs)
    labels[1] = 1
    nlabels = 1
    return _map_notpure(f, x, 2, invpool, pool, labels, nlabels)
end

function _map_notpure(f, xs::PooledArray, start,
                      invpool::Dict{T,I}, pool::Vector{T},
                      labels::AbstractArray{I}, nlabels::Int) where {T, I<:Integer}
    for i in start:length(xs)
        vi = f(xs[i])
        lbl = get(invpool, vi, zero(I))
        if lbl != zero(I)
            labels[i] = lbl
        else
            if nlabels == typemax(I) || !(vi isa T)
                I2 = nlabels == typemax(I) ? Int : I
                T2 = vi isa T ? T : Base.promote_typejoin(T, typeof(vi))
                nlabels += 1
                invpool2 = convert(Dict{T2, I2}, invpool)
                invpool2[vi] = nlabels
                pool2 = convert(Vector{T2}, pool)
                push!(pool2, vi)
                labels2 = convert(AbstractArray{I2}, labels)
                labels2[i] = nlabels
                return _map_notpure(f, xs, i + 1, invpool2, pool2,
                                    labels2, nlabels)
            end
            nlabels += 1
            labels[i] = nlabels
            invpool[vi] = nlabels
            push!(pool, vi)
        end
    end
    return PooledArray(RefArray(labels), invpool, pool)
end

function _map_pure(f, x::PooledArray)
    ks = collect(keys(x.invpool))
    vs = collect(values(x.invpool))
    ks1 = map(f, ks)
    uks = Set(ks1)
    if length(uks) < length(ks1)
        # this means some keys have repeated
        newinvpool = Dict{eltype(ks1), eltype(vs)}()
        translate = Dict{eltype(vs), eltype(vs)}()
        i = 1
        for (k, k1) in zip(ks, ks1)
            if haskey(newinvpool, k1)
                translate[x.invpool[k]] = newinvpool[k1]
            else
                newinvpool[k1] = i
                translate[x.invpool[k]] = i
                i+=1
            end
        end
        refarray = map(x->translate[x], x.refs)
    else
        newinvpool = Dict(zip(ks1, vs))
        refarray = copy(x.refs)
    end
    return PooledArray(RefArray(refarray), newinvpool)
end

##############################################################################
##
## Sorting can use the pool to speed things up
##
##############################################################################

function groupsort_indexer(x::AbstractVector, ngroups::Integer, perm)
    # translated from Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).

    # count group sizes, location 0 for NA
    n = length(x)
    # counts = x.invpool
    counts = fill(0, ngroups + 1)
    @inbounds for i = 1:n
        counts[x[i] + 1] += 1
    end
    counts[2:end] = counts[perm.+1]

    # mark the start of each contiguous group of like-indexed data
    where = fill(1, ngroups + 1)
    @inbounds for i = 2:ngroups+1
        where[i] = where[i - 1] + counts[i - 1]
    end

    # this is our indexer
    result = fill(0, n)
    iperm = invperm(perm)

    @inbounds for i = 1:n
        label = iperm[x[i]] + 1
        result[where[label]] = i
        where[label] += 1
    end
    result, where, counts
end

function Base.sortperm(pa::PooledArray; alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE,
                       lt::Function=isless, by::Function=identity,
                       rev::Bool=false, order=Base.Sort.Forward,
                       _ord = Base.ord(lt, by, rev, order),
                       poolperm = sortperm(pa.pool, alg=alg, order=_ord))

    groupsort_indexer(pa.refs, length(pa.pool), poolperm)[1]
end

Base.sort(pa::PooledArray; kw...) = pa[sortperm(pa; kw...)]

#type FastPerm{O<:Base.Sort.Ordering,V<:AbstractVector} <: Base.Sort.Ordering
#    ord::O
#    vec::V
#end
#Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ForwardOrdering,V}) = x[sortperm(o.vec)]
#Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ReverseOrdering,V}) = x[reverse(sortperm(o.vec))]
#Perm{O<:Base.Sort.Ordering}(o::O, v::PooledVector) = FastPerm(o, v)

##############################################################################
##
## conversions
##
##############################################################################

function Base.convert(::Type{PooledArray{S,R1,N}}, pa::PooledArray{T,R2,N}) where {S,T,R1<:Integer,R2<:Integer,N}
    invpool_conv = convert(Dict{S,R1}, pa.invpool)
    @assert invpool_conv !== pa.invpool

    if R1 === R2
        refs_conv = pa.refs
    else
        refs_conv = convert(Array{R1,N}, pa.refs)
        @assert refs_conv !== pa.refs
    end

    return PooledArray(RefArray(refs_conv), invpool_conv)
end

Base.convert(::Type{PooledArray{T,R,N}}, pa::PooledArray{T,R,N}) where {T,R<:Integer,N} = pa
Base.convert(::Type{PooledArray{S,R1}}, pa::PooledArray{T,R2,N}) where {S,T,R1<:Integer,R2<:Integer,N} =
    convert(PooledArray{S,R1,N}, pa)
Base.convert(::Type{PooledArray{S}}, pa::PooledArray{T,R,N}) where {S,T,R<:Integer,N} =
    convert(PooledArray{S,R,N}, pa)
Base.convert(::Type{PooledArray}, pa::PooledArray{T,R,N}) where {T,R<:Integer,N} = pa

Base.convert(::Type{PooledArray{S,R,N}}, a::AbstractArray{T,N}) where {S,T,R<:Integer,N} =
    PooledArray(convert(Array{S,N}, a), R)
Base.convert(::Type{PooledArray{S,R}}, a::AbstractArray{T,N}) where {S,T,R<:Integer,N} =
    PooledArray(convert(Array{S,N}, a), R)
Base.convert(::Type{PooledArray{S}}, a::AbstractArray{T,N}) where {S,T,N} =
    PooledArray(convert(Array{S,N}, a))
Base.convert(::Type{PooledArray}, a::AbstractArray) =
    PooledArray(a)

function Base.convert(::Type{Array{S, N}}, pa::PooledArray{T, R, N}) where {S, T, R, N}
    res = Array{S}(undef, size(pa))
    for i in 1:length(pa)
        if pa.refs[i] != 0
            res[i] = pa.pool[pa.refs[i]]
        end
    end
    return res
end

Base.convert(::Type{Vector}, pv::PooledVector{T, R}) where {T, R} = convert(Array{T, 1}, pv)

Base.convert(::Type{Matrix}, pm::PooledMatrix{T, R}) where {T, R} = convert(Array{T, 2}, pm)

Base.convert(::Type{Array}, pa::PooledArray{T, R, N}) where {T, R, N} = convert(Array{T, N}, pa)

##############################################################################
##
## indexing
##
##############################################################################

# We need separate functions due to dispatch ambiguities

Base.@propagate_inbounds function Base.getindex(A::PooledArray, I::Int)
    idx = DataAPI.refarray(A)[I]
    iszero(idx) && throw(UndefRefError())
    return @inbounds DataAPI.refpool(A)[idx]
end

# we handle fast only the case when the first index is an abstract vector
# this is to make sure other indexing synraxes use standard dispatch from Base
# the reason is that creation of DataAPI.refarray(A) is unfortunately slow
Base.@propagate_inbounds function Base.getindex(A::PooledArrOrSub,
                                                I1::AbstractVector,
                                                I2::Union{Real, AbstractVector}...)
    # make sure we do not increase A.refcount in case creation of newrefs fails
    newrefs = DataAPI.refarray(A)[I1, I2...]
    @assert newrefs isa AbstractArray
    Threads.atomic_add!(refcount(A), 1)
    return PooledArray(RefArray(newrefs), DataAPI.invrefpool(A), DataAPI.refpool(A), refcount(A))
end

Base.@propagate_inbounds function Base.isassigned(pa::PooledArrOrSub, I::Int...)
    !iszero(DataAPI.refarray(pa)[I...])
end

##############################################################################
##
## setindex!() definitions
##
##############################################################################

function getpoolidx(pa::PooledArray{T,R}, val::Any) where {T,R}
    val::T = convert(T,val)
    pool_idx = get(pa.invpool, val, zero(R))
    if pool_idx == zero(R)
        pool_idx = unsafe_pool_push!(pa, val)
    end
    return pool_idx
end

function unsafe_pool_push!(pa::PooledArray{T,R}, val) where {T,R}
    # Warning - unsafe_pool_push! may not be used in any multithreaded context
    _pool_idx = length(pa.pool) + 1
    if _pool_idx > typemax(R)
        throw(ErrorException(string(
            "You're using a PooledArray with ref type $R, which can only hold $(Int(typemax(R))) values,\n",
            "and you just tried to add the $(typemax(R)+1)th reference.  Please change the ref type\n",
            "to a larger int type, or use the default ref type ($DEFAULT_POOLED_REF_TYPE)."
           )))
    end
    pool_idx = convert(R, _pool_idx)
    if pa.refcount[] > 1
        pa.invpool = copy(pa.invpool)
        pa.pool = copy(pa.pool)
        Threads.atomic_sub!(pa.refcount, 1)
        pa.refcount = Threads.Atomic{Int}(1)
    end
    pa.invpool[val] = pool_idx
    push!(pa.pool, val)
    pool_idx
end

# assume PooledArray is only used with Arrays as this is what _label does
# this simplifies code below
Base.IndexStyle(::Type{<:PooledArray}) = IndexLinear()

Base.@propagate_inbounds function Base.setindex!(x::PooledArray, val, ind::Int)
    x.refs[ind] = getpoolidx(x, val)
    return x
end

##############################################################################
##
## growing and shrinking
##
##############################################################################

function Base.push!(pv::PooledVector, v) # this function is not thread safe
    push!(pv.refs, getpoolidx(pv, v))
    return pv
end

function Base.insert!(pv::PooledVector, i::Integer, v) # this function is not thread safe
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    if !(1 <= i <= length(pv.refs) + 1)
        throw(BoundsError("attempt to insert to a vector with length $(length(pv)) at index $i"))
    end
    insert!(pv.refs, i, getpoolidx(pv, v))
    return pv
end

function Base.append!(pv::PooledVector, items::AbstractArray)
    itemindices = eachindex(items)
    l = length(pv)
    n = length(itemindices)
    resize!(pv.refs, l+n)
    copyto!(pv, l+1, items, first(itemindices), n)
    return pv
end

Base.pop!(pv::PooledVector) = pv.pool[pop!(pv.refs)]

function Base.pushfirst!(pv::PooledVector{S,R}, v::T) where {S,R,T}
    pushfirst!(pv.refs, getpoolidx(pv, v))
    return pv
end

Base.popfirst!(pv::PooledVector) = pv.pool[popfirst!(pv.refs)]

Base.empty!(pv::PooledVector) = (empty!(pv.refs); pv)

Base.deleteat!(pv::PooledVector, inds) = (deleteat!(pv.refs, inds); pv)

function _vcat!(c, a, b)
    copyto!(c, 1, a, 1, length(a))
    return copyto!(c, length(a)+1, b, 1, length(b))
end

function Base.vcat(a::PooledArray{<:Any, <:Integer, 1}, b::AbstractArray{<:Any, 1})
    output = similar(b, promote_type(eltype(a), eltype(b)), length(b) + length(a))
    return _vcat!(output, a, b)
end

function Base.vcat(a::AbstractArray{<:Any, 1}, b::PooledArray{<:Any, <:Integer, 1})
    output = similar(a, promote_type(eltype(a), eltype(b)), length(b) + length(a))
    return _vcat!(output, a, b)
end

function Base.vcat(a::PooledArray{T, <:Integer, 1}, b::PooledArray{S, <:Integer, 1}) where {T, S}
    ap = a.invpool
    bp = b.invpool

    U = promote_type(T,S)

    poolmap = Dict{Int, Int}()
    l = length(ap)
    newlabels = Dict{U, Int}(ap)
    for (x, i) in bp
        j = if x in keys(ap)
            poolmap[i] = ap[x]
        else
            poolmap[i] = (l+=1)
        end
        newlabels[x] = j
    end
    types = [UInt8, UInt16, UInt32, UInt64]
    tidx = findfirst(t->l < typemax(t), types)
    refT = types[tidx]
    refs2 = map(r->convert(refT, poolmap[r]), b.refs)
    newrefs = Base.typed_vcat(refT, a.refs, refs2)
    return PooledArray(RefArray(newrefs), convert(Dict{U, refT}, newlabels))
end

fast_sortable(y::PooledArray) = _fast_sortable(y)
fast_sortable(y::PooledArray{T}) where {T<:Integer} = isbitstype(T) ? y : _fast_sortable(y)

function _fast_sortable(y::PooledArray)
    poolranks = invperm(sortperm(y.pool))
    newpool = Dict(j=>convert(eltype(y.refs), i) for (i,j) in enumerate(poolranks))
    PooledArray(RefArray(y.refs), newpool)
end

_perm(o::F, z::V) where {F, V} = Base.Order.Perm{F, V}(o, z)

Base.Order.Perm(o::Base.Order.ForwardOrdering, y::PooledArray) = _perm(o, fast_sortable(y))

function Base.repeat(x::PooledArray, m::Integer...)
    Threads.atomic_add!(x.refcount, 1)
    PooledArray(RefArray(repeat(x.refs, m...)), x.invpool, x.pool, x.refcount)
end

function Base.repeat(x::PooledArray; inner = nothing, outer = nothing)
    Threads.atomic_add!(x.refcount, 1)
    PooledArray(RefArray(repeat(x.refs; inner = inner, outer = outer)),
                                x.invpool, x.pool, x.refcount)
end

end


text:
__precompile__()
module Utils

export convert_tmpfile, defaultinputdict!, printparams, initlogo

using SeisIO, Printf, Dates, JLD2, FileIO

"""
convert_tmpfile(InputDict::Dict)

convert temporal file in "./seisdownload_tmp" to prescribed format.
It has salvage mode, which allows to compile the temporal files in the case of failing during the download.
"""
function convert_tmpfile(InputDict::Dict; salvage::Bool=false)

    println("-------START CONVERTING-------------")

	# save data to fopath file
	fopath = InputDict["fopath"]

	t = jldopen(InputDict["finame"])

	file = jldopen(fopath, "w")

	#!!!should be debuged because of Isocomponents!!!#
	file["info/stationlist"]     = t["info/stationlist"];

	if InputDict["IsStartendtime"]
		file["info/DLtimestamplist"] = InputDict["DLtimestamplist_selected"];
		file["info/starttime"]       = InputDict["starttime"];
		file["info/endtime"]         = InputDict["endtime"];
	else
		file["info/DLtimestamplist"] = t["info/DLtimestamplist"];
		file["info/starttime"]       = t["info/starttime"];
		file["info/endtime"]         = t["info/endtime"];
	end

	JLD2.close(t)

	# find all temporal files
    paths_all = ls(InputDict["tmppath"])
    fmt = InputDict["outputformat"]

    stationlist     = []
    DLtimestamplist = []
    varnamelist     = []

	if InputDict["IsIsolateComponents"]
		# isolate components based on priority dictionary
		paths = isolate_components(paths_all, InputDict)
	else
		paths = paths_all
	end

    for path in paths
        #println(path)

		# split path
		# tmpname = split(path, "/")[end]

		# ftmpname = split(tmpname, ".")
		# if occursin("-", ftmpname[3])
		# 	# format would be y, jd, T00:00:00, sta, loc, cha
		# 	y, d, tmpT, net, sta, loc, cha = split(tmpname, ".")
		# 	#iso_stationinfo = (join([y, d, net, sta, loc], "-"), cha)
		# 	iso_stationinfo = (join([y, d, net, sta], "-"), cha)
		#
		# elseif !occursin("-", ftmpname[3]) && IsIsolateComponents == true
		# 	@warn "format of tmp file is not y, jd, time, sta, loc, cha. Turn off IsIsolateComponents."
		# 	IsIsolateComponents = false
		# 	iso_stationinfo = []
		# end

		#println(iso_stationinfo)

        S = try
				SeisIO.rseis(path)[1]
			catch y
				println(y)
				@warn("cannot read tmpfile in seisremoveeq_tmp_sample. skipping")
				continue;
			end
        #println(S)

        for ii = 1:S.n #loop at each seis channel

			# if IsIsolateComponents
			# 	# check channel isolation: if y, d, net, sta and components are same, we discard one of them.
			# 	@show isostationlist
			# 	@show iso_stationinfo
			# 	conflictsta = findfirst(x -> x[1]==iso_stationinfo[1] && string(x[2][end])==string(iso_stationinfo[2][end]), isostationlist)
			# else
			# 	conflictsta = []
			# end

			#println(conflictsta)
			# if !isempty(conflictsta)
			# 	# here this channel has conflicting channel such as same day, same components but different channel.
			# 	# if IsIsolateComponents == true, skip to avoid multiple stations for the purpose of cross-correlation.
			# 	if IsIsolateComponents
			# 		sta1 = join([iso_stationinfo[1], iso_stationinfo[2]], "-")
			# 		sta2 = join([conflictsta[1][1], conflictsta[1][2]], "-")
			# 		txt = @sprintf("Two identical location but different channel were found.\n%s :%s: Discard %s.",
			# 					sta1, sta2, sta1)
			# 		println(txt)
			# 		#rm(path)
			# 		continue;
			# 	end
			# end

			# if IsIsolateComponents
			# 	# update isostationlist
			# 	push!(isostationlist, iso_stationinfo)
			# end

            # make station list
            staid = S[ii].id

            # save data (checking whether it's already in the jld2 because it causes an error)
            #parse info
            s_str = string(u2d(S[ii].t[1,2]*1e-6))

            # select output format
            if fmt == "JLD2"
                yj = parse(Int64, s_str[1:4])
                mj = parse(Int64, s_str[6:7])
                dj = parse(Int64, s_str[9:10])
                tj = string(s_str)[11:19]

                djm2j = md2j(yj, mj, dj)
                groupname = string(yj)*"."*string(djm2j)*"."*tj #Year_Julianday_Starttime
                varname = joinpath(groupname, staid)

                if isempty(filter(x -> x==varname, varnamelist))
                    push!(varnamelist, varname)
                    file[varname] = S[ii]
                end

            else
                error("output format in $fmt is not implemented yet.")
            end
        end

		# remove tmpfile
		rm(path)

    end

    JLD2.close(file)

	rm(InputDict["tmppath"], recursive=true, force=true)

    return nothing
end

"""
isolate_components(paths_all::AbstractArray, InputDict::Dict)

Isolate components based on InputDict["priority_channles"].

I. What is the purpose?
	Some stations have multiple stations with different channels
	(e.g. BP.LCCB..BP1 and BP.LCCB.40.SP1 during transition.)
	This increases number of xcorr processes, while the phisicall
	meaning is auto correlation using different stations; thus we duplicate
	auto-correlation in this case (e.g. BP1-BP1 and SP1-BP1 at same place).
	To avoid that, we implemented this function to isolate components with
	each stations based on priority, so that we can perform long-term dv/v
	analysis without such duplications.

II. Process flow

	1. Try to find channels which has same network.station and same component,
	but different channel

	2. If found, search priority with `haskey(InputDict["priority_channles"])`.
	If not found, arbitrually pick up the first one.

	3. If the second one has priority, search the order of it in Dictionary;
	(e.g. findfirst("BP1", InputDict["priority_channles"]["BP"])). if not,
	arbitrually pick up the first one.

	4. Then search the priority for first one;
	(e.g. findfirst("SP1", InputDict["priority_channles"]["BP"])).
	If not found, arbitrually pick up the second one as it is listed in
	priority dictionary.

	5. Compare the priority between first and second one, and pick up the
	earlier one as representive channel at this station.

III. Potential issue

	We assume that the physical measurement at same station does not change so much
	that the cross-correlation is not influenced by the replacement of channels.
	If it is not satisfied, it causes diference in the cross-correlation result.
	Please check the consistency between channels at same station if you find
	some discontinuous result before and after the switching of channel.

"""
function isolate_components(paths_all::AbstractArray, InputDict::Dict)

	iso_list = []

	for path in paths_all
		tmp = split(path, "/")[end]
		# read meta data from file name
		ftmpname = split(tmp, ".")
		println(ftmpname)
		if occursin("-", ftmpname[3])
			# format would be y, jd, T00-00-00, sta, loc, cha
			y, d, tmpT, net, sta, loc, cha = ftmpname
			#iso_stationinfo = (join([y, d, net, sta, loc], "-"), cha)
			# (time, net, sta for find the station, channel name and component)
			push!(iso_list, (path, join([y, d, tmpT, net, sta], "-"), net, sta, cha[1:2], cha[3]))

		else
			@warn "Format of tmp file is not y, jd, time, sta, loc, cha. Turn off IsIsolateComponents."
			return paths
		end
	end

	@show iso_list

	isocomp_idlist = []

	for (ista, current_sta) in enumerate(iso_list)
		for jsta = ista:length(iso_list)
			compared_sta = iso_list[jsta]
			if current_sta[2] != compared_sta[2] || current_sta[6] != compared_sta[6]
				# there is no conflict in channel, so add the current one to isocomp list
				if ista ∉ isocomp_idlist
					push!(isocomp_idlist, ista)
				end
			else
				# here current_sta[2] == compared_sta[2] && current_sta[6] == compared_sta[6]
				# perform process 2. we currently take into account priority with
				# all stations in certain network

				iso_net = compared_sta[3]
				current_cha = current_sta[5]
				compared_cha = compared_sta[5]

				if !haskey(InputDict["priority_channles"],iso_net)
					# second compared station does not have priority. take the current one
					if ista ∉ isocomp_idlist
						push!(isocomp_idlist, ista)
					end

				else
					# this has priority; perform process 3. e.g. compared_sta[5] = "SP"
					# @show compared_cha
					priority_compared = findfirst(x -> x == compared_cha, InputDict["priority_channles"][iso_net])
					if isempty(priority_compared)
						# this channel has no priority. take the current one
						if ista ∉ isocomp_idlist
							push!(isocomp_idlist, ista)
						end

					else
						# this has priority so that search priority for current one.
						priority_current = findfirst(x -> x == current_cha, InputDict["priority_channles"][iso_net])
						if isempty(priority_current)
							# current channel has no priority. take the compared one
							if jsta ∉ isocomp_idlist
								push!(isocomp_idlist, jsta)
							end

							# filter out ista if it's in iscomp_list
							filter!(x -> x != ista, isocomp_idlist)

						else
							# both current and compared one has priority. compare the order, and
							# add or replace ista in isocomp_idlist
							if priority_current > priority_compared
								# take current one
								if ista ∉ isocomp_idlist
									push!(isocomp_idlist, ista)
								end
							else
								# take compared one
								if jsta ∉ isocomp_idlist
									push!(isocomp_idlist, jsta)
								end
								# filter out ista if it's in iscomp_list
								filter!(x -> x != ista, isocomp_idlist)
							end
						end
					end
				end
			end
		end
	end

	#DEBUG:
	for id in isocomp_idlist
		temppath = paths_all[id]
		tmp = split(temppath, "/")[end]
		println(tmp)
	end

	return paths_all[isocomp_idlist]

end

"""
printparams(param::Dict)

print parameters
"""
function printparams(param::Dict)
    printstyled("-----------Input Parameters-----------\n"; color=:cyan, bold=true)
    for key in keys(param)

		try
	        if length(string(param["$key"])) > 60
	            param_str = string(param["$key"])[1:30]*"..."
	        else
	            param_str = string(param["$key"])
	        end
	        println(@sprintf("%-24s = %-10s", key, param_str))
		catch
		end
    end
end

"""
defaultinputdict(InputDict::Dict)

default if input parameter is missing.
"""
function defaultinputdict!(InputDict::Dict)

	# default values
	def = Dict()
	def["finame"]				 	= "./input.jld2"
	def["IsKurtosisRemoval"] 		= true
	def["max_edgetaper_duration"] 	= 60 * 5
	def["kurtosis_tw_sparse"] 		= 60
	def["kurtosis_timewindow"] 		= 60*3
	def["kurtosis_threshold"] 		= 2.0
	def["IsSTALTARemoval"] 			= true
	def["stalta_longtimewindow"] 	= 60*60*2
	def["stalta_threshold"] 		= 1.2
	def["stalta_absoluteclip"] 		= 1e20
	def["max_wintaper_duration"] 	= 60 * 3
	def["removal_shorttimewindow"] 	= 60 * 3
	def["overlap"] 					= 60
	def["IsIsolateComponents"] 		= false
	def["priority_channles"] 		= Dict()
	def["IsSaveFig"] 				= false
	def["plot_kurtosis_α"] 			= 1.2
	def["plot_boxheight	"] 			= 1.5
	def["plot_span"] 				= 100
	def["outputformat"]				= "JLD2"
	def["IsStartendtime"] 			= false
	def["fodir"] 					= "./dataset"
	def["foname"] 					= "eq_removed.jld2"

	#For dump everytrace
	#This is always false. If you need to plot some figures for
	#details of kurtosis and STA/LTA, please turn it true; it slowdowns down computation.
	def["dumptraces"] 				= false
	def["dumppath"]					= InputDict["fodir"]*"/dumptraces"

	for key in keys(def)
		if !haskey(InputDict, key)
			InputDict["$key"] = def["$key"]
		end
	end

end

"""
initlogo()

print initial logo
"""
function initlogo()

    print("
    _____      _      ______
   /  ___|    (_)     | ___ \\
   \\ `--.  ___ _ ___  | |_/ /___ _ __ ___   _____   _____
    `--. \\/ _ \\ / __| |    // _ \\ '_ ` _ \\ / _ \\ \\ / / _ \\
   /\\__/ /  __/ \\__ \\ | |\\ \\  __/ | | | | | (_) \\ V /  __/
   \\____/ \\___|_|___/ \\_| \\_\\___|_| |_| |_|\\___/ \\_/ \\___|
    ______ ____
   |  ____/ __ \\       |
   | |__ | |  | |      |
   |  __|| |  | |      |  v1.0 (Last update 07/07/2019)
   | |___| |__| |      |  © Kurama Okubo
   |______\\___\\_\\      |

")

    println("Job start running at "*string(now())*"\n")

end

end

max_stars_repo_name: kura-okubo/SeisRemoveEQ

max_issues_repo_name: kura-okubo/SeisRemoveEQ

max_forks_repo_name: kura-okubo/SeisRemoveEQ


text:
using StingerGraphs
using Base.Test

include("bfstest.jl")

s = Stinger(5, 0, [0 for i=1:6], [i%5 for i=1:5], [2 for i=1:5], Int64[], Int64[], -2)
@test insert_edge!(s, 0, 1, 5, 2, 2) == 1
@test outdegree(s, 1) == 1
@test indegree(s, 5) == 1
@test getsuccessors(s, 1) == [5]
@test insert_edge!(s, 0, 5, 2, 2, 2) == 1
@test insert_edge!(s, 0, 5, 4, 2, 2) == 1
@test insert_edge!(s, 0, 2, 0, 2, 2) == 1
@test insert_edge!(s, 0, 3, 4, 2, 2) == 1
@test insert_edge!(s, 0, 4, 0, 2, 2) == 1
@test bfs(s, 1, 6) == [2, -1, 5, -2, 5,  1]
@test remove_edge!(s, 0, 1, 5) == 1
@test outdegree(s, 1) == 0
@test indegree(s, 5) == 0
@test consistency_check(s, 100) == true
s[max_netypes] = 6
@test s[max_netypes] == 6

bfsutil(5, 10)


text:
using Metalhead
using Metalhead: classify

vgg = VGG19()
img = load("goma.jpeg")

@show classify(vgg, img)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module TypeFreeze

using InteractiveUtils  # For examining the outputs

# This is the idea we want to implement:
function widened_typeof_helper(::Type{T}) where {T}
    println("typed_helper(::$T)")
    widen(T)
end
function widened_typeof(::T) where {T}
    val = widened_typeof_helper(T)
    @eval widened_typeof(::$T) = $val
    return val
end

# Before it's called and memoized, the function is several lines long:
length((@code_typed widened_typeof(6))[1].code) > 1
# Calling the function memoizes its result as a newly defined method.
widened_typeof(Int64(3))
# Now, the entire body of the function is simply `return Int128`
(@code_typed widened_typeof(6))[1].code == [:(return $(Int128))]

include("utilities.jl")

"""
    @typefreeze function f(x) ... end

"Freezes" the output of a function for given parameter types.

Provides function output memoization, keyed by input types. This is implemented using
Julia's built-in dispatch mechanisms, by using `@eval` to generate a new method for each
unique input types tuple, whose function body is simply the return value.

This memoization occurs at the end of the first invocation, such that the function has a
one-time cost not unlike the compilation penalty incurred when type inference attempts to
compile a new instantiation of a method for new input parameters.
"""
macro typefreeze(funcexpr)
    typefreeze_helper(funcexpr, __module__)
end
function typefreeze_helper(funcexpr, __module__)
    dump(funcexpr)
    @assert(funcexpr.head == :function || Base.is_short_function_def(funcexpr))
    signature = funcexpr.args[1]
    func_callexpr = call_expr(signature)
    name = func_callexpr.args[1]

    escname = esc(name)
    args = func_callexpr.args[2:end]
    function_params = paramexprs(args)
    names = argnames(function_params)
    escnames = [esc(n) for n in argnames(function_params)]

    @show function_params

    # Special case for Type{T} types:
    #if isa(v, Expr) && v.head == :curly && v.args[1] == :Type
    #    v = v.args[2]
    #end
    #v

    esctypes = [esc(t) for t in argtypes(function_params)]
    escargs = [esc(a) for a in args]
    body = funcexpr.args[2]
    escbody = esc(body)

    main_signature = deepcopy(signature)
    call_expr(main_signature).args[2:end] = function_params

    # TODO: There is a bug in the julia macro-expander, preventing escaping the entire signature.
    # Therefore, for now, we have to manually escape its internals
    methods_helpername = gensym("$(name)_methods")
    methods_signature = deepcopy(main_signature)
    call_expr(methods_signature).args[1] = esc(methods_helpername)
    call_expr(methods_signature).args[2:end] = [esc(p) for p in function_params]
    # Escape the type parameters in where clauses (TODO: remove when bug is fixed)
    sigexpr = methods_signature
    while sigexpr.head != :call
        sigexpr.args[2:end] = [esc(e) for e in sigexpr.args[2:end]]
        sigexpr = sigexpr.args[1]
    end

    generator_helpername = gensym("$(name)_generator")
    generator_signature = deepcopy(signature)
    call_expr(generator_signature).args[1] = generator_helpername

    quote
        # This is the user's original function:
        $(esc(Expr(:function, generator_signature, body)))
        # This is the new implementation, which tacks new methods onto itself for each new
        # argument types tuple.
        # TODO: when the macro expander bug is fixed, esc the signature here instead.
        $(Expr(:function, (methods_signature), quote
            # Get actual return value from user's function:
            val = $(esc(generator_helpername))($(escnames...))
            # Compute the actual runtime types of the input arguments:
            types = [$(esctypes...)]
            typedargs = Tuple(Expr(Symbol("::"), t) for t in types)
            @eval $__module__ @inline $methods_helpername($(Expr(:$, :(typedargs...)))) = $(Expr(:$, :val))
            #@eval @show methods($methods_helpername)
            return val
        end))
        # The actual function simply delegates to the helper above
        # TODO: when the macro expander bug is fixed, esc the signature here instead.
        $(esc(Expr(:function, (main_signature), quote
            $methods_helpername($(names...))
        end)))
    end
end

#out1 = typefreeze_helper(:(myzero(x::Type{T}) where T = x), @__MODULE__)
@macroexpand @typefreeze myzero(x::Type{T}) where T = x

# -------- Illustration: ----------
# As an example, `@typefreeze function tzero(t::Tuple) Tuple(zero(x) for x in t) end` would
# produce the following code:
tzero_generator(t::Tuple) = Tuple(zero(x) for x in t)
function tzero_methods_helper(t::Tuple)
    val = tzero_generator(t)
    @eval @inline tzero_methods_helper(::$(typeof(t))) = $val
    return val
end
tzero(t::Tuple) = tzero_methods_helper(t)

tzero((1, 3, 2)) == (0, 0, 0)
tzero((1.0,)) == (0.0,)
# And the function is entirely memoized away:
@code_typed(tzero((1.0,)))[1].code == [:(return $((0.0,)))]

# ------- Tests: --------

@typefreeze function foo(x::Int8, y, ::Number)
    return x*y
end
methods(foo)
@code_typed foo(Int8(2), 5.0, 3)
foo(Int8(2), 5.0, 3)
methods(foo)
@code_typed foo(Int8(2), 5.0, 3)
foo(Int8(2), 4, 3)
foo(Int8(2), 6, 3)

@code_typed(foo(Int8(2), 5, 3))[1].code == [:(return $(8))]

# Real example:

# FixedPointDecimals.max_exp10(::Type{T}) where T
# In that package, this is manually "frozen" after its definition via manual @evals:
#   @eval max_exp10(::Type{Int128}) = $(max_exp10(Int128))  # Freeze for Int128, since it doesn't fold.

@typefreeze function max_exp10(::Type{T}) where T
    W = widen(T)
    type_max = W(typemax(T))

    powt = one(W)
    ten = W(10)
    exponent = 0

    while type_max > powt
        powt *= ten
        exponent += 1
    end

    exponent - 1
end

max_exp10(Int8)
@code_typed max_exp10(Int8)
# Even though this creates a BigInt, it can still be const-folded!
# This process is basically "manual (const) folding", where you're demanding folding even
# when the input might not be (provably) Const.
max_exp10(Int128)
@code_typed max_exp10(Int128)

end


text:
using Plasmons
using LinearAlgebra
using HDF5
using Test
using CUDA
using Adapt

CUDA.allowscalar(false)

@testset "Plasmons.jl" begin
    @testset "fermidirac" begin
        @test Plasmons.fermidirac(0.52193; mu = 0.4, kT = 0.1) ≈ 0.228059661030488549677
        @test Plasmons.fermidirac(2.0; mu = -0.8, kT = 0.13) ≈ 4.425527104542660880979e-10
        @test Plasmons.fermidirac(-2.0; mu = -0.8, kT = 0.13) ≈ 0.9999020317590553138168788
        @test Plasmons.fermidirac(100.0; mu = 0.82, kT = 0.1) ≈ 0
        @test Plasmons.fermidirac(-100.0; mu = 0.87, kT = 0.1) ≈ 1
    end

    @testset "g_matrix" begin
        for T in (Float32, Float64)
            E = rand(T, 10)
            μ = T(0.45)
            kT = T(0.01)
            ħω = complex(T)(0.5 + 1e-2im)
            kwargs = (mu = μ, kT = kT, blocks = false)
            G = Plasmons.g_matrix(complex(T), ħω, E; kwargs...)
            @test eltype(G) == complex(T)
            @test size(G) == (length(E), length(E))
            @test Diagonal(G) ≈ UniformScaling(0.0)
            if CUDA.functional()
                G₂ = Plasmons.g_matrix(complex(T), ħω, CuVector(E); kwargs...)
                @test eltype(G₂) == complex(T)
                @test size(G₂) == (length(E), length(E))
                @test typeof(G₂) <: CuArray
                @test G ≈ Array(G₂)
            end
        end
    end

    @testset "_analyze_top_left" begin
        G = rand(50, 50)
        fill!(view(G, 1:5, 1:7), zero(eltype(G)))
        G[6, 2] = 0.1
        @test Plasmons._analyze_top_left(G, 1e-2) == 5

        G = rand(50, 50)
        fill!(view(G, 1:5, 1:5), zero(eltype(G)))
        G[5, 6] = -0.1
        @test Plasmons._analyze_top_left(G, 1e-2) == 5

        G = rand(50, 50)
        fill!(view(G, :, 1:5), zero(eltype(G)))
        G[4, 6] = 0.1
        @test Plasmons._analyze_top_left(G, 1e-2) == 5

        G = rand(50, 50)
        fill!(view(G, :, 1:5), zero(eltype(G)))
        G[1, 1] = 0.1
        @test Plasmons._analyze_top_left(G, 1e-2) == 0
    end

    @testset "_analyze_bottom_right" begin
        G = rand(50, 50)
        fill!(view(G, 25:50, 27:50), zero(eltype(G)))
        G[26, 48] = 0.1
        @test Plasmons._analyze_bottom_right(G, 1e-2) == 24

        G = rand(50, 50)
        fill!(view(G, 27:50, 27:50), zero(eltype(G)))
        G[48, 27] = 0.1
        @test Plasmons._analyze_bottom_right(G, 1e-2) == 23
    end

    @testset "ThreeBlockMatrix" begin
        G = rand(50, 60) .- 0.5
        fill!(view(G, 1:3, 1:11), zero(eltype(G)))
        fill!(view(G, 35:50, 37:60), zero(eltype(G)))
        @test_throws DimensionMismatch Plasmons.ThreeBlockMatrix(G)

        for T in (Float32, Float64)
            H = rand(T, 50, 50)
            H = Hermitian(H + transpose(H))
            E = eigvals(H)
            μ = T(0.45)
            for kT in T[0.001, 0.01, 0.1, 1.0, 10.0]
                for _ω in -1.0:0.5:1.0
                    ħω = complex(T)(_ω + 1e-3im)
                    G₁ = Plasmons._build_g(ħω, E; mu = μ, kT = kT)
                    # full matrix
                    G₂ = Plasmons.ThreeBlockMatrix(G₁)
                    @test G₁ ≈ convert(Array, G₂)
                    # real and imaginary parts
                    G₂ = Plasmons.ThreeBlockMatrix(real(G₁))
                    @test real(G₁) ≈ convert(Array, G₂)
                    G₂ = Plasmons.ThreeBlockMatrix(imag(G₁))
                    @test imag(G₁) ≈ convert(Array, G₂)
                end
            end
        end
    end

    @testset "mul! for ThreeBlockMatrix" begin
        G = rand(50, 50) .- 0.5
        fill!(view(G, 1:3, 1:11), zero(eltype(G)))
        fill!(view(G, 35:50, 37:50), zero(eltype(G)))
        B = Plasmons.ThreeBlockMatrix(G)
        A = rand(38, 50)
        C = rand(38, 50)
        D = deepcopy(C)
        @test Plasmons._mul_with_strides!(C, A, B, 0.183, -0.482) ≈
              0.183 .* (A * G) .+ (-0.482) .* D

        for T in (Float32, Float64)
            for n in [1, 2, 3, 10, 50]
                for use_cuda in (false, true)
                    if use_cuda && !CUDA.functional()
                        continue
                    end
                    G = rand(T, n, n) .- 0.5
                    n₁ = rand(0:div(n, 2))
                    n₂ = rand(0:div(n, 2))
                    fill!(view(G, 1:n₁, 1:n₁), zero(eltype(G)))
                    fill!(view(G, (n - n₂ + 1):n, (n - n₂ + 1):n), zero(eltype(G)))
                    B = Plasmons.ThreeBlockMatrix(G)
                    use_cuda && (G = CuArray(G))
                    use_cuda && (B = adapt(CuArray, B))

                    A = rand(T, 7, n)
                    use_cuda && (A = CuArray(A))
                    C = rand(T, 7, n)
                    use_cuda && (C = CuArray(C))
                    D = deepcopy(C)
                    Plasmons._mul_with_strides!(C, A, B, T(0.183), T(-0.482))
                    Plasmons._mul_with_strides!(D, A, G, T(0.183), T(-0.482))
                    @test C ≈ D
                end
            end
        end
    end

    @testset "g_matrix with blocks" begin
        for T in (Float32, Float64)
            for use_cuda in (false, true)
                if use_cuda && !CUDA.functional()
                    continue
                end
                E = rand(T, 10)
                use_cuda && (E = CuArray(E))
                μ = T(0.45)
                kT = T(0.01)
                ħω = complex(T)(0.5 + 1e-2im)
                (Gᵣ, Gᵢ) = Plasmons.g_matrix(T, ħω, E; mu = μ, kT = kT, blocks = true)
                @test eltype(Gᵣ) == T
                @test size(Gᵣ) == (length(E), length(E))
                @test eltype(Gᵢ) == T
                @test size(Gᵢ) == (length(E), length(E))
            end
        end
    end

    @testset "dot_batched" begin
        reference(x, y) = sum(x .* y, dims = 2)
        for T in (Float32, Float64)
            for n in [1, 2, 3, 30, 50, 5000]
                for m in [1, 2, 3, 30, 50, 5000]
                    A = rand(n, m) .- 0.5
                    B = rand(n, m) .- 0.5
                    expected = reference(A, B)
                    @test Plasmons.dot_batched(A, B) ≈ expected
                    if CUDA.functional()
                        A = CuArray(A)
                        B = CuArray(B)
                        @test Array(Plasmons.dot_batched(A, B)) ≈ expected
                    end
                end
            end
        end
    end

    @testset "polarizability_batched" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            H = Hermitian(rand(T, 50, 50) .- T(0.5))
            E, ψ = eigen(H)
            μ = real(T)(0.45)
            for kT in real(T)[0.01, 0.1, 1.0, 10.0]
                for _ω in -1.0:0.5:1.0
                    ħω = complex(T)(_ω + 1e-3im)
                    args = (ħω, E, ψ)
                    kwargs = (method = :batched, mu = μ, kT = kT)
                    χ₁ = polarizability(args...; blocks = false, kwargs...)
                    χ₂ = polarizability(args...; blocks = true, kwargs...)
                    χ₃ = polarizability(args...; method = :simple, mu = μ, kT = kT)
                    @test χ₁ ≈ χ₃
                    @test χ₂ ≈ χ₃
                    if CUDA.functional()
                        args = (ħω, CuArray(E), CuArray(ψ))
                        χ₁ = Array(polarizability(args...; blocks = false, kwargs...))
                        χ₂ = Array(polarizability(args...; blocks = true, kwargs...))
                        @test χ₁ ≈ χ₃
                        @test χ₂ ≈ χ₃
                    end
                end
            end
        end
    end

    @testset "dimension checks" begin
        ħω = 0.23 + 0.02im
        μ = 0.2
        kT = 0.5
        kwargs = (method = :simple, mu = μ, kT = kT)
        @test_throws DimensionMismatch polarizability(ħω, rand(10), rand(10, 9); kwargs...)
        @test_throws DimensionMismatch polarizability(ħω, rand(10), rand(9, 9); kwargs...)
        @test_throws DimensionMismatch dielectric(rand(ComplexF64, 10, 9), rand(10, 10))
        @test_throws DimensionMismatch dielectric(rand(ComplexF64, 9, 9), rand(10, 10))
        @test_throws DimensionMismatch dielectric(rand(ComplexF64, 10, 9), rand(10, 9))
    end
end