text:
# Level-3 LinearAlgebra.BLAS interface.
#

# Level3 BLAS size checker.
bli_check_lv3(transa::BliTrans,
              transb::BliTrans,
              transc::BliTrans,
              m ::Integer, k ::Integer,
              k_::Integer, n ::Integer,
              m_::Integer, n_::Integer) = begin
    # Assumption:
    # m, k * k_, n -> m_, n_

    ((transa.enum & BLIS_TRANS_BIT) != 0) && ((m , k ) = (k , m ))
    ((transb.enum & BLIS_TRANS_BIT) != 0) && ((k_, n ) = (n , k_))
    ((transc.enum & BLIS_TRANS_BIT) != 0) && ((m_, n_) = (n_, m_))

    (m == m_ &&
     n == n_) || throw(DimensionMismatch("Target buffer size mismatch."))
    (k == k_) || throw(DimensionMismatch("Contracted size mismatch."))

    nothing
end

# Julia has a specific `gemm_wrapper!` for spetial GEMM calls.
# In case of BLIS these separations are not needed so this function
#   calls directly `gemm!`.
macro blis_interface_linalg_lv3_gemm_wrapper(T1, T2, T3)

    wrapperfunc = :( gemm_wrapper! )

    return quote

        $(esc(wrapperfunc))(C::StridedVecOrMat{$T1},
                            tA::AbstractChar,
                            tB::AbstractChar,
                            A::StridedVecOrMat{$T2},
                            B::StridedVecOrMat{$T3},
                            _coeffs = LinearAlgebra.MulAddMul()) = begin

            α = 1.0 * _coeffs.alpha
            β = 1.0 * _coeffs.beta

            # TODO: 2 small cases.
            gemm!(tA, tB, α, A, B, β, C)

        end
    end
end

macro blis_interface_linalg_lv3_gemm(Tc1, T1, T2, Tc2, T3, targetfunc, bliname)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(tA::AbstractChar,
                           tB::AbstractChar,
                           α::$Tc1,
                           A::StridedVecOrMat{<:$T1},
                           B::StridedVecOrMat{<:$T2},
                           β::$Tc2,
                           C::StridedVecOrMat{<:$T3}) = begin

            bli_tA = char_to_trans[tA]
            bli_tB = char_to_trans[tB]

            bli_check_lv3(bli_tA, bli_tB, BLIS_NO_TRANSPOSE,
                          size(A)...,
                          size(B)...,
                          size(C)...)

            oα = BliObj(α)
            oA = BliObj(A)
            oB = BliObj(B)
            oβ = BliObj(β)
            oC = BliObj(C)

            ObjectBackend.bli_obj_set_onlytrans!(bli_tA, oA.obj)
            ObjectBackend.bli_obj_set_onlytrans!(bli_tB, oB.obj)
            $blifunc(oα, oA, oB, oβ, oC)
            C

        end
    end
end

if "gemm" ∉ blacklist
@doc """
        gemm!(tA, tB, α, A, B, β, C)
    BLIS-based GEMM with generic strides & mixed precision directly supported:
    ```math
    C = β C + α A^{tA} B^{tB}
    ```
    """ gemm!
@blis_interface_linalg_lv3_gemm(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                gemm!, gemm!)
# To avoid obfuscating with LinearAlgebra.gemm!,
#  instantiate more specifically.
@blis_interface_linalg_lv3_gemm Float32    Float32    Float32    Float32    Float32    gemm! gemm!
@blis_interface_linalg_lv3_gemm Float64    Float64    Float64    Float64    Float64    gemm! gemm!
@blis_interface_linalg_lv3_gemm ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 gemm! gemm!
@blis_interface_linalg_lv3_gemm ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 gemm! gemm!

# Do this for gemm_wrapper! as well.
# Direct overriding triggers warning.
@blis_interface_linalg_lv3_gemm_wrapper Float32    Float32    Float32
@blis_interface_linalg_lv3_gemm_wrapper Float64    Float64    Float64
@blis_interface_linalg_lv3_gemm_wrapper ComplexF32 ComplexF32 ComplexF32
@blis_interface_linalg_lv3_gemm_wrapper ComplexF64 ComplexF64 ComplexF64
end

macro blis_interface_linalg_lv3_hemm(Tc1, T1, T2, Tc2, T3, targetfunc, bliname, bli_struc)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(si::AbstractChar,
                           ul::AbstractChar,
                           α::$Tc1,
                           A::StridedMatrix{<:$T1},
                           B::StridedMatrix{<:$T2},
                           β::$Tc2,
                           C::StridedMatrix{<:$T3}) = begin

            bli_si = char_to_side[si]
            bli_ul = char_to_uplo[ul]

            if bli_si == BLIS_LEFT
                bli_check_lv3(BLIS_NO_TRANSPOSE, BLIS_NO_TRANSPOSE, BLIS_NO_TRANSPOSE,
                              size(A)...,
                              size(B)...,
                              size(C)...)
            else
                bli_check_lv3(BLIS_NO_TRANSPOSE, BLIS_NO_TRANSPOSE, BLIS_NO_TRANSPOSE,
                              size(B)...,
                              size(A)...,
                              size(C)...)
            end

            oα = BliObj(α)
            oA = BliObj(A)
            oB = BliObj(B)
            oβ = BliObj(β)
            oC = BliObj(C)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oA.obj)
            ObjectBackend.bli_obj_set_struc!($bli_struc, oA.obj)
            $blifunc(bli_si, oα, oA, oB, oβ, oC)
            C

        end
    end
end

if "hemm" ∉ blacklist
@doc """
        hemm!(side, uplo, α, A, B, β, C)
    BLIS-based HEMM with generic strides & mixed precision directly supported.
    `A` expresses a `Hermitian` matrix with its `uplo` triangle.
    """ hemm!
@blis_interface_linalg_lv3_hemm(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                hemm!, hemm!,
                                BLIS_HERMITIAN)
@blis_interface_linalg_lv3_hemm Float32    Float32    Float32    Float32    Float32    hemm! hemm! BLIS_HERMITIAN
@blis_interface_linalg_lv3_hemm Float64    Float64    Float64    Float64    Float64    hemm! hemm! BLIS_HERMITIAN
@blis_interface_linalg_lv3_hemm ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 hemm! hemm! BLIS_HERMITIAN
@blis_interface_linalg_lv3_hemm ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 hemm! hemm! BLIS_HERMITIAN
end

if "symm" ∉ blacklist
@doc """
        symm!(side, uplo, α, A, B, β, C)
    BLIS-based SYMM with generic strides & mixed precision directly supported.
    `A` expresses a `Symmetric` matrix with its `uplo` triangle.
    """ symm!
@blis_interface_linalg_lv3_hemm(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                symm!, symm!,
                                BLIS_SYMMETRIC)
@blis_interface_linalg_lv3_hemm Float32    Float32    Float32    Float32    Float32    symm! symm! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_hemm Float64    Float64    Float64    Float64    Float64    symm! symm! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_hemm ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 symm! symm! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_hemm ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 symm! symm! BLIS_SYMMETRIC
end

macro blis_interface_linalg_lv3_her2k(Tc1, T1, T2, Tc2, T3, targetfunc, bliname, bli_struc)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(ul::AbstractChar,
                           tAB::AbstractChar,
                           α::$Tc1,
                           A::StridedMatrix{<:$T1},
                           B::StridedMatrix{<:$T2},
                           β::$Tc2,
                           C::StridedMatrix{<:$T3}) = begin

            bli_ul = char_to_uplo[ul]
            bli_tAB = char_to_trans[tAB]

            bli_check_lv3(bli_tAB, bli_tAB, BLIS_NO_TRANSPOSE,
                          size(A)...,
                          size(B')...,
                          size(C)...)

            oα = BliObj(α)
            oA = BliObj(A)
            oB = BliObj(B)
            oβ = BliObj(β)
            oC = BliObj(C)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oC.obj)
            ObjectBackend.bli_obj_set_struc!($bli_struc, oC.obj)
            ObjectBackend.bli_obj_set_onlytrans!(bli_tAB, oA.obj)
            ObjectBackend.bli_obj_set_onlytrans!(bli_tAB, oB.obj)
            $blifunc(oα, oA, oB, oβ, oC)
            C

        end
    end
end

if "her2k" ∉ blacklist
@doc """
        her2k!(uplo, tAB, α, A, B, β, C)
    BLIS-based HER2K with generic strides & mixed precision directly supported.
    Performs rank-2k update on `Hermitian` matrix `C` (expressed by `uplo`-triangle):
    ```math
    C = β C + (α A B^† + \\bar α B A^†)^{tAB}
    ```
    """ her2k!
@blis_interface_linalg_lv3_her2k(BliCompatibleType,
                                 BliCompatibleType,
                                 BliCompatibleType,
                                 BliCompatibleType,
                                 BliCompatibleType,
                                 her2k!, her2k!,
                                 BLIS_HERMITIAN)
@blis_interface_linalg_lv3_her2k Float32    Float32    Float32    Float32 Float32    her2k! her2k! BLIS_HERMITIAN
@blis_interface_linalg_lv3_her2k Float64    Float64    Float64    Float64 Float64    her2k! her2k! BLIS_HERMITIAN
@blis_interface_linalg_lv3_her2k ComplexF32 ComplexF32 ComplexF32 Float32 ComplexF32 her2k! her2k! BLIS_HERMITIAN
@blis_interface_linalg_lv3_her2k ComplexF64 ComplexF64 ComplexF64 Float64 ComplexF64 her2k! her2k! BLIS_HERMITIAN
end

if "syr2k" ∉ blacklist
@doc """
        syr2k!(uplo, tAB, α, A, B, β, C)
    BLIS-based SYR2K with generic strides & mixed precision directly supported.
    Performs rank-2k update on `Symmetric` matrix `C` (expressed by `uplo`-triangle):
    ```math
    C = β C + (α A B^T + \\bar α B A^T)^{tAB}
    ```
    """ syr2k!
@blis_interface_linalg_lv3_her2k(BliCompatibleType,
                                 BliCompatibleType,
                                 BliCompatibleType,
                                 BliCompatibleType,
                                 BliCompatibleType,
                                 syr2k!, syr2k!,
                                 BLIS_SYMMETRIC)
@blis_interface_linalg_lv3_her2k Float32    Float32    Float32    Float32    Float32    syr2k! syr2k! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_her2k Float64    Float64    Float64    Float64    Float64    syr2k! syr2k! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_her2k ComplexF32 ComplexF32 ComplexF32 ComplexF32 ComplexF32 syr2k! syr2k! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_her2k ComplexF64 ComplexF64 ComplexF64 ComplexF64 ComplexF64 syr2k! syr2k! BLIS_SYMMETRIC
end

macro blis_interface_linalg_lv3_herk(Tc1, T1, Tc2, T2, targetfunc, bliname, bli_struc)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(ul::AbstractChar,
                           tA::AbstractChar,
                           α::$Tc1,
                           A::StridedMatrix{<:$T1},
                           β::$Tc2,
                           C::StridedMatrix{<:$T2}) = begin

            bli_ul = char_to_uplo[ul]
            bli_tA = char_to_trans[tA]

            bli_check_lv3(bli_tA, bli_tA, BLIS_NO_TRANSPOSE,
                          size(A)...,
                          size(A')...,
                          size(C)...)

            oα = BliObj(α)
            oA = BliObj(A)
            oβ = BliObj(β)
            oC = BliObj(C)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oC.obj)
            ObjectBackend.bli_obj_set_struc!($bli_struc, oC.obj)
            ObjectBackend.bli_obj_set_onlytrans!(bli_tA, oA.obj)
            $blifunc(oα, oA, oβ, oC)
            C

        end
    end
end

if "herk" ∉ blacklist
@doc """
        herk!(uplo, tA, α, A, β, C)
    BLIS-based HERK with generic strides & mixed precision directly supported.
    Performs rank-k update on `Hermitian` matrix `C` (expressed by `uplo`-triangle).
    """ herk!
@blis_interface_linalg_lv3_herk(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                herk!, herk!,
                                BLIS_HERMITIAN)
@blis_interface_linalg_lv3_herk Float32 Float32    Float32 Float32    herk! herk! BLIS_HERMITIAN
@blis_interface_linalg_lv3_herk Float64 Float64    Float64 Float64    herk! herk! BLIS_HERMITIAN
@blis_interface_linalg_lv3_herk Float32 ComplexF32 Float32 ComplexF32 herk! herk! BLIS_HERMITIAN
@blis_interface_linalg_lv3_herk Float64 ComplexF64 Float64 ComplexF64 herk! herk! BLIS_HERMITIAN
end

if "syrk" ∉ blacklist
@doc """
        syrk!(uplo, tA, α, A, β, C)
    BLIS-based SYRK with generic strides & mixed precision directly supported.
    Performs rank-k update on `Symmetric` matrix `C` (expressed by `uplo`-triangle).
    """ syrk!
@blis_interface_linalg_lv3_herk(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                syrk!, syrk!,
                                BLIS_SYMMETRIC)
@blis_interface_linalg_lv3_herk Float32    Float32    Float32    Float32    syrk! syrk! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_herk Float64    Float64    Float64    Float64    syrk! syrk! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_herk ComplexF32 ComplexF32 ComplexF32 ComplexF32 syrk! syrk! BLIS_SYMMETRIC
@blis_interface_linalg_lv3_herk ComplexF64 ComplexF64 ComplexF64 ComplexF64 syrk! syrk! BLIS_SYMMETRIC
end

macro blis_interface_linalg_lv3_trmm(Tc1, T1, T2, targetfunc, bliname)

    # Get method for the object API backend.
    blifuncname = Symbol("bli_", bliname)
    blifunc = getproperty(ObjectBackend, blifuncname)

    return quote

        $(esc(targetfunc))(si::AbstractChar,
                           ul::AbstractChar,
                           tA::AbstractChar,
                           dA::AbstractChar,
                           α::$Tc1,
                           A::StridedMatrix{<:$T1},
                           B::StridedMatrix{<:$T2}) = begin

            bli_si = char_to_side[si]
            bli_ul = char_to_uplo[ul]
            bli_dA = char_to_diag[dA]
            bli_tA = char_to_trans[tA]

            if bli_si == BLIS_LEFT
                bli_check_lv3(bli_tA, BLIS_NO_TRANSPOSE, BLIS_NO_TRANSPOSE,
                              size(A)...,
                              size(B)...,
                              size(B)...)
            else
                bli_check_lv3(BLIS_NO_TRANSPOSE, bli_tA, BLIS_NO_TRANSPOSE,
                              size(B)...,
                              size(A)...,
                              size(B)...)
            end

            oα = BliObj(α)
            oA = BliObj(A)
            oB = BliObj(B)

            ObjectBackend.bli_obj_set_uplo!(bli_ul, oA.obj)
            ObjectBackend.bli_obj_set_diag!(bli_dA, oA.obj)
            ObjectBackend.bli_obj_set_onlytrans!(bli_tA, oA.obj)
            ObjectBackend.bli_obj_set_struc!(BLIS_TRIANGULAR, oA.obj)
            $blifunc(bli_si, oα, oA, oB)
            B

        end
    end
end

if "trmm" ∉ blacklist
@doc """
        trmm!(side, uplo, tA, dA, α, A, B)
    BLIS-based TRMM with generic strides & mixed precision directly supported.
    """ trmm!
@blis_interface_linalg_lv3_trmm(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                trmm!, trmm!)
@blis_interface_linalg_lv3_trmm Float32    Float32    Float32    trmm! trmm!
@blis_interface_linalg_lv3_trmm Float64    Float64    Float64    trmm! trmm!
@blis_interface_linalg_lv3_trmm ComplexF32 ComplexF32 ComplexF32 trmm! trmm!
@blis_interface_linalg_lv3_trmm ComplexF64 ComplexF64 ComplexF64 trmm! trmm!
end

if "trsm" ∉ blacklist
@doc """
        trsm!(side, uplo, tA, dA, α, A, B)
    BLIS-based TRSM with generic strides & mixed precision directly supported.
    """ trsm!
@blis_interface_linalg_lv3_trmm(BliCompatibleType,
                                BliCompatibleType,
                                BliCompatibleType,
                                trsm!, trsm!)
@blis_interface_linalg_lv3_trmm Float32    Float32    Float32    trsm! trsm!
@blis_interface_linalg_lv3_trmm Float64    Float64    Float64    trsm! trsm!
@blis_interface_linalg_lv3_trmm ComplexF32 ComplexF32 ComplexF32 trsm! trsm!
@blis_interface_linalg_lv3_trmm ComplexF64 ComplexF64 ComplexF64 trsm! trsm!
end


text:
abstract type EasingType end

struct Easing{T <: EasingType}
    easing::T
    n::Int
    yoyo::Bool
    prewait::Float64
    postwait::Float64
end

struct FuncEasing <: EasingType
    f::Function
    args::Tuple

    function FuncEasing(f, args)
        f0 = f(0, args...)
        f1 = f(1, args...)
        if !(f0 == 0 && f1 == 1)
            error("The given easing function is not 0 and 1 at t = 0 and t = 1, but $f0 and $f1")
        end
        new(f, args)
    end
end

FuncEasing(f::Function) = FuncEasing(f, ())

"""
Mixes two easings with a constant ratio r, so that final = e1 * r + e2 * (1 - r)
"""
struct MixedEasing{S <: Easing, T <: Easing} <: EasingType
    e1::S
    e2::T
    mix::Float64

    function MixedEasing(e1::S, e2::T, mix::Float64=0.5) where {S <: Easing, T <: Easing}
        if !(0 <= mix <= 1)
            error("Mix has to be between 0 and 1, but is $mix")
        end
        new{S, T}(e1, e2, mix)
    end
end

"""
Multiplies two easings with each other, so final = e1 * e2
"""
struct MultipliedEasing{S <: Easing, T <: Easing} <: EasingType
    e1::S
    e2::T
end

"""
Mixes two easings with a factor decided with a third easing, so that
final = (1 - easing()) * e1 + easing() * e2
"""
struct EasedEasing{S <: Easing, T <: Easing, U <: Easing} <: EasingType
    e1::S
    e2::T
    easing::U

    function EasedEasing(e1, e2, easing=Easing(easing=LinearEasing()))
        new{typeof(e1), typeof(e2), typeof(easing)}(e1, e2, easing)
    end
end

Easing(easing = FuncEasing(f_linear); n=1, yoyo=false, prewait=0.0, postwait=0.0) = Easing(easing, n, yoyo, prewait, postwait)

function opposite(f)
    function opfunc(fraction, args...)
        1 - f(1 - fraction, args...)
    end
end

funcease(f, args...; kwargs...) = Easing(FuncEasing(f, args); kwargs...)

f_noease(fraction) = fraction == 1 ? 1 : 0
noease(;kwargs...) = funcease(f_noease; kwargs...)

f_stepease(fraction, stepfraction) = fraction < stepfraction ? 0 : 1
stepease(stepfraction; kwargs...) = funcease(f_stepease, stepfraction; kwargs...)

f_sineio(fraction) = sin(pi * fraction - 0.5pi) * 0.5 + 0.5
sineio(;kwargs...) = funcease(f_sineio; kwargs...)

f_saccadic(fraction, power) = -(sin((-fraction + 1) ^ power * pi - pi/2) * 0.5 + 0.5) + 1
saccadic(power; kwargs...) = funcease(f_saccadic, power; kwargs...)

f_linear(fraction) = fraction
linear(;kwargs...) = funcease(f_linear; kwargs...)

f_polyin(fraction, power) = fraction ^ power
polyin(power; kwargs...) = funcease(f_polyin, power; kwargs...)

polyout(power; kwargs...) = funcease(opposite(f_polyin), power; kwargs...)

f_polyio(fraction, power) = (
    fraction < 0.5 ?
    f_polyin(2 * fraction, power) / 2 :
    opposite(f_polyin)(2 * (fraction - 0.5), power) / 2 + 0.5)
polyio(power; kwargs...) = funcease(f_polyio, power; kwargs...)

f_expin(fraction, exponent) = ((exponent ^ fraction) - 1) / (exponent - 1)
expin(exponent; kwargs...) = funcease(f_expin, exponent; kwargs...)

expout(exponent; kwargs...) = funcease(opposite(f_expin), exponent; kwargs...)

mixed(e1, e2; kwargs...) = Easing(MixedEasing(e1, e2); kwargs...)
multiplied(e1, e2; kwargs...) = Easing(MultipliedEasing(e1, e2); kwargs...)
eased(e1, e2, with_e; kwargs...) = Easing(EasedEasing(e1, e2, with_e); kwargs...)

Base.Broadcast.broadcastable(e::Easing) = Ref(e)

ext: jl


text:
using PyCall,Distributions,DataStructures
@pyimport networkx as nx
@pyimport matplotlib.pyplot as plt
function generator_graphs(L,m=2)
    g=nx.Graph()
    h=nx.barabasi_albert_graph(L,m)
    d=Normal(1,0.33)
    edges=h[:edges]()
    for edge in edges
        g[:add_edge](edge[1]+1,edge[2]+1,weight=abs(rand(d)))
    end
    return g
end

function _all_pairs_dijkstra_path(g)
    path=nx.all_pairs_dijkstra_path(g)
    return path
end

function _all_pairs_shortest_path_length(g)
    pathlength=nx.all_pairs_shortest_path_length(g)
    return pathlength
end

type Package
    size_p::Float64
    src_p::Int
    dis_p::Int
    pos_p::Int
    time_p::Int
    step_p::Int
    priority::Int
end

function package_brith(g,package,packagequeue,P,N,t)
    for i =1:P
        d=Normal(1,0.33)
        #n=nx.number_of_nPriorityes(g)
        sizep=abs(rand(d))
        srcp=sample(1:N)
        disp=sample(1:N)
        priorityp=sample(1:5)
        package[(P*(t-1)+i)]=Package(sizep,srcp,disp,srcp,0,0,priorityp)
        packagequeue[package[(P*(t-1)+i)].pos_p][P*(t-1)+i]=package[P*(t-1)+i].priority
    end
    return package,packagequeue
end

function _load(load_arr,package)
    fill!(load_arr,0)
    load=0
    for (k,v) in package
        k_pos=package[k].pos_p
        k_size=package[k].size_p
        load_arr[(k_pos)]+=k_size
        package[k].time_p+=1 #update package time
        load+=k_size
    end
    return load_arr,load
end

function effective_distance(g,k,j,path,load_arr,h)
    ekj=h*length(path[k][j])+(1-h)*load_arr[k]
    return ekj
end

function traffic_aware_routing_strategy(i,g,package,package_life,package_step,packagequeue,path,load_arr,h)
    pos=package[i].pos_p
    dis=package[i].dis_p
    if pos == dis
        package_life[i]=package[i].time_p
        package_step[i]=package[i].step_p
        delete!(package,i)
    else
        nebor_arr=nx.neighbors(g,pos)
        min_effecitve_distance=10000
        best_k=0
        for j=1:length(nebor_arr)
            ekj=effective_distance(g,nebor_arr[j],dis,path,load_arr,h)
            if min_effecitve_distance>ekj
                min_effecitve_distance=ekj
                best_k=nebor_arr[j]
            end
        end
        package[i].pos_p=best_k
        package[i].step_p+=1
        if best_k==dis
            package[i].time_p+=1
            package_life[i]=package[i].time_p
            package_step[i]=package[i].step_p
            delete!(package,i)
        else
            packagequeue[best_k][i]=package[i].priority
        end
    end
end

function initiate_queue(packagequeue,N)
    for i=1:N
        packagequeue[i]=Collections.PriorityQueue()
    end
end

function Package_move(packagequeue,move_package,N)
    for i=1:N
        if !isempty(packagequeue[i])
            push!(move_package,Collections.dequeue!(packagequeue[i]))
        end
    end
end

function update_package(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,h)
    while !isempty(move_package)
        k=pop!(move_package)
        traffic_aware_routing_strategy(k,g,package,package_life,package_step,packagequeue,path,load_arr,h)
    end
end

function initiate_graph(L)
    # initiate a graph
    g=generator_graphs(L)
    #g=remove_edges(g)
    #number of nodes
    N=nx.number_of_nodes(g)
    D=nx.diameter(g)
    path=_all_pairs_dijkstra_path(g)
    pathlength=_all_pairs_shortest_path_length(g)
    return g,N,D,path,pathlength
end

function initiate_package(L,N)
    # initiate a package
    package=Dict{Int,Package}()
    package_life=Dict{Int,Int}()
    package_step=Dict{Int,Int}()
    # initiate the moving package
    move_package=Set{Int}()
    # initiage the queue of nodes
    packagequeue=Array(Collections.PriorityQueue,N)
    initiate_queue(packagequeue,N)
    # initiate a load array
    load_arr=zeros(N)
    total_load=Float64[]
    return package,package_life,package_step,move_package,packagequeue,load_arr,total_load
end

function reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
    for t=1:T
        # generate a package
        package,packagequeue=package_brith(g,package,packagequeue,P,N,t)
        Package_move(packagequeue,move_package,N)
        # update the load array
        load_arr,load=_load(load_arr,package)
        # update the package
        update_package(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,h)
        #calculate the order parameter
        push!(total_load,load)
    end
end

function free_state_nodes(s,load_arr,m,N)
    empty!(s)
    S=0
    for i=1:N
        if load_arr[i]<=m
            push!(s,i)
            S+=1
        end
    end
    return s,S
end

function largest_cluster(g,s,load_arr,m,N)
    s,S=free_state_nodes(s,load_arr,m,N)
    G=nx.subgraph(g,s)
    SG=nx.connected_component_subgraphs(G)
    largest=length(SG[1])
    second_largest= (length(SG)>1)? length(SG[2]) : 0
    return S,largest,second_largest
end

function giant_component(g,s,load_arr,mc,N)
    s,S=free_state_nodes(s,load_arr,mc,N)
    G=nx.subgraph(g,s)
    SG=nx.connected_component_subgraphs(G)
    return SG[1]
end
function find_pc(p,sl,s,g,N,load_arr,m_min,m_max)
    empty!(p)
    empty!(sl)
    Δm1=0.2
    Δm2=0.02
    m=m_min
    while m<=m_max
        S,largest,second_largest=largest_cluster(g,s,load_arr,m,N)
        push!(p,S/N)
        push!(sl,second_largest/N)
        m += (0.3<S/N<0.7) ?  Δm2 : Δm1
    end
    max_index=findmax(sl)[2]
    pc=p[max_index]
    return pc
end

function find_mc(mm,sl,s,g,N,load_arr,m_min,m_max)
    empty!(mm)
    empty!(sl)
    Δm1=0.2
    Δm2=0.02
    m=m_min
    while m<=m_max
        S,largest,second_largest=largest_cluster(g,s,load_arr,m,N)
        push!(mm,m)
        push!(sl,second_largest/N)
        m += (0.3<S/N<0.7) ?  Δm2 : Δm1
    end
    max_index=findmax(sl)[2]
    mc=mm[max_index]
    return mc
end

function _critical_p(p,sl,s,g,N,load_arr)
    # find the minimum m and minmax m
    m_min=minimum(load_arr)
    m_max=maximum(load_arr)
    # critical threshold
    pc=find_pc(p,sl,s,g,N,load_arr,m_min,m_max)
    return pc
end

function _ρ(total_load,P,Δt,T)
    ρ=0
    T0=T-0.1*T+1
    for t=T0:T
        ρ+=(total_load[t]-total_load[t-Δt])/(Δt*P)
    end
    return ρ/(T-T0)
end

function rs(pathlength,N,D)
    R=Array(Vector{Tuple{Int,Int}},D)
    for i =1:D
        R[i]=Vector{Tuple{Int,Int}}[]
    end
    for i=1:N,j=1:N
        if pathlength[i][j] != 0
        append!(R[pathlength[i][j]],[(i,j)])
        end
    end
    return R
end

function _cor(σ2,cor_data,R,r)
    C=0
    if !isempty(R[r])
        s=0
        for i=1:length(R[r])
            xi=R[r][i][1]
            xj=R[r][i][2]
            s+=cor_data[xi]*cor_data[xj]
        end
        C=s/(length(R[r])*σ2)
    else
        C=0
    end
    return C
end

function _ζ(load_arr,pathlength,N,D)
    average_load=mean(load_arr)
    σ2=var(load_arr)
    cor_data=load_arr-average_load
    r=1
    ζ=0
    R=rs(pathlength,N,D)
    while r<=D
        C=_cor(σ2,cor_data,R,r)
        if C<=0
            ζ=r
            break
        end
        r+=1
    end
    return ζ
end

function average_bw(mm,sl,s,g,N,load_arr)
    m_min=minimum(load_arr)
    m_max=maximum(load_arr)
    mc=find_mc(mm,sl,s,g,N,load_arr,m_min,m_max)
    IIC=giant_component(g,s,load_arr,mc,N)
    bw=nx.betweenness_centrality(IIC)  #Dict
    average_C=mean(collect(values(bw)))
    return average_C
end
function reach_stable_time(P=150,h=0.3,L=100,T=10000)
    outfreachstabletime=open("/home/SuZ/Data/Priority/ScaleFree_ReachStableTime_($h)_($L)_($P).txt","w")
    g,N,D,path,pathlength=initiate_graph(L)
    package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
    # reach stable state
    for t=1:T
        # generate a package
        package,packagequeue=package_brith(g,package,packagequeue,P,N,t)
        Package_move(packagequeue,move_package,N)
        # update the load array
        load_arr,load=_load(load_arr,package)
        # update the package
        update_package(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,h)
        #calculate the order parameter
        push!(total_load,load)
    end
    #ρ+=_ρ(total_load,P,Δt,T)
    writedlm(outfreachstabletime,total_load)
    close(outfreachstabletime)
end
function order_parameter(Pstat=10,Pend=100,ΔP=2,samples=5,h=0.3,L=100,T=10000,Δt=10)
    outforderparameter=open("/home/SuZ/Data/Priority/ScaleFree_Orderparameter_($h)_($L).txt","w")
    g,N,D,path,pathlength=initiate_graph(L)
    P=Pstat
    while P<=Pend
        ρ=0
        for i=1:samples
            #initiate
            package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
            # reach stable state
            reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
            ρ+=_ρ(total_load,P,Δt,T)
        end
        writedlm(outforderparameter,transpose([P,ρ/samples]))
        P+=ΔP
    end
    close(outforderparameter)
end

function percolation(P=150,h=0.3,L=100,T=10000)
    #outfile
    outfpercolation=open("/home/SuZ/Data/Priority/ScaleFree_Percolation_($P)_($h)_($L).txt","w")
    #parameters array
    s=Int[]
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    sizehint!(s,N)
    package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
    # reach stable state
    reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
    # find the minimum m and minmax m
    m_min=minimum(load_arr)
    m_max=maximum(load_arr)
    m=m_min
    Δm=0.2
    while m<=m_max
        S,largest,second_largest=largest_cluster(g,s,load_arr,m,N)
        writedlm(outfpercolation,transpose([S/N,largest/N,second_largest/N]))
        m+=Δm
    end
    close(outfpercolation)
end

function critical_threshold(Pstat=10,Pend=100,ΔP=2,samples=5,h=0.3,L=100,T=10000)
    outfcriticalthreshold=open("/home/SuZ/Data/Priority/ScaleFree_CriticalThreshold_($h)_($L).txt","w")
    #parameters array
    p=Float64[]
    sl=Float64[]
    s=Int[]
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    sizehint!(s,N)
    P=Pstat
    while P<=Pend
        pc=0
        for i=1:samples
            #initiate
            package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
            # reach stable state
            reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
            pc+=_critical_p(p,sl,s,g,N,load_arr)
        end
        writedlm(outfcriticalthreshold,transpose([P,pc/samples]))
        P+=ΔP
    end
    close(outfcriticalthreshold)
end

function correlation(P=150,h=0.3,L=100,T=10000)
    #outfile
    outfcorrelaiton=open("/home/SuZ/Data/Priority/ScaleFree_Correlation_($P)_($h)_($L).txt","w")
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
    # reach stable state
    reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
    average_load=mean(load_arr)
    σ2=var(load_arr)
    cor_data=load_arr-average_load
    R=rs(pathlength,N,D)
    r=1
    while r<=D
        C=_cor(σ2,cor_data,R,r)
        writedlm(outfcorrelaiton,transpose([r,C]))
        r+=1
    end
    close(outfcorrelaiton)
end

function correlation_length(Pstat=10,Pend=100,ΔP=2,samples=5,h=0.3,L=100,T=10000)
    outfcorrelationlength=open("/home/SuZ/Data/Priority/ScaleFree_Correlationlength_($h)_($L).txt","w")
    g,N,D,path,pathlength=initiate_graph(L)
    P=Pstat
    while P<=Pend
        ζ=0
        for i=1:samples
            #initiate
            package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
            # reach stable state
            reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
            ζ+=_ζ(load_arr,pathlength,N,D)
        end
        writedlm(outfcorrelationlength,transpose([P,ζ/samples]))
        P+=ΔP
    end
    close(outfcorrelationlength)
end

function betweenness_distribution(P=150,h=0.3,L=100,T=10000)
    outfbetweenness=open("/home/SuZ/Data/Priority/ScaleFree_Betweenness_($P)_($h)_($L).txt","w")
    #parameters array
    mm=Float64[]
    sl=Float64[]
    s=Int[]
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    sizehint!(s,N)
    package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
    # reach stable state
    reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
    m_min=minimum(load_arr)
    m_max=maximum(load_arr)
    mc=find_mc(mm,sl,s,g,N,load_arr,m_min,m_max)
    IIC=giant_component(g,s,load_arr,mc,N)
    bw=nx.betweenness_centrality(IIC)
    writedlm(outfbetweenness,bw)
    close(outfbetweenness)
end

function average_betweenness(Pstat=10,Pend=100,ΔP=2,samples=5,h=0.3,L=100,T=10000)
    outfaveragebetweenness=open("/home/SuZ/Data/Priority/ScaleFree_AverageBetweenness_($h)_($L).txt","w")
    #parameters array
    mm=Float64[]
    sl=Float64[]
    s=Int[]
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    sizehint!(s,N)
    P=Pstat
    while P<=Pend
        average_C=0
        for i=1:samples
            #initiate
            package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
            # reach stable state
            reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
            average_C+=average_bw(mm,sl,s,g,N,load_arr)
        end
        writedlm(outfaveragebetweenness,transpose([P,average_C/samples]))
        P+=ΔP
    end
    close(outfaveragebetweenness)
end

function pathlength_and_travelingtime_distribution(P=150,h=0.3,L=100,T=10000)
    #outfile
    outfpathlength=open("/home/SuZ/Data/Priority/ScaleFree_PathLengthDistribution_($P)_($h)_($L).txt","w")
    outftravelingtime=open("/home/SuZ/Data/Priority/ScaleFree_TravelingTimeDistribution_($P)_($h)_($L).txt","w")
    #parameters array
    s=Int[]
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    sizehint!(s,N)
    package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
    # reach stable state
    reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
    writedlm(outfpathlength,values(package_step))
    writedlm(outftravelingtime,values(package_life))
    close(outfpathlength)
    close(outftravelingtime)
end
function average_pathlength_and_travelingtime(Pstat=10,Pend=100,ΔP=2,samples=5,h=0.3,L=100,T=10000)
    outfaveragepathlength=open("/home/SuZ/Data/Priority/ScaleFree_AveragePathlength_($h)_($L).txt","w")
    outfaveragetravelingtime=open("/home/SuZ/Data/Priority/ScaleFree_AverageTravelingtime_($h)_($L).txt","w")
    #parameters array
    mm=Float64[]
    sl=Float64[]
    s=Int[]
    #initiate
    g,N,D,path,pathlength=initiate_graph(L)
    sizehint!(s,N)
    P=Pstat
    while P<=Pend
        average_pathlength=0
        average_travelingtime=0
        for i=1:samples
            #initiate
            package,package_life,package_step,move_package,packagequeue,load_arr,total_load=initiate_package(L,N)
            # reach stable state
            reach_stable_state(g,package,package_life,package_step,packagequeue,move_package,path,load_arr,total_load,h,P,N,T)
            average_pathlength+=mean(values(package_step))
            average_travelingtime+=mean(values(package_life))
        end
        writedlm(outfaveragepathlength,transpose([P,average_pathlength/samples]))
        writedlm(outfaveragetravelingtime,transpose([P,average_travelingtime/samples]))
        P+=ΔP
    end
    close(outfaveragepathlength)
    close(outfaveragetravelingtime)
end


text:
#=
We want to simulate a number of machines that each produce a specific product and
put it in a warehouse (a `Store`). There is a combination process that needs one
of each generated products and combines it into another. The simulation ends when
the fictive container of combined goods is full.

We create our own type `Product` with two fields that allow to identify the kind of
product and to identify its "creator" by means of a serial number.
We also create a type `Machine` that has an associated process. This machine has
an ID an make one type of products. All generated product are put into the same
store.

The combiner process works as follows:
* generate the events matching the availability of one of each product
* wait until all of these events are realised. *Note*: if you only have two events,
you can simply use `ev1` & `ev2`.
* generate the event of putting a fictive combined product in a container
* @yield the event (i.e. time-out until done)
* verify if the container is full and if so stop the simulation on this event.

This application illustrates how you can wait on multiple other events before
continuing the simulation. Keep in mind that this requires ALL events to be
processed.
=#

using Logging
using SimJulia
using ResumableFunctions
import Base.show

struct Product
    kind::Symbol
    serial_number::String
end

show(io::Core.IO, p::Product) =  print(io, "$(p.kind) with serial N° $(p.serial_number)")

mutable struct Machine
    id::String
    prod_kind::Symbol
    prod_time::Int
    n_prod::Int
    prod_store::Store{Product}
    proc::Process
    function Machine(id::String,
                     prod_kind::Symbol, prod_time::Int,
                     prod_store::Store{Product}, n_prod::Int=0)
        m = new()
        m.id = id
        m.prod_kind = prod_kind
        m.prod_time = prod_time
        m.n_prod = n_prod
        m.prod_store = prod_store
        m.proc = @process produce(m.prod_store.env, m)
        return m
    end
end

@resumable function produce(env::Environment, m::Machine)
    while true
        @yield timeout(env, m.prod_time)
        m.n_prod += 1
        newprod = Product(m.prod_kind, "$(m.id)-"*lpad("$(m.n_prod)", 4, "0"))
        @debug "$(env.time) - Machine $(m.id) produced a $(newprod)"
        @yield put(m.prod_store, newprod)
    end
end

@resumable function combiner(env::Environment, prod_store::Store{Product},
                             products::Array{Symbol, 1}, result_container::Container)
    while true
        # obtain a request for each product kind
        requests = map(x::Symbol -> get(prod_store, p::Product -> p.kind == x), products)
        @debug "$(env.time) - Combiner requests: $(requests)"
        # all requests must be matched
        @yield Operator(SimJulia.eval_and, requests...)
        @debug "$(env.time) - Combiner got all required products!"
        newcomb = put(result_container, 1)
        @yield newcomb
        # stop the simulation after the final combination was made
        if result_container.level == result_container.capacity
            @info "Combiner made all products in $(env.time) time units"
            SimJulia.stop_simulation(newcomb)
        end
        @debug "$(env.time) - Current container level $(result_container.level)"
    end

end

function mysim()
    # setup simulation
    @info "\n$("-"^70)\nMachine application\n$("-"^70)\n"
    sim = Simulation()
    # products
    prod_store = Store{Product}(sim)
    products = [:nut, :bolt, :rivet, :beam]
    prod_times = [1;2;4;6]
    # combined results
    result_container = Container(sim, 10)
    # machines
    machines = [Machine("Machine $(i)",products[i], prod_times[i],prod_store) for i in 1:length(products)]
    @process combiner(sim, prod_store, products, result_container)
    run(sim)
end

# logging settings
Logging.disable_logging(LogLevel(-1001)) # (de-)activate debug messages
logger = Logging.ConsoleLogger(stdout, Logging.Debug)

# run simulation
with_logger(logger) do
    mysim()
end


text:
"""
    Read a block of fastq reads into Vector{String}[]
"""
function read_block!(io, data::Vector{Vector{String}}, blocksize=length(data))
    r = 0
    for i = 1:blocksize
        eof(io) && break
        data[i][1] = readline(io)
        data[i][2] = readline(io)
        data[i][3] = readline(io)
        data[i][4] = readline(io)
        r += 1
    end
    r
end

"""
    Write a block of fastq reads from Vector{String}[]
    Writes trimmed sequence and quality string
"""
function write_block(io, data, rn, trim_data)
    for i = 1:rn
        ind = 1:trim_data[2, i]
        write(io, data[i][1], "\n")
        write(io, data[i][2][ind], "\n")
        write(io, data[i][3], "\n")
        write(io, data[i][4][ind], "\n")
    end
end

"""
    Count the number of Ambiguous bases in a sequence
    i.e. count number of N's
"""
count_n(seq) = count(isambiguous, seq)

"""
    Find reverse complement alignment of two sequences
    find
"""
function align_rc(seq1, seq2, align_length=20, max_distance=1, max_n=5)
    rseq2 = reverse_complement(seq2)

    (seq1 == rseq2)         && return 1:0
    rc2 = rseq2[(end-align_length+1):end]
    (count_n(rc2) >= max_n) && return 1:0

    cand_ind = approxrsearch(seq1, rc2, max_distance)

    isempty(cand_ind) && return 1:0
    (cand_ind.stop == length(seq1)) && return 1:0

    return 1:cand_ind.stop
end

"""
    Increment a pwm with a sequence
"""
function add_pwm!(PWM, seq)
    for i = 1:min(size(PWM, 2), length(seq))
        PWM[5, i] += 1
        if seq[i] == DNA_A
            PWM[1, i] += 1
        elseif seq[i] == DNA_C
            PWM[2, i] += 1
        elseif seq[i] == DNA_G
            PWM[3, i] += 1
        elseif seq[i] == DNA_T
            PWM[4, i] += 1
        end
    end
end

"""
    Estimate over respresented sequences in the portion surrounding a reverse complement match between read 1 and 2
"""
function estimate_libtype(fq1, fq2, num_reads=5, rl = 50, sl=0)

    f1 = GzipDecompressorStream(open(fq1))
    f2 = GzipDecompressorStream(open(fq2))

    PWM_L = zeros(5, rl)
    PWM_R = zeros(5, rl)

    nr = 0
    tls = 0
    while !eof(f1) && !eof(f2) && (nr != (num_reads + sl))
        id1   = readline(f1)
        read1 = readline(f1)
        qid1  = readline(f1)
        qs1   = readline(f1)

        id2   = readline(f2)
        read2 = readline(f2)
        qid2  = readline(f2)
        qs2   = readline(f2)

        nr += 1
        (nr <= sl) && continue

        seq1 = LongDNASeq(read1)
        seq2 = LongDNASeq(read2)
        ind = align_rc(seq1, seq2)

        if ind != 1:0
            tls += 1
            adapter1 = seq1[(ind.stop+1):end]
            adapter2 = seq2[(ind.stop+1):end]

            add_pwm!(PWM_L, adapter1)
            add_pwm!(PWM_R, adapter2)
        end
    end

    close(f1)
    close(f2)

    L_consen, L_adapter_type, L_ic = pwm_summary(PWM_L, "[TFR]\tLeft  $(tls) sequences :")
    R_consen, R_adapter_type, R_ic = pwm_summary(PWM_R, "[TFR]\tRight $(tls) sequences :")

    L_consen, R_consen
end

"""
    Get consensus base from a vector of probabilites
"""
function getbase(v, seq=[DNA_A, DNA_C, DNA_G, DNA_T], thresh=0.25)
    m, mi = findmax(v)
    if m < thresh
        return DNA_Gap
    else
        seq[mi]
    end
end

"""
    Get consensus from PWM
"""
consensus(PWM, l=size(PWM, 2)) = LongDNASeq([getbase(view(PWM, 1:4, i)) for i = 1:l])

"""
    Determine is seq `a` is an instance of seq `b` upto mismatches
"""
function isinstance(a, b, mm = 3)
    tmm = 0
    for i = 1:min(length(a), length(b))
        tmm += ifelse(a[i] == b[i], 0, 1)
    end
    tmm ≤ mm
end

"""
    classify_consensus sequences
"""
function classify_consensus(cl)
    isinstance(cl, dna"AGATCGGAAGAGC")        && return :Illumina
    isinstance(cl, dna"CTGTCTCTTATACACATCT")  && return :Nextera
    return :Unknown
end

"""
    Print pwm summary
"""
function pwm_summary(PWM, label="Adapt")
    consen = consensus(PWM)
    adapter_type = classify_consensus(consen)
    ic = info_content(PWM)
    println(label, " consensus: ", consen, "\tType: ", adapter_type, "\tIC: ", ic, " bits")
    consen, adapter_type, ic
end

"""
    Calculate information content of a PWM in bits
"""
function info_content(PWM)
    ts = PWM[5, :]
    en = 3.0./(2*ts*log(2))
    F = PWM[1:4, :]./ts'
    H = -sum(F.*log2.(F), dims=1)
    R = log2(4) .- (H + en')
    sum(filter(r -> !isnan(r) && !isinf(r), R))
end

"""
    Find the adapter location within a sequence upto a mismatch threshold
"""
function find_adapt_location(seq1, seq2, L_match, R_match, thresh=3/13, tail=7)
    n = length(seq1)
    amm = 0
    bmm = 0

    min_p = min(length(L_match), length(R_match))

    for i = 1:(n-tail)
        amm = 0
        bmm = 0

        loop_len = min(min_p, n - i + 1)
        for j = 1:loop_len
            amm += ifelse((seq1[i + j - 1] != L_match[j]), 1, 0)
            bmm += ifelse((seq2[i + j - 1] != R_match[j]), 1, 0)
        end

        if 0.5*(amm + bmm)/loop_len < thresh
            return i - 1, amm, bmm
        end
    end
    return 0, 0, 0
end

"""
    Determine if seqA ≈ rc(seqB) upto mismatches
"""
function rcmm(seqA, seqB)
    rb = reverse_complement(seqB)
    mm = mismatches(seqA, rb) - count(isambiguous, seqA, rb)
    max(0, mm)
end

"""
    Trim Reads function
    Takes pairs of reads and trimm

    ## Trimming procedure
    1. Attempt to find L_match and R_Match adapters in R1 and R2 respectively
    2. Check remaining sequence is RC
    3. If adapter discovery fails attempt to find RC region
    4. Check that remaining region contains adapters

"""
function trim_reads(readR1, readR2, L_match, R_match, align_length, max_distance, mismatch_rate, max_n)

    # Convert string to DNA sequence
    seq1 = LongDNASeq(readR1[2])
    seq2 = LongDNASeq(readR2[2])

    ### Find the adapter location
    adapt_location, r1mm, r2mm = find_adapt_location(seq1, seq2, L_match, R_match, mismatch_rate)
    trim_ind = 1:adapt_location

    if adapt_location > 0 && (rcmm(seq1[trim_ind], seq2[trim_ind]))/length(trim_ind) < mismatch_rate
        return 1, adapt_location
    else
        ### find RC
        trim_ind = align_rc(seq1, seq2, align_length, max_distance, max_n)

        if !isempty(trim_ind)
            adaptind = (trim_ind.stop + 1):length(seq1)

            mm1 = mismatches(seq1[adaptind], L_match)
            mm2 = mismatches(seq1[adaptind], R_match)

            compare_length = min(length(adaptind), length(L_match))
            if (mm1 + mm2)/(2*compare_length) < mismatch_rate
                return 2, trim_ind.stop
            end
        end
    end
    0, length(seq1)
end

"""
    Threaded fastq trimmer

"""
function trim_fastq_threads(fq_R1, fq_R2, align_length=20, max_distance=1, mismatch_rate=1/5, max_n=5,  maxreads=-1, blocksize=100, adapt_length=13, method=:estimate, est_consenus_reads=100000)

    # Set output
    outfile_R1 = replace(fq_R1, ".fastq.gz"   =>  ".trim.fastq.gz")
    outfile_R2 = replace(fq_R2, ".fastq.gz"   =>  ".trim.fastq.gz")
    outfile_R1 = replace(outfile_R1, ".fq.gz" => ".trim.fq.gz")
    outfile_R2 = replace(outfile_R2, ".fq.gz" => ".trim.fq.gz")

    (fq_R1 == outfile_R1) && error("In/Out Equal: $fq_R1 $outfile_R1")
    (fq_R2 == outfile_R2) && error("In/Out Equal: $fq_R2 $outfile_R2")

    ### Print parameters
    println("[TFR]\tFrag align trim       :\t$fq_R1, $fq_R2")
    println("[TFR]\talign_length          :\t", align_length)
    println("[TFR]\tmax_distance          :\t", max_distance)
    println("[TFR]\tmismatch_rate         :\t", mismatch_rate)
    println("[TFR]\tmax_n                 :\t", max_n)
    println("[TFR]\tmaxreads              :\t", maxreads)
    println("[TFR]\tblocksize             :\t", blocksize)
    println("[TFR]\tnum_threads           :\t", Threads.nthreads())
    println("[TFR]\test_consensus_reads   :\t", est_consenus_reads)
    println("[TFR]\tadapt_length          :\t", adapt_length)

    ### Find consensus sequences to trim from
    if method == :estimate # Estimate the libtype from overepresented sequences
        println("[TFR]\tEstimating adapters   :    ")
        L_consen, R_consen = estimate_libtype(fq_R1, fq_R2, est_consenus_reads)
    elseif method == :Illumina
        println("[TFR]\tAdapter Search        :\t", method)
        L_consen = dna"AGATCGGAAGAGC"
        R_consen = dna"AGATCGGAAGAGC"
    elseif method == :Nextera
        println("[TFR]\tAdapter Search        :\t", method)
        L_consen = dna"CTGTCTCTTATACACATCT"
        R_consen = dna"CTGTCTCTTATACACATCT"
    else
        println("[TFR]\tEstimating adapters   :    ")
        L_consen, R_consen = estimate_libtype(fq_R1, fq_R2, est_consenus_reads)
    end

    L_match = L_consen[1:min(length(L_consen), adapt_length)]
    R_match = R_consen[1:min(length(R_consen), adapt_length)]

    println("[TFR]\tL_match               :\t", L_match)
    println("[TFR]\tR_match               :\t", R_match)
    println("[TFR]\tWriting               :\t$outfile_R1, $outfile_R2")
    flush(stdout) ### Flushes IO so can read the output parameters whilst trimming
    starttime = time()

    #### Open readers and writers
    reader_r1 = open(fq_R1) |> GzipDecompressorStream
    reader_r2 = open(fq_R2) |> GzipDecompressorStream

    writer_r1 = open(outfile_R1, "w") |> GzipCompressorStream
    writer_r2 = open(outfile_R2, "w") |> GzipCompressorStream

    #### Datastructure to read fastq files and trim them
    datablock_R1 = [Vector{String}(undef, 4) for i = 1:blocksize]
    datablock_R2 = [Vector{String}(undef, 4) for i = 1:blocksize]

    ### Trimming stats
    temp_trim_data = Vector{Tuple{Int, Int}}(undef, blocksize)
    trim_data  = zeros(Int, 2, blocksize)
    trim_stats = counter(Tuple{Int, Int})

    total_reads = 0
    max_trim_length = 0

    while true
        ### Read block of reads
        rn1 = read_block!(reader_r1, datablock_R1, blocksize)
        rn2 = read_block!(reader_r2, datablock_R2, blocksize)
        (rn1 != rn2) && error("FASTQ mismatch")
        (rn1 == 0) && break

        ### Threaded trimming of the block
        #Threads.@threads
        for i = 1:rn1
            trim_data[1, i], trim_data[2, i] = trim_reads(datablock_R1[i], datablock_R2[i], L_match, R_match, align_length, max_distance, mismatch_rate, max_n)
        end

        #### Accumulate trim stats
        for i = 1:rn1
            if trim_data[1, i] > 0
                max_trim_length = max(trim_data[2, i], max_trim_length)
            end
            push!(trim_stats, (trim_data[1, i], trim_data[2, i]))
        end

        #### Write Reads
        write_block(writer_r1, datablock_R1, rn1, trim_data)
        write_block(writer_r2, datablock_R2, rn1, trim_data)

        total_reads += rn1
        (maxreads != -1) && (total_reads > maxreads) && break
    end

    ### Close readers and writers
    close(reader_r1)
    close(reader_r2)
    close(writer_r1)
    close(writer_r2)

    ### Print stats
    total_trimmmed = 0

    F = zeros(Int, 3, max_trim_length)

    for ((trim_type, frag), trim_count) ∈ trim_stats
        (trim_type == 0) && continue
        total_trimmmed += trim_count

        F[1, frag]             += trim_count
        F[trim_type + 1, frag] += trim_count
    end

    println("[TFR]\tComplete in           :\t", time() - starttime, " seconds")
    println("[TFR]\ttotal_reads           :\t", total_reads)
    println("[TFR]\ttotal_trimmed         :\t", total_trimmmed, "\t", total_trimmmed/total_reads)
    println("[TFR]\tTrimmed frag dist     :\t")

    for i = 1:max_trim_length
        println(i, "\t", F[1, i], "\t", F[2, i], "\t", F[3, i])
    end

end


text:
# dummy type wrapping a Float64 used in tests
struct F64 <: AbstractFloat
    x::Float64
end
F64(x::F64) = x

# operations
for op in (:+, :-, :sin, :cos, :asin, :acos)
    @eval Base.$op(a::F64) = F64($op(a.x))
end
for op in (:+, :-, :*, :/, :rem, isdefined(Base, :atan2) ? :atan2 : :atan)
    @eval Base.$op(a::F64, b::F64) = F64($op(a.x, b.x))
end
for op in (:zero, :one,)
    @eval Base.$op(::Type{F64}) = F64($op(Float64))
end

if VERSION >= v"0.7-"
    Random.rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{F64}}) = F64(rand(rng))
else
    Base.rand(rng::AbstractRNG, ::Type{F64}) = F64(rand())
end
Base.sqrt(a::F64) = F64(sqrt(a.x))
Base.:^(a::F64, b::Number) = F64(a.x^b)
Base.:^(a::F64, b::Int) = F64(a.x^b)
Base.:^(a::F64, b::F64) = F64(a.x^b.x)
Base.:^(a::Number, b::F64) = a^b.x
Base.log(a::F64) = F64(log(a.x))
Base.isfinite(a::F64) = isfinite(a.x)
Base.float(a::F64) = a.x
Base.rtoldefault(a::Type{F64}, b::Type{F64}) = Base.rtoldefault(Float64, Float64)
Base.mod(a::F64, b::F64) = mod(a.x, b.x)
# comparison
Base.isapprox(a::F64, b::F64) = isapprox(a.x, b.x)
Base.:<(a::F64, b::F64) = a.x < b.x
Base.:<=(a::F64, b::F64) = a.x <= b.x
Base.eps(::Type{F64}) = eps(Float64)

# promotion
Base.promote_type(::Type{Float32}, ::Type{F64}) = Float64 # for eig
Base.promote_type(::Type{Float64}, ::Type{F64}) = Float64 # for vecnorm
Base.promote(a::F64, b::T) where {T <: Number} = a, F64(float(b))
Base.promote(a::T, b::F64) where {T <: Number} = F64(float(a)), b

Base.convert(::Type{F64}, a::F64) = a
Base.convert(::Type{Float64}, a::F64) = a.x
Base.convert(::Type{F64}, a::T) where {T <: Number} = F64(a)

# conversion
Base.Float64(a::F64) = Float64(a.x)
Base.Int64(a::F64) = Int64(a.x)
Base.Int32(a::F64) = Int32(a.x)


text:
export print_histogram

# This is essentially Base.show(io::IO, ::MIME"text/plain", t::Trial)
# in BenchmarkTools.jl

function bindata(sorteddata, nbins, min, max)
    Δ = (max - min) / nbins
    bins = zeros(nbins)
    lastpos = 0
    for i ∈ 1:nbins
        pos = searchsortedlast(sorteddata, min + i * Δ)
        bins[i] = pos - lastpos
        lastpos = pos
    end
    return bins
end

function asciihist(bins, height=1)
    histbars = ['▁', '▂', '▃', '▄', '▅', '▆', '▇', '█']
    if minimum(bins) == 0
        barheights = 2 .+ round.(Int, (height * length(histbars) - 2) * bins ./ maximum(bins))
        barheights[bins .== 0] .= 1
    else
        barheights = 1 .+ round.(Int, (height * length(histbars) - 1) * bins ./ maximum(bins))
    end
    heightmatrix = [min(length(histbars), barheights[b] - (h-1) * length(histbars))
                    for h in height:-1:1, b in 1:length(bins)]
    return map(height -> if height < 1; ' ' else histbars[height] end, heightmatrix)
end

function print_histogram(v::Vector{T}) where {T<:AbstractFloat}
    medv = median(v)
    avgv = mean(v)
    stdv = std(v)
    minv = minimum(v)
    maxv = maximum(v)

    medstr = string(@sprintf("%.3f", medv))
    avgstr = string(@sprintf("%.3f", avgv))
    stdstr = string(@sprintf("%.3f", stdv))
    minstr = string(@sprintf("%.3f", minv))
    maxstr = string(@sprintf("%.3f", maxv))

    lmaxwidth = maximum(length.((medstr, avgstr, minstr)))
    rmaxwidth = maximum(length.((stdstr, maxstr)))

    print(" Range ")
    printstyled("("; color=:light_black)
    printstyled("min"; color=:cyan, bold=true)
    print(" … ")
    printstyled("max"; color=:magenta)
    printstyled("):  "; color=:light_black)
    printstyled(lpad(minstr, lmaxwidth); color=:cyan, bold=true)
    print(" … ")
    printstyled(lpad(maxstr, rmaxwidth); color=:magenta)
    print("  ")
    printstyled("┊"; color=:light_black)

    print("\n Value " )
    printstyled("("; color=:light_black)
    printstyled("median"; color=:blue, bold=true)
    printstyled("):     "; color=:light_black)
    printstyled(lpad(medstr, lmaxwidth), rpad(" ", rmaxwidth + 5); color=:blue, bold=true)
    printstyled("┊"; color=:light_black)

    print("\n Value " )
    printstyled("("; color=:light_black)
    printstyled("mean"; color=:green, bold=true)
    print(" ± ")
    printstyled("σ"; color=:green)
    printstyled("):   "; color=:light_black)
    printstyled(lpad(avgstr, lmaxwidth); color=:green, bold=true)
    print(" ± ")
    printstyled(lpad(stdstr, rmaxwidth); color=:green)
    print("  ")
    printstyled("┊"; color=:light_black)

    # The height and width of the printed histogram in characters.
    histheight = 8
    histwidth = 42 + lmaxwidth + rmaxwidth

    histv = sort(v)

    bins = bindata(histv, histwidth - 1, minv, maxv)
    # if median size of (bins with >10% average data/bin) is less than 5% of max bin size, log the bin sizes
    if median(filter(b -> b > 0.1 * length(v) / histwidth, bins)) / maximum(bins) < 0.05
        bins, logbins = log.(1 .+ bins), true
    else
        logbins = false
    end
    hist = asciihist(bins, histheight)
    hist[:,end-1] .= ' '
    maxbin = maximum(bins)

    delta1 = (maxv - minv) / (histwidth - 1)
    if delta1 > 0
        medpos = 1 + round(Int, (histv[length(v) ÷ 2] - minv) / delta1)
        avgpos = 1 + round(Int, (avgv - minv) / delta1)
    else
        medpos, avgpos = 1, 1
    end

    print("\n")
    for r in axes(hist, 1)
        print("\n  ")
        for (i, bar) in enumerate(view(hist, r, :))
            color = :default
            if i == avgpos color = :green end
            if i == medpos color = :blue end
            printstyled(bar; color=color)
        end
    end

    print("\n  ", minstr)
    caption = "Histogram: " * ( logbins ? "log(frequency)" : "frequency" )
    if logbins
        printstyled(" " ^ ((histwidth - length(caption)) ÷ 2 - length(minstr)); color=:light_black)
        printstyled("Histogram: "; color=:light_black)
        printstyled("log("; bold=true, color=:light_black)
        printstyled("frequency"; color=:light_black)
        printstyled(")"; bold=true, color=:light_black)
    else
        printstyled(" " ^ ((histwidth - length(caption)) ÷ 2 - length(minstr)), caption; color=:light_black)
    end
    print(lpad(maxstr, ceil(Int, (histwidth - length(caption)) / 2) - 1), " ")
    printstyled("<"; bold=true)
    print("\n")
    return nothing
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# workspace()
using Indicators
using PyPlot
using Base.Dates
using Temporal

aapl = yahoo("AAPL")
aapl = aapl["2015"]

t = aapl.index
aapl = aapl.values
op = aapl[:,1]
hi = aapl[:,2]
lo = aapl[:,3]
cl = aapl[:,4]
vo = aapl[:,5]
hlc = [hi lo cl]

subplot(411)
plot(t, cl, lw=2, c="k", label="AAPL")
plot(t, kama(cl), c="b", label="Kaufman AMA")
plot(t, trima(cl), c="g", label="Triangula MA")
plot(t, hma(cl), c="r", label="Hull MA")
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

subplot(412)
plot(t, kst(cl), c="m", label="KST")
plot(t, sma(kst(cl),n=9), c="c", label="Signal")
plot([t[1],t[end]], [0,0], ls="--", c=[0.4,0.4,0.4])
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

subplot(413)
plot(t, wpr(hlc), c=[1,0.5,0], label="Williams %R")
plot([t[1],t[end]], [-20,-20], c="r", ls="--")
plot([t[1],t[end]], [-80,-80], c="g", ls="--")
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

subplot(414)
plot(t, cci(hlc), c="c", label="CCI")
plot([t[1],t[end]], [-100,-100], c="g", ls="--")
plot([t[1],t[end]], [100,100], c="r", ls="--")
grid(ls="-", c=[0.8,0.8,0.8])
legend(loc="best", frameon=false)

tight_layout()


text:
import Pkg
Pkg.activate("..")

using Distributions: Normal, LogNormal, truncated
using CrimsonSkyline
using Random: seed!
using Logging
using StatsBase: mean, std
using PrettyPrint: pprintln

seed!(2021)
const Maybe{T} = Union{T, Nothing}

# version for inference -- no return value, so lower memory reqd

function random_walk(t :: Trace, data :: Array{Float64, 2})
    loc = sample(t, :loc, Normal())
    scale = sample(t, :scale, LogNormal())
    shape = size(data)
    for n in 1:shape[1]
        x_n = sample(t, (:ic, n), Normal())
        for t_ix in 1:shape[2]
            x_n = observe(t, (:data, n, t_ix), Normal(loc + x_n, scale), data[n, t_ix])
        end
    end
end

# version for prediction and forecasting that returns sampled paths

function random_walk(t :: Trace, shape :: Tuple{Int64, Int64})
    loc = sample(t, :loc, Normal())
    scale = sample(t, :scale, LogNormal())
    out = Vector{Vector{Float64}}()
    for n in 1:shape[1]
        path = Float64[]
        x_n = sample(t, (:ic, n), Normal())
        push!(path, x_n)
        for t_ix in 1:shape[2]
            x_n = sample(t, (:data, n, t_ix), Normal(loc + x_n, scale))
            push!(path, x_n)
        end
        push!(out, path)
    end
    out
end

function ic_proposal_factory(data :: Array{Float64, 2})
    proposals = []
    n = size(data, 1)
    for k in 1:n
        push!(proposals, (t0, t1, data) -> propose(t1, (:ic, k), Normal(t0[(:ic, k)].value, 1.0)))
    end
    proposals
end
loc_proposal(t0, t1, data) = propose(t1, :loc, Normal(t0[:loc].value, 1.0))
scale_proposal(t0, t1, data) = propose(t1, :scale, truncated(Normal(t0[:scale].value, 1.0), 0.0, Inf))

function fake_data()
    @info "Generating fake data"
    n = 2
    t = 20
    t_pred = 10

    loc = 1.0
    scale = 1.5
    noise = randn(n, t + t_pred)
    data = cumsum(loc .+ scale .* noise, dims=2)
    test_data = data[:, t:end]
    data = data[:, 1:t]
    ic_proposals = ic_proposal_factory(data)
    proposals = append!([loc_proposal, scale_proposal], ic_proposals)

    @info "True loc = $loc, true scale = $scale"
    @time results = mh(random_walk, proposals; params=(data,), burn=250, thin=20, num_iterations=3000)
    @info "Estimated loc = ($(mean(results, :loc)) ± $(2.0 * std(results, :loc)))"
    @info "Estimated scale = ($(mean(results, :scale)) ± $(2.0 * std(results, :scale)))"

    n_forecast = 10
    posterior_model = update(random_walk, results)
    evidence = Dict((:ic, k) => data[k, end] for k in 1:n)
    forecast_model = condition(posterior_model, evidence)
    forecast = [forecast_model(trace(), (n, t_pred))[end][end] for _ in 1:n_forecast]
    @info "Forecasts: $forecast"
end

fake_data()