text:
# represent the   basis x^(mbar) = x⋅(x-1)⋯(x-m+1) ∈ Πm
"""
    FallingFactorial{T}

Construct  a  polynomial with   respect to the basis `x⁰̲,  x¹̲, x²̲, …` where
`xⁱ̲ = x  ⋅  (x-1) ⋅  (x-2)  ⋯ (x-i+1)` is the falling Pochhammer  symbol.  See
[Falling factorial](https://en.wikipedia.org/wiki/Falling_and_rising_factorials)  for several  facts
about this  polynomial basis.

In [Koepf and Schmersau](https://arxiv.org/pdf/math/9703217.pdf)
connection coefficients between the falling factorial polynomial
system and classical discrete orthogonal polynomials are given.

## Examples

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> p = basis(FallingFactorial, 3)
FallingFactorial(1.0⋅x³̲)

julia> x = variable(Polynomial)
Polynomials.Polynomial(1.0*x)

julia> p(x) ≈ x*(x-1)*(x-2)
true
```

"""
struct FallingFactorial{T <: Number,X} <: AbstractSpecialPolynomial{T,X}
    coeffs::Vector{T}
    function FallingFactorial{T,X}(coeffs::Vector{T}) where {T <: Number,X}
        length(coeffs) == 0 && return new{T,X}(zeros(T, 1))
        last_nz = findlast(!iszero, coeffs)
        last = max(1, last_nz === nothing ? 0 : last_nz)
        return new{T,X}(coeffs[1:last])
    end
    function FallingFactorial{T}(coeffs::Vector{T}, var::Polynomials.SymbolLike=:x) where {T <: Number}
        FallingFactorial{T,Symbol(var)}(coeffs)
    end

end

Polynomials.@register FallingFactorial
export  FallingFactorial

# modified from  Polynomials.unicode_exponent
function unicode_exponent(io, var, j)
    print(io, var)
    a = ("⁻","","","⁰","¹","²","³","⁴","⁵","⁶","⁷","⁸","⁹")
    for i in string(j)
        print(io, a[Int(i)-44])
    end
end

# show as 1⋅(x)₀ + 2⋅(x)₁ + 3⋅(x)₂
function Polynomials.showterm(io::IO, ::Type{P}, pj::T, var, j, first::Bool, mimetype) where {N, T, P <: FallingFactorial}
    iszero(pj) && return false
    !first &&  print(io, " ")
    print(io, Polynomials.hasneg(T)  && Polynomials.isneg(pj) ? "- " :  (!first ? "+ " : ""))
    print(io, "$(abs(pj))⋅")
#    print(io,"$(var)")
    unicode_exponent(io, var, j)
    print(io,"̲")
    return true
end

function  Polynomials.evalpoly(x, p::FallingFactorial)
    d = degree(p)
    d <= 0 &&  return  p[0]*one(x)

    xⁱ̲ =  one(x)
    tot =  p[0]*xⁱ̲
    for i in  1:d
        xⁱ̲ *= x-(i-1)
        tot = muladd(xⁱ̲, p[i], tot)
    end

    tot
end

function Base.one(::Type{P}) where {P<:FallingFactorial}
    T,X = eltype(P), Polynomials.indeterminate(P)
    ⟒(P){T,X}(ones(T,1))
end

function Polynomials.variable(::Type{P}) where {P<:FallingFactorial}
    T,X = eltype(P), Polynomials.indeterminate(P)
    ⟒(P){T,X}([zero(T),one(T)])
end

k0(P::Type{<:FallingFactorial}) = one(eltype(P))

Base.:*(p::FallingFactorial, q::FallingFactorial) = convert(FallingFactorial,  convert(Polynomial,p)*convert(Polynomial,q))

## Connect FallingFactorial with AbstractCDOP
function Base.convert(::Type{P}, q::Q) where {
    P <: FallingFactorial,
    Q <: AbstractCDOP}
    _convert_cop(P,q)
end
function Base.convert(::Type{P}, q::Q) where {
    P <: AbstractCDOP,
    Q <: FallingFactorial}
    _convert_cop(P,q)
end

##  Koepf and Schmersa Thm 6 connection for P, Q=x^n̲
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: AbstractCDOP,
    Q <: FallingFactorial}

    a,b,c,d,e = abcde(P)

    c₀ = (a*n + a*m - a + d) * (n - m)

    c₁ = (m+1) * (a*n^2 - 2a*m^2 - a*n - a*m  + n*d - 2d*m  - b*m - d - e)

    c₂ = -(m+1) * (m+2) * (a*m^2 +  2a*m + d*m + b*m + a + d + b + c + e)

    (c₀,c₁,c₂)
end

##  Koepf and Schmersa Thm 7, p25 connection for P=x^n̲, Q
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: FallingFactorial,
    Q <: AbstractCDOP}

    ā,b̄,c̄,d̄,ē = abcde(Q)

    c₀ =  (2m*ā + ā + d̄)  * (2m*ā +  3ā + d̄) * (2m*ā + 2ā + d̄)^2 * (2m*ā + d̄) * (n - m)

    c₁ =  (2m*ā + ā + d̄)  * (2m*ā +  3ā + d̄) * (2m*ā + 2ā + d̄)  * (m+1)
    c₁a = 2m^2*n*ā^2  - 2m^2*ā^2
    c₁a += m^2*ā*d̄ + 2m^2*ā*b̄ + 2m*n*ā^2 + 2m*n*ā*d̄ - 2m*ā^2 - m*ā*d̄ + 2m*ā*b̄ + m*d̄^2
    c₁a += 2m*d̄*b̄ + n*ā*d̄ +  2n*ā*ē - n*d̄*b̄ -  ā*d̄ + d̄*b̄  + d̄*ē
    c₁ *= c₁a

    c₂ = (m+1) * (2m*ā + d̄)
    c₂a = m^4*ā^3 + 4m^3*ā^3 + 2m^3*ā^2*d̄ + 6m^2*ā^3 + 6m*ā^2*d̄
    c₂a += 4m^2*ā^2*c̄ + 2m^2*ā^2*ē + m^2*ā*d̄^2 - m^2*ā*d̄*b̄ - m^2*ā*b̄^2 + 4m*ā^3 + 6m^2*ā^2*d̄
    c₂a += 8m*ā^2*c̄ + 4m*ā^2*ē + 2m*ā*d̄^2 - 2m*ā*d̄*b̄ + 4m*ā*d̄*c̄ + 2m*ā*d̄*ē
    c₂a += -2m*ā*b̄^2 - m*d̄^2*b̄ - m*d̄*b̄^2 + ā^3 + 2ā^2*d̄ +  4ā^2*c̄ + 2ā^2*ē +  ā*d̄^2
    c₂a += -ā*d̄*b̄  + 4ā*d̄*c̄ + 2*ā*d̄*ē - ā*b̄^2 + ā*ē^2 - d̄^2*b̄ + d̄^2*c̄
    c₂a += -d̄*b̄^2 - d̄*b̄*ē
    c₂ *= c₂a * (m+2) *  (m*ā + n*ā + ā + d̄)

    (c₀,c₁,c₂)
end

function Base.convert(::Type{P}, q::Q) where {
    P <: Polynomials.StandardBasisPolynomial,
    Q <: FallingFactorial}
    q(variable(P))
end
function Base.convert(::Type{P}, q::Q) where {
    P <: FallingFactorial,
    Q <: Polynomials.StandardBasisPolynomial}
    connection(P,q)
end

# stirling of 2nd kind
@memoize function sterling²(n::Int, k::Int)
    (iszero(n) && iszero(k)) &&  return 1
    (iszero(n) || iszero(k)) && return 0

    k*sterling²(n-1,k) + sterling²(n-1, k-1)
end

# xʲ = ∑ᵢ sterling²(j,i) (x)ᵢ [wikipedia](https://en.wikipedia.org/wiki/Falling_and_rising_factorials)
function Base.iterate(o::Connection{P, Q}, state=nothing) where {P<:FallingFactorial, Q<:Polynomials.StandardBasisPolynomial}
    n,k = o.n, o.k
    if state  == nothing
        j = k
    else
        j = state
        j += 1
    end
    j > n && return nothing
    val = sterling²(j,k)
    (j,val), j
end


text:
"""
PyQuante in Julia
Experimenting with writing quantum chemistry in Julia

"""

function test()
  test_utils()
  test_pgbf()
  test_cgbf()
  test_overlap()
  test_kinetic()
  test_a_terms()
  test_gamma()
  test_na()
  test_fgamma()
  test_one()
  test_na2()
  test_two_terms()
  test_coul1()
  test_vrr()
  test_hrr()
  test_geo_basis()
  test_h2()
  test_lih()
  test_h2o()
end

test()

ext: jl


text:
"""
    struct DiagnosticIDs <: Any
Hold a pointer to a `clang::DiagnosticIDs` object.
"""
struct DiagnosticIDs
    ptr::CXDiagnosticIDs
end

Base.unsafe_convert(::Type{CXDiagnosticIDs}, x::DiagnosticIDs) = x.ptr
Base.cconvert(::Type{CXDiagnosticIDs}, x::DiagnosticIDs) = x


text:
### MCTunerState subtypes hold the samplers' temporary output used for tuning the sampler

abstract type MCTunerState end

mutable struct BasicMCTune <: MCTunerState
  step::Real # Stepsize of MCMC iteration (for ex leapfrog in HMC or drift stepsize in MALA)
  accepted::Integer # Number of accepted MCMC samples during current tuning period
  proposed::Integer # Number of proposed MCMC samples during current tuning period
  totproposed::Integer # Total number of proposed MCMC samples during burnin
  rate::Real # Observed acceptance rate over current tuning period

  function BasicMCTune(step::Real, accepted::Integer, proposed::Integer, totproposed::Integer, rate::Real)
    @assert (step > 0 || isnan(step)) "Stepsize of MCMC iteration should be positive or not a number (NaN)"
    @assert accepted >= 0 "Number of accepted MCMC samples should be non-negative"
    @assert proposed >= 0 "Number of proposed MCMC samples should be non-negative"
    @assert totproposed >= 0 "Total number of proposed MCMC samples should be non-negative"
    if !isnan(rate)
      @assert 0 <= rate <= 1 "Observed acceptance rate should be in [0, 1]"
    end
    new(step, accepted, proposed, totproposed, rate)
  end
end

BasicMCTune(step::Real=1., accepted::Integer=0, proposed::Integer=0, totproposed::Integer=0) =
  BasicMCTune(step, accepted, proposed, totproposed, NaN)

function reset_burnin!(tune::MCTunerState)
  tune.totproposed += tune.proposed
  (tune.accepted, tune.proposed, tune.rate) = (0, 0, NaN)
end

rate!(tune::MCTunerState) = (tune.rate = tune.accepted/tune.proposed)

### Root Monte Carlo Tuner

abstract type MCTuner end


text:
using StaticArrays, Distributions, ProgressMeter

export neutral_growth, neutral_growth!, q

"""
    Haplotype(mutations, n, t_birth, parent)

Type used by the `neutral_growth` simulation. Fields are `mutations` an ordered list of mutations defining the haplotype,
`n` the number of cells of that type, `t_birth` the time of birth of the first cell of this type and `parent` the rowindex of the parental haplotype
within the tumor Array.
"""
mutable struct Haplotype
    mutations :: SVector{L, Int} where L
    n :: Int
    t_birth :: Float64
    parent :: Int
end

"""
    neutral_growth!(tumor, obs, endsize, T=Inf; b, d, μ, t=0., Nthresh=endsize, Tthresh = Inf )
    neutral_growth!(tumor, endsize, T=Inf; b, d, μ, t=0., Nthresh=endsize, Tthresh = Inf )

Neutral (exponential) growth simulation of `Haplotypes`. Receives a Vector of `Haplotypes`, optionally an according `obs` Matrix, an `endsize` and final time `T` value.
If an `obs` array is handed over the following observables are recorded during recording:

* obs[m,1] records the parents haplotype n at birth of mutation m

* obs[m,2] records the sum of n of all haplotypes containing the parents last mutation while not
         containing m

         this includes estranged parents while avoiding offspring of m

* obs[m,3] records tumor size at birth of m

Kwargs are `b` birth-rate, `d` death-rate, `μ` mutation-rate, `t` current time, `Nthresh` and `Tthresh` size and time after which `μ` is set to zero.

Output is a named Tuple `(:tumor, :obs, :N, :t)` containing the final `Haplotype` Vector, observables, endsize and final time.
"""
function neutral_growth!(tumor::Vector{Haplotype}, obs::Matrix, endsize::Int, T=Inf; b, d, μ, t=0., Nthresh=endsize, Tthresh = T, showprogress=true )

    mutID = length(tumor)-1
    obsarr = collect.(eachrow(obs))
    obstemp = Vector{Int}(undef, size(obs,2))
    F = cumsum(getfield.(tumor,:n))
    N = last(F)

    if showprogress prog = ProgressUnknown("Progress: N ") end
    while 0 < N < endsize && t < T
        μp = N < Nthresh && t < Tthresh ? μ : 0.

        p = rand()*(b+d)
        m = rand(1:N)
        row = findfirst(i-> i >= m, F)
        parent = tumor[row]

        t += rand(Exponential(1.))/((b+d)*N)
        if p > d

            mnew =  (rand()<μp/2) + (rand()<μp/2)
            for n=1:mnew
                mutID += 1
                push!(tumor, Haplotype(push(parent.mutations, mutID), 1, t, row) )
                push!(F, N+n)

                obstemp[1] = parent.n + 1 - mnew
                obstemp[2] = isempty(parent.mutations) ? N : sum( getfield.(filter(h-> last(parent.mutations) in h.mutations, tumor[row : end-n]),:n) )
                obstemp[3] = N
                push!(obsarr, copy(obstemp) )
            end
            if mnew == 0
                parent.n += 1
                F[row:end] .+= 1
            elseif mnew == 2
                if isone(parent.n) && !isone(row)
                    deleteat!(tumor, row)
                    deleteat!(F, row)
                    deleteat!(obsarr, row-1)
                else
                    parent.n -= 1
                end
                F[row:end] .-= 1
            end
            N += 1

        else
            if isone(parent.n) && !isone(row)
                deleteat!(tumor, row)
                deleteat!(F, row)
                deleteat!(obsarr, row-1)
            else
                parent.n -= 1
            end
            F[row:end] .-= 1
            N -= 1
        end

        if showprogress ProgressMeter.update!(prog, N) end
    end
    obs = Matrix{Int}(undef, length(obsarr), size(obs,2))
    for i=1:size(obs,2)
        obs[:, i] = getindex.(obsarr, i)
    end

    return (tumor=tumor, obs=obs, N=N, t=t)
end

function neutral_growth!(tumor::Vector{Haplotype}, endsize::Int, T=Inf; b, d, μ, t=0., Nthresh=endsize, Tthresh = T, showprogress=true )

    mutID = length(tumor)-1
    F = cumsum(getfield.(tumor,:n))
    N = last(F)

    if showprogress prog = ProgressUnknown("Progress: N ") end
    while 0 < N < endsize && t < T
        μp = N < Nthresh && t < Tthresh ? μ : 0.

        p = rand()*(b+d)
        m = rand(1:N)
        row = findfirst(i-> i >= m, F)
        parent = tumor[row]

        t += rand(Exponential(1.))/((b+d)*N)
        if p > d

            mnew =  (rand()<μp/2) + (rand()<μp/2)
            for n=1:mnew
                mutID += 1
                push!(tumor, Haplotype(push(parent.mutations, mutID), 1, t, row) )
                push!(F, N+n)
            end
            if mnew == 0
                parent.n += 1
                F[row:end] .+= 1
            elseif mnew == 2
                if isone(parent.n)
                    deleteat!(tumor, row)
                    deleteat!(F, row)
                else
                    parent.n -= 1
                end
                F[row:end] .-= 1
            end
            N += 1

        else
            if isone(parent.n)
                deleteat!(tumor, row)
                deleteat!(F, row)
            else
                parent.n -= 1
            end
            F[row:end] .-= 1
            N -= 1
        end

        if showprogress ProgressMeter.update!(prog, N) end
    end

    return (tumor=tumor, N=N, t=t)

end

"""
    neutral_growth(endsize, T; return_obs::Bool, args... )

Neutral (exponential) growth simulation of `Haplotypes`. Receives `endsize` and final time `T`. See documentation of `neutral_growth!`
for kwargs.

Output is a named Tuple `(:tumor, :N, :t)` containing the final `Haplotype` Vector, endsize and final time.

If `return_obs` is `true` observables are recorded during simulation and additionally returned under the key `:obs`.
"""
function neutral_growth(endsize, T=Inf; return_obs::Bool, args... )
    mutID = 0
    tumor = [Haplotype(SVector{0, Int}(), 1, 0., 0 )]
    obs = Matrix{Int}(undef, 0, 3)

    out =  return_obs ? neutral_growth!(tumor, obs, endsize, T; args... ) : neutral_growth!(tumor, endsize, T; args... )
    return iszero(out[:N]) ? neutral_growth(endsize, T; return_obs=return_obs, args... ) : out
end

"""
    q(d; b, μ )
Function returning the extinction probability for given `d` death-rate, `b` birth-rate and `μ` mutation-rate.

**Branchingprocess `v2`**: independent mutation of offspring pair at division. `μ` is the probability of mutation per cell division.
The mutation probability per cell division per cell is `μ/2`.

"""
q(d; b, μ ) = (d + b*(μ/2)^2) / (b*(1-(μ/2))^2)

;

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
module Tokens

using Compat
using AbstractTrees

import AbstractTrees: children, AbstractShadowTree, make_zip, first_tree, second_tree
import Base: convert

export ¬, ⨳, ⪥, ⤄, √, AbstractToken, Token, SourceLocToken, SourceRange

abstract AbstractToken

const ASTVerbatim = Union{Symbol,String,LineNumberNode,Integer,Void,Char,Bool,AbstractFloat}
const ASTExprs = Union{Expr, QuoteNode, GlobalRef}

# Defensive definitions
Base.isequal(x::AbstractToken, y::ASTVerbatim) = error("Comparing token to raw AST Node")

global √

immutable Token <: AbstractToken
  val::ASTVerbatim
end
val(t::Token) = t.val
val(t::Union{ASTVerbatim,ASTExprs}) = t
√(t::Union{ASTVerbatim, ASTExprs, Expr, Symbol}) = SourceRange()
√(t::Token) = SourceRange()

immutable SourceRange
    offset::UInt32
    length::UInt32
    file::UInt8
end
SourceRange() = SourceRange(-1 % UInt32,-1 % UInt16,-1 % UInt8)

immutable SourceLocToken <: AbstractToken
  val::ASTVerbatim
  loc::SourceRange
end
SourceLocToken(tok::SourceLocToken) = tok
SourceLocToken(val, offset, length, file) =
    SourceLocToken(val, SourceRange(offset, length, file))
SourceLocToken(val) = SourceLocToken(val, SourceRange())
val(t::SourceLocToken) = t.val

const ¬ = val
√(tok::SourceLocToken) = tok.loc
√(r::SourceRange) = r

# Expr and source tree Construction
immutable SourceNode
    loc::Any # SourceRange
    children::Vector
end
Base.copy(x::SourceNode) = SourceNode(x.loc,copy(x.children))
SourceNode(_::Void) = SourceNode(SourceRange(),SourceNode[])
convert(::Type{SourceNode},loc::SourceRange) = SourceNode(loc,SourceNode[])
SourceNode(x::SourceNode) = copy(x)
children(node::SourceNode) = node.children
Base.show(io::IO, node::SourceNode) = print(io, node.loc)

immutable SourceExpr <: AbstractShadowTree
    expr::Any
    loc::SourceNode
end
make_zip(x::SourceExpr) = AbstractTrees.zip_min(children(x.expr),children(x.loc))
first_tree(x::SourceExpr) = x.expr
second_tree(x::SourceExpr) = x.loc
val(expr::SourceExpr) = expr.expr
√(expr::SourceExpr) = expr.loc

normalize(loc::SourceNode) = normalize(loc.loc)
normalize(x) = x
typealias NodeOrRange Union{SourceNode,SourceRange}
function merge(loc1::NodeOrRange, loc2::NodeOrRange)
    loc1, loc2 = normalize(loc1), normalize(loc2)
    loc1 == SourceRange() && return loc2
    loc2 == SourceRange() && return loc1
    @assert loc1.file == loc2.file
    start = min(loc1.offset, loc2.offset)
    stop = max(loc1.offset + loc1.length, loc2.offset + loc2.length)
    res = SourceRange(start, stop-start, loc1.file)
    res
end

const ⤄ = merge
function ⤄(ex::SourceExpr, x::SourceRange)
    SourceExpr(ex.expr,SourceNode(ex.loc.loc ⤄ x, ex.loc.children))
end
⤄(ex::SourceExpr, x::SourceExpr) = ⤄(ex, √x)
⤄(ex::SourceExpr, x::SourceNode) = ⤄(ex, x.loc)
⤄(x::Void, y::SourceRange) = SourceExpr(x,SourceNode(y))
⤄(ex::Union{ASTVerbatim, ASTExprs}, x::NodeOrRange) = SourceExpr(ex,SourceNode(x))
⤄(ex::Union{ASTVerbatim, ASTExprs}, x::Union{SourceExpr,SourceLocToken}) = SourceExpr(ex,SourceNode(normalize(√x)))
⤄(ex::ASTVerbatim, x::Union{SourceExpr,SourceLocToken}) = SourceLocToken(ex,normalize(√x))
⤄(ex::SourceExpr, x::SourceLocToken) = ⤄(ex,√x)
⤄(tok::SourceLocToken, x::SourceRange) = SourceLocToken(tok.val, tok.loc ⤄ x)
⤄(tok::SourceLocToken, x::SourceLocToken) = SourceLocToken(tok.val, tok.loc ⤄ x.loc)
⤄(x::Union{ASTExprs, ASTVerbatim, Token},y::Union{ASTVerbatim,ASTExprs,Void,Token}) = x

function sortedcomplement(of::SourceRange, set)
    complement = SourceRange[]
    lastend = of.offset
    laststart = lastend
    for x in set
        x = normalize(x)
        @assert x.offset >= laststart
        if x.offset >= lastend
            push!(complement,SourceRange(lastend,x.offset-lastend,0))
        end
        laststart = max(laststart, x.offset)
        lastend = max(lastend, x.offset + x.length)
    end
    if lastend < of.offset + of.length
        push!(complement,SourceRange(lastend,of.offset+of.length-lastend-1,0))
    end
    complement
end

⨳(sym::Symbol) = Expr(sym)
⨳(sym::Symbol, args::LineNumberNode...) = Expr(sym::Symbol, args...)
⨳(sym::Symbol, args::Union{ASTVerbatim}...) = Expr(sym,args...)
function ⨳(sym::Symbol, args::Union{Token,ASTVerbatim,ASTExprs}...)
    ex = Expr(sym)
    for arg in args
        push!(ex.args, ¬arg)
    end
    ex
end
function ⨳(sym::Symbol, args::Union{ASTVerbatim,SourceLocToken,LineNumberNode,SourceExpr}...)
    accum = SourceRange()
    children = SourceNode[]
    ex = Expr(sym)
    for arg in args
        loc = √arg
        accum = accum ⤄ loc
        push!(children, SourceNode(loc))
        push!(ex.args, ¬arg)
    end
    SourceExpr(ex,SourceNode(normalize(accum),children))
end
⨳(sym::SourceLocToken,args...) = (SourceExpr(Expr(¬sym), SourceRange()) ⪥ args) ⤄ √sym
⨳(sym::Token,args...) = ⨳(¬sym,args...)

_push!(ex::Expr, t) = push!(ex.args, ¬t)
function _push!(ex::SourceExpr, t)
    push!((¬ex).args, ¬t); push!(ex.loc.children, SourceNode(√t))
end

function expr_append!(ex::SourceExpr, args::Union{Array,Tuple})
    !isempty(args) || return ex
    append!(ex.expr.args,[¬x for x in args])
    append!(ex.loc.children,[SourceNode(√x) for x in args])
    ex = SourceExpr(ex.expr,
        SourceNode(normalize(reduce(⤄,ex.loc.children)),ex.loc.children))
    ex
end
function expr_append!(ex::SourceExpr, new::SourceExpr)
    expr_append!(ex, collect(map(x->SourceExpr(x[1],SourceNode(x[2])),zip(new.expr.args,children(new.loc)))))
end
function expr_append!(ex::Expr, new::SourceExpr)
    expr_append!(SourceExpr(ex, SourceRange()), new)
end
expr_append!(ex::Expr, args::Array) = (for t in args; ex = ex ⤄ t; _push!(ex, t); end; ex)
expr_append!(ex::Expr, args::Tuple) = (for t in args; ex = ex ⤄ t; _push!(ex, t); end; ex)
expr_append!(ex::Expr, new::Expr) = expr_append!(ex, new.args)

const ⪥ = expr_append!

end


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

include("model_ApplicationGatewayBackendAddress.jl")
include("model_ApplicationGatewayBackendAddressPool.jl")
include("model_ApplicationGatewayBackendAddressPoolPropertiesFormat.jl")
include("model_ApplicationSecurityGroup.jl")
include("model_ApplicationSecurityGroupPropertiesFormat.jl")
include("model_BackendAddressPool.jl")
include("model_BackendAddressPoolPropertiesFormat.jl")
include("model_CloudError.jl")
include("model_CloudErrorBody.jl")
include("model_CustomDnsConfigPropertiesFormat.jl")
include("model_DdosSettings.jl")
include("model_Delegation.jl")
include("model_FlowLog.jl")
include("model_FlowLogFormatParameters.jl")
include("model_FlowLogPropertiesFormat.jl")
include("model_FrontendIPConfiguration.jl")
include("model_FrontendIPConfigurationPropertiesFormat.jl")
include("model_IPAllocationMethod.jl")
include("model_IPConfiguration.jl")
include("model_IPConfigurationProfile.jl")
include("model_IPConfigurationProfilePropertiesFormat.jl")
include("model_IPConfigurationPropertiesFormat.jl")
include("model_IPVersion.jl")
include("model_InboundNatRule.jl")
include("model_InboundNatRulePropertiesFormat.jl")
include("model_IpTag.jl")
include("model_LoadBalancerBackendAddress.jl")
include("model_LoadBalancerBackendAddressPropertiesFormat.jl")
include("model_NetworkInterface.jl")
include("model_NetworkInterfaceDnsSettings.jl")
include("model_NetworkInterfaceIPConfiguration.jl")
include("model_NetworkInterfaceIPConfigurationPrivateLinkConnectionProperties.jl")
include("model_NetworkInterfaceIPConfigurationPropertiesFormat.jl")
include("model_NetworkInterfacePropertiesFormat.jl")
include("model_NetworkInterfaceTapConfiguration.jl")
include("model_NetworkInterfaceTapConfigurationPropertiesFormat.jl")
include("model_NetworkSecurityGroup.jl")
include("model_NetworkSecurityGroupPropertiesFormat.jl")
include("model_PrivateEndpoint.jl")
include("model_PrivateEndpointProperties.jl")
include("model_PrivateLinkServiceConnection.jl")
include("model_PrivateLinkServiceConnectionProperties.jl")
include("model_PrivateLinkServiceConnectionState.jl")
include("model_ProvisioningState.jl")
include("model_PublicIPAddress.jl")
include("model_PublicIPAddressDnsSettings.jl")
include("model_PublicIPAddressPropertiesFormat.jl")
include("model_PublicIPAddressSku.jl")
include("model_Resource.jl")
include("model_ResourceNavigationLink.jl")
include("model_ResourceNavigationLinkFormat.jl")
include("model_RetentionPolicyParameters.jl")
include("model_Route.jl")
include("model_Route2.jl")
include("model_RouteListResult.jl")
include("model_RouteNextHopType.jl")
include("model_RouteNextHopType2.jl")
include("model_RoutePropertiesFormat.jl")
include("model_RoutePropertiesFormat2.jl")
include("model_RouteTable.jl")
include("model_RouteTable2.jl")
include("model_RouteTableListResult.jl")
include("model_RouteTablePropertiesFormat.jl")
include("model_RouteTablePropertiesFormat2.jl")
include("model_SecurityRule.jl")
include("model_SecurityRuleAccess.jl")
include("model_SecurityRuleDirection.jl")
include("model_SecurityRulePropertiesFormat.jl")
include("model_ServiceAssociationLink.jl")
include("model_ServiceAssociationLinkPropertiesFormat.jl")
include("model_ServiceDelegationPropertiesFormat.jl")
include("model_ServiceEndpointPolicy.jl")
include("model_ServiceEndpointPolicyDefinition.jl")
include("model_ServiceEndpointPolicyDefinitionPropertiesFormat.jl")
include("model_ServiceEndpointPolicyPropertiesFormat.jl")
include("model_ServiceEndpointPropertiesFormat.jl")
include("model_SubResource.jl")
include("model_Subnet.jl")
include("model_SubnetPropertiesFormat.jl")
include("model_TagsObject.jl")
include("model_TrafficAnalyticsConfigurationProperties.jl")
include("model_TrafficAnalyticsProperties.jl")
include("model_TransportProtocol.jl")
include("model_VirtualNetworkTap.jl")
include("model_VirtualNetworkTapPropertiesFormat.jl")


text:
using FFTW, LinearAlgebra

export MicMac, solve

mutable struct MicMac

    data      :: DataSet
    ntau      :: Int64
    ktau      :: Vector{Float64}
    matr      :: Array{ComplexF64,2}
    conjmatr  :: Array{ComplexF64,2}
    A1        :: Vector{Float64}
    A2        :: Vector{Float64}
    sigma     :: Int64
    llambda   :: Int64
    epsilon   :: Float64
    u         :: Array{ComplexF64,1}
    v         :: Array{ComplexF64,1}
    ut        :: Array{ComplexF64,2}
    vt        :: Array{ComplexF64,2}
    z         :: Array{ComplexF64,2}

    function MicMac( data :: DataSet, ntau :: Int64 )

        nx        = data.nx
        T         = data.T
        kx        = data.kx

        epsilon   = data.epsilon

        llambda   = data.llambda
        sigma     = data.sigma

        tau       = zeros(Float64, ntau)
        tau      .= T * collect(0:ntau-1) / ntau
        ktau      = similar(tau)
        ktau     .= 2π / T * vcat(0:ntau÷2-1,-ntau÷2:-1)
        ktau[1]   = 1.0

        matr      = zeros(ComplexF64,(1,ntau))
        conjmatr  = zeros(ComplexF64,(1,ntau))
        matr     .= transpose(exp.( 1im * tau))
        conjmatr .= transpose(exp.(-1im * tau))

        A1 = zeros(Float64, nx)
        A2 = zeros(Float64, nx)

        if epsilon > 0
            A1 .= (sqrt.(1 .+ epsilon * kx .^2) .- 1) / epsilon
            A2 .= (1 .+ epsilon * kx .^2) .^ (-1/2)
        else
            A1 .= 0.5 * kx .^ 2
            A2 .= 1.0
        end

        u  = zeros(ComplexF64,nx)
        v  = zeros(ComplexF64,nx)
        ut = zeros(ComplexF64,(nx, ntau))
        vt = zeros(ComplexF64,(nx, ntau))
        z  = zeros(ComplexF64,(nx, ntau))

        new( data, ntau, ktau, matr, conjmatr, A1 , A2,
             sigma, llambda, epsilon, u, v, ut, vt, z)

    end

end

include("ftau.jl")
include("C1.jl")
include("dtftau.jl")
include("duftau.jl")
include("init_2.jl")
include("champs_2.jl")

function reconstr(u, t, T, ntau)

    v    = fft(u,1)
    w    = vcat(0:ntau÷2-1, -ntau÷2:-1)
    v  .*= exp.(1im * 2π / T * w * t)

    vec(sum(v, dims=1) / ntau)

end

function solve(self, dt)

    Tfinal = self.data.Tfinal

    T  = self.data.T
    kx = self.data.kx

    epsilon = self.data.epsilon

    nx   = self.data.nx
    ntau = self.ntau

    t = 0.0
    iter = 0

    fft_u = copy(self.data.u)
    fft_v = copy(self.data.v)

    fft!(fft_u,1)
    fft!(fft_v,1)

    fft_ubar = copy(fft_u)
    fft_vbar = copy(fft_v)

    fft_ug = copy(fft_u)
    fft_vg = copy(fft_v)

    ichampgu = zeros(ComplexF64,(ntau,nx))
    ichampgv = zeros(ComplexF64,(ntau,nx))

    init_2!(ichampgu, ichampgv, self, t, fft_ubar, fft_vbar, fft_ug, fft_vg )

    champubaru = similar(fft_ubar)
    champubarv = similar(fft_vbar)
    champmoyu  = similar(fft_ug)
    champmoyv  = similar(fft_vg)

    while t < Tfinal

        iter = iter + 1
        dt   = min(Tfinal-t, dt)
        hdt  = dt / 2

        champubaru .= fft_ubar
        champubarv .= fft_vbar
        champmoyu  .= fft_ug
        champmoyv  .= fft_vg

        champs_2!(ichampgu, ichampgv,
                  self, t,
                  champubaru, champubarv,
                  champmoyu, champmoyv )

        champubaru .*= hdt
        champubaru .+= fft_ubar
        champubarv .*= hdt
        champubarv .+= fft_vbar

        champmoyu  .*= hdt
        champmoyu  .+= fft_ug
        champmoyu  .+= epsilon .* reconstr(ichampgu, (t+hdt) / epsilon, T, ntau)
        champmoyu  .-= epsilon .* reconstr(ichampgu,  t      / epsilon, T, ntau)

        champmoyv  .*= hdt
        champmoyv  .+= fft_vg
        champmoyv  .+= epsilon .* reconstr(ichampgv, (t+hdt) / epsilon, T, ntau)
        champmoyv  .-= epsilon .* reconstr(ichampgv,  t      / epsilon, T, ntau)

        champs_2!(ichampgu, ichampgv,
                  self, t + hdt,
                  champubaru, champubarv,
                  champmoyu, champmoyv )

        fft_ubar .+= dt .* champubaru
        fft_vbar .+= dt .* champubarv

        fft_ug .+= epsilon .* reconstr(ichampgu, (t + dt) / epsilon, T, ntau)
        fft_ug .-= epsilon .* reconstr(ichampgu,  t       / epsilon, T, ntau)
        fft_ug .+= dt .* champmoyu

        fft_vg .+= epsilon .* reconstr(ichampgv, (t + dt) / epsilon, T, ntau)
        fft_vg .-= epsilon .* reconstr(ichampgv,  t       / epsilon, T, ntau)
        fft_vg .+= dt .* champmoyv

        t = t + dt

        C1!(ichampgu, ichampgv, self, t, fft_ubar, fft_vbar)

        fft_u .= reconstr(ichampgu, t / epsilon, T, ntau)
        fft_v .= reconstr(ichampgv, t / epsilon, T, ntau)

        fft_u .+= fft_ug
        fft_v .+= fft_vg

    end

    fft_u .*= exp.(1im .* sqrt.(1 .+ epsilon * kx .^ 2) .* t / epsilon)
    fft_v .*= exp.(1im .* sqrt.(1 .+ epsilon * kx .^ 2) .* t / epsilon)

    ifft!(fft_u,1)
    ifft!(fft_v,1)

    fft_u, fft_v

end


text:
using DataFrames

mdata = readtable("original_data/mdata.csv")
qdata = readtable("original_data/qdata.csv")

# y/y log change in GDP, CPI, and IPindex
GDP = (log(qdata[:GDPlev][5:end]) - log(qdata[:GDPlev][1:(end - 4)]))*100
INFL = (log(mdata[:CPI][13:end]) - log(mdata[:CPI][1:(end - 12)]))*100
IP = (log(mdata[:INDPRO][13:end]) - log(mdata[:INDPRO][1:(end - 12)]))*100

mdata_processed = DataFrame(date = mdata[:date][13:end], IP = IP, INFL = INFL, UNEMP = mdata[:UNEMP][13:end])
qdata_processed = DataFrame(date = qdata[:date][5:end], GDP = GDP)

writetable("importable_data/mdata_processed.csv", mdata_processed)
writetable("importable_data/qdata_processed.csv", qdata_processed)

max_stars_repo_name: lnsongxf/MF-VAR

max_issues_repo_name: lnsongxf/MF-VAR

max_forks_repo_name: lnsongxf/MF-VAR


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using FEMBase, Materials, FEMMaterials, Test

@testset "Test FEMMaterials.jl" begin
    @testset "test simulator" begin
        include("test_simulator.jl")
    end
    @testset "test mechanical material solver" begin
        include("test_mecamatso.jl")
    end
    @testset "test material expression" begin
        include("test_mecamatso_expr.jl")
    end
    @testset "test ideal plastic material model with shear strain" begin
        include("test_idealplastic_shear.jl")
    end
    @testset "test chaboche material model" begin
        include("test_chaboche.jl")
    end
    @testset "test chaboche material model with shear strain" begin
        include("test_chaboche_shear.jl")
    end
    """
    @testset "test viscoplastic" begin
        include("test_viscoplastic.jl")
    end
    """
end


text:
function rpn(s)
    stack = Any[]
    for op in map(eval, map(parse, split(s)))
        if isa(op, Function)
            arg2 = pop!(stack)
            arg1 = pop!(stack)
            push!(stack, op(arg1, arg2))
        else
            push!(stack, op)
        end
        println("$op: ", join(stack, ", "))
    end
    length(stack) != 1 && error("invalid RPN expression $s")
    return stack[1]
end
rpn("3 4 2 * 1 5 - 2 3 ^ ^ / +")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using InverseIsing
using LinearAlgebra:Symmetric
using Test
using Dates
using PyPlot

y = [] # processing time
x = [2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 50, 70, 90, 100] # num_unit
for n_unit in x
    samples = ones(1000, n_unit)
    model = GBM(n_unit)
    print("\nn_unit:", n_unit)
    start = time_ns()
    fit(model, samples, delta=0.0)
    push!(y, Float64(time_ns() - start))
    d = decode(infer(model))
    flg = collect(values(d)) == ones(binomial(n_unit, 2))
    print(" -> ", flg)
end
plot(x, y, ".-")
xlabel("num_unit")
ylabel("time")

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
"""
**data_conversion.jl** contains functions to convert imported data from e.g. CSV or netCDF files to the GeoData format
that is used internally to store all data related to a certain data set. For importing files, use standard methods, such
as CSV import using *readdlm* (see the [DelimitedFiles.jl](https://docs.julialang.org/en/v1/stdlib/DelimitedFiles/) package)
or netCDF import (see the [NetCDF.jl](https://github.com/JuliaGeo/NetCDF.jl) package).
"""

"""
Using *readdlm* on CSV files will provide two output arguments, one containing the header as a Array{AbstractString, 2} and
the other one containing the data as Array{Float64, 2}.
"""

############ CONVERT DLM DATA TO GEO DATA
function DLM2Geo(hdr::Array{AbstractString, 2},data::Array{Float64, 2},DepthCon::AbstractString)

    # initialize array of structures to store the data
    # while doing so, separate the unit from the variable name
    ndata   = size(data,1) # number of entries
    nfields = size(data,2) # number of fields

    # declare some variables as local, otherwise they are not known outside of the following for loop
    local LonData
    local LatData
    local DepthData
    local vals_range
    vals_range = zeros(Int64,nfields-3)
    ivals = 1;
    # get the fields for lon/lat/depth
    for ifield = 1:nfields
        if occursin("lon",hdr[ifield])
            lon_ind = ifield;
            varname = GetVariableName(hdr[ifield])# get variable name
            varunit = GetVariableUnit(hdr[ifield])# get variable unit
            LonData = data[1:end,ifield]
        elseif occursin("lat",hdr[ifield])
            lat_ind = ifield;
            varname = GetVariableName(hdr[ifield])# get variable name
            varunit = GetVariableUnit(hdr[ifield])# get variable unit
            LatData = data[1:end,ifield]
        elseif occursin("depth",hdr[ifield])
            # ISSUE: WE DEFINE DEPTH AS NEGATIVE, BUT HOW DO WE SET THAT?
            # WE COULD ADD A FLAG THAT INDICATES THE DEPTH CONVENTION AND
            # TREAT IT ACCORDINGLY
            depth_ind = ifield;
            varname = GetVariableName(hdr[ifield])# get variable name
            varunit = GetVariableUnit(hdr[ifield])# get variable unit

            # take care of positive or negative convection for depth (here we use negative)
            if cmp(DepthCon,"positive")==0 # if depth is given as positive values, convert to negative
                DepthData = -1*data[1:end,ifield]
            elseif cmp(DepthCon,"negative")==0
                DepthData = data[1:end,ifield]
            else # default behaviour assumes that dpeth is negative
                DepthData = data[1:end,ifield]
            end

            # if depth is given in m, convert to km
            if cmp(varunit,"m")==0
                DepthData = DepthData./1e3;
            end

        else
            vals_range[ivals] = ifield
            ivals = ivals+1
        end
    end

    # check if the CSV data is actually given on a regular lon/lat/depth grid and rearrange to create a 3D matrix
    lat_unique = unique(LatData)
    lon_unique = unique(LonData)
    depth_unique = unique(DepthData)

    if size(lat_unique,1)*size(lon_unique,1)*size(depth_unique,1) == size(data,1)
            println("Data has to be converted to 3D matrix")
    end

    # create named tuple for additional data
    tmp_hdr  = hdr[vals_range];
    tmp_data = data[1:end,vals_range];

    nhdr = size(tmp_hdr,1)
    tmp_vec = Vector{Vector{Float64}}(undef, nhdr) # this is used for later tuple creation, I haven't found a better way around

    for ihdr = 1:nhdr

        # take care of the header strings
        varname = GetVariableName(tmp_hdr[ihdr])# get variable name
        varunit = GetVariableUnit(tmp_hdr[ihdr])# get variable unit
        if cmp(varunit,"%")==0
            tmp_hdr[ihdr] = string(varname,"_percentage")
        else
            tmp_hdr[ihdr] = string(varname,"_",varunit)
        end

        # take care of the matrix columns
        tmp_vec[ihdr] = tmp_data[1:end,ihdr];
    end

    hdr_tpl  = Tuple(Symbol(x) for x in tmp_hdr) # convert header to tuple
    data_tpl = Tuple.(tmp_vec for i in size(tmp_vec,1)) # convert data to tuple
    tmp = NamedTuple{hdr_tpl}(data_tpl)

    println(typeof(tmp))

    # initialize data structure
    importdata = GeoData(LonData,LatData,DepthData,tmp)

    # assign data to output
    return importdata

end

########### REARRANGE DATA TO OBTAIN A 3D MATIX IF NECESSARY ##########
function RearrangeToMatrix()

end
########### CONVERT NETCDF DATA TO GEO DATA ########
"""
Converting netCDF data to GeoData is fairly easy. One can read in Data from a netCDF file using ncread("filename","variable")
(the contents of the netcdf file can be queried beforehand using ncinfo("filename")). Via *ncread*, one then reads in all the
desired variables. NetCDF2Geo then takes care of converting this data to GeoData.
"""
function NetCDF2Geo()

end

#############################################################################################
############################ SOME USEFUL FUNCTIONS FOR STRING MANIPULATION ##################
function GetVariableName(inputstring::SubString{String})
    # convert input to normal String
    inputstring = String(inputstring)
    # assume that if the string contains a unit, it is given in brackets (),[],{}
    indfirst = nothing
    iloop     = 1
    str2find = ["(","[","{"]
    # find first occurence of one of the brackets
    while isnothing(indfirst)
        indfirst = findfirst(str2find[iloop],inputstring)
        iloop = iloop + 1
        if iloop>length(str2find)
            break
        end
    end

    # either return the whole inputstring or only the part before the unit
    if isnothing(indfirst)
        return inputstring
    else
        indfirst = indfirst[1]-1
        return inputstring[1:indfirst]
    end
end

function GetVariableUnit(inputstring::SubString{String})
    # convert input to normal String
    inputstring = String(inputstring)
    # assume that if the string contains a unit, it is given in brackets (),[],{}
    indfirst = nothing
    iloop     = 1;
    firststr2find = ["(","[","{"]
    laststr2find = [")","]","}"]
    while isnothing(indfirst)
        indfirst = findfirst(firststr2find[iloop],inputstring)
        iloop = iloop + 1
        if iloop>length(firststr2find)
            break
        end
    end

    if isnothing(indfirst)
        return "none"
    else
        indlast = findfirst(laststr2find[iloop-1],inputstring)
        indfirst = indfirst[1]+1
        indlast = indlast[1]-1
        return inputstring[indfirst:indlast]
    end

end


text:
rand_glorot(_s) = (CuArrays.rand(_s...) .- 0.5f0) .* sqrt(24f0 / sum(_s))

apply!(opt::Optimizer, Δ::GlobalArrayType, df::Function) = opt.f(Δ, df, opt.params)

forward!(node::NetworkNode, x::GlobalArrayType) = begin
    y, d = node.forward(x, node.params)
    for (k, v) in d
        node.derivs[k] = v
    end
    return y
end
backward!(node::NetworkNode, Δ::GlobalArrayType, optim::Optimizer) = begin
    node.apply_derivs!(node, Δ, optim)

    @assert haskey(node.derivs, "dydx") "layer has no 'dydx' key"

    return node.derivs["dydx"](Δ)
end

forward!(nodes::Array{NetworkNode}, x::GlobalArrayType) = begin
    for n in nodes
        x = forward!(n, x)
    end
    return x
end
backward!(nodes::Array{NetworkNode}, Δ::GlobalArrayType, opt::Optimizer) = begin
    for n in nodes
        Δ = backward!(n, Δ, opt)
    end
    return Δ
end

# Convolution
function flat_kernel_indecies(input_size, K, D)
    height = input_size[1]
    width = prod(input_size[1:2])
    chan = prod(input_size[1:3])

    h_i = [ki * D[1] for ki in 0:K[1]-1]
    w_i = [height * ki * D[2] for ki in 0:K[2]-1]
    c_i = [width * ki * D[3] for ki in 0:K[3]-1]
    b_i = [chan * ki * D[4] for ki in 0:K[4]-1]

    h_w = [h_i .+ w for w in w_i] |> Iterators.flatten |> collect
    h_w_c = [h_w .+ c for c in c_i] |> Iterators.flatten |> collect
    h_w_c_b = [h_w_c .+ b for b in b_i] |> Iterators.flatten |> collect

    return h_w_c_b
end
function pad_same(x, p)
    for _ in 1:p[1]
        s = size(x)
        p_d1 = reshape(x[end, :, :, :], (1, s[2], s[3], s[4]))
        x = cat(x, p_d1, dims=1)
    end

    for _ in 1:p[2]
        s = size(x)
        p_d2 = reshape(x[:, end, :, :], (s[1], 1, s[3], s[4]))
        x = cat(x, p_d2, dims=2)
    end

    for _ in 1:p[3]
        s = size(x)
        p_d3 = reshape(x[:, :, end, :], (s[1], s[2], 1, s[4]))
        x = cat(x, p_d3, dims=3)
    end

    for _ in 1:p[4]
        s = size(x)
        p_d4 = reshape(x[:, :, :, end], (s[1], s[2], s[3], 1))
        x = cat(x, p_d4, dims=4)
    end

    return x
end
function get_valid_indecies(input_size, output_size, S) # get indecies for each valid position in image based on the kernels shape

    idim1 = [ai * S[1] + 1 for ai in 0:output_size[1]-1]
    idim2 = [ai * input_size[1] * S[2] + 1 for ai in 0:output_size[2]-1]
    idim3 = [ai * input_size[1] * input_size[2] * S[3] + 1 for ai in 0:output_size[3]-1]
    idim4 = [ai * input_size[1] * input_size[2] * input_size[3] * S[4] + 1 for ai in 0:output_size[4]-1]

    I = idim1

    if length(idim2) > 1
        append!(I, [id2 .+ I .- 1 for id2 in idim2[2:end]] |> Iterators.flatten |> collect)
    end
    if length(idim3) > 1
        append!(I, [id3 .+ I .- 1 for id3 in idim3[2:end]] |> Iterators.flatten |> collect)
    end
    if length(idim4) > 1
        append!(I, [id4 .+ I .- 1 for id3 in idim4[2:end]] |> Iterators.flatten |> collect)
    end
    I
end


text:
export upperbound, lowerbound

function upperbound(a, b::Curve{D}, β::Interval) where {D}
    evalpt = b(mid(β))
    return minimum_distance(a, evalpt, @SVector(ones(D)))
end

function upperbound(a::Curve, b::Curve, α::Interval, β::Interval)
    norm(a(mid(α)) - b(mid(β)))
end

function lowerbound(a, b::Curve{D}, β::Interval) where {D}
    approx = cvxhull(b, β)
    return minimum_distance(a, approx, @SVector(ones(D)))
end

function lowerbound(a::Curve{D}, b::Curve{D}, α::Interval, β::Interval) where {D}
    bapprox = cvxhull(b, β)
    aapprox = cvxhull(a, α)
    return minimum_distance(aapprox, bapprox, @SVector(ones(D)))
end

function getbounds(a, b, θ...)
    lb = lowerbound(a, b, θ...)
    ub = upperbound(a, b, θ...)
    return lb, ub
end

import ConvexBodyProximityQueries: support
support(pt::SVector{D, T}, dir::SVector{D, T}) where {D, N, T} = pt
function support(vertices::SVector{N, SVector{D, T}}, dir::SVector{D, T}) where {D, N, T}
    @inbounds vertices[argmax(Ref(dir').*vertices)]
end

function cvxhull(b::Curve{2, T}, β::Interval) where {T}
    fA = b(β.lo)
    fB = b(β.hi)
    center = (fA + fB)/2.0
    fdist = norm(fA - fB)
    rotate = SMatrix{2,2}(1.0I)
    if fdist > 0.0
        dir = (fB-fA)/fdist
        rotate = SMatrix{2,2}(dir[1],dir[2],-dir[2],dir[1])
    end
    a = arclength(b, β)/2.0
    if (a - fdist/2.0) > eps(T)
        b = sqrt(a^2 - (fdist/2)^2)
        octagon = SVector{8}(
                   SVector{2}(a*(√2 - 1), b),
                   SVector{2}(a, b*(√2 - 1)),
                   SVector{2}(a, -b*(√2 - 1)),
                   SVector{2}(a*(√2 - 1), -b),
                   SVector{2}(-a*(√2 - 1), -b),
                   SVector{2}(-a, -b*(√2 - 1)),
                   SVector{2}(-a, b*(√2 - 1)),
                   SVector{2}(-a*(√2 - 1), b),
                  )
        vertices = Ref(rotate).*octagon .+ Ref(center)
    else
        vertices = SVector{2}(fA, fB)
    end
end

function cvxhull(b::Curve{3, T}, β::Interval) where {T}
    fA = b(β.lo)
    fB = b(β.hi)
    center = (fA + fB)/2.0
    fdist = norm(fA - fB)
    rotate = SMatrix{3,3}(1.0I)
    if fdist > 0.0
        dir = (fB-fA)/fdist
        aux = SVector{3}(1.0, 0.0, 0.0)
        aux = abs(dir'*aux) < 0.99 ? aux : SVector{3}(0.0, 1.0, 0.0)
        auxperp = normalize(aux × dir)
        rotate = hcat(dir, aux, auxperp)
    end
    a = arclength(b, β)/2.0
    if (a - fdist/2.0) > eps(T)
        b = sqrt(a^2 - (fdist/2)^2)
        tetradecahedron = SVector{16}(
                   SVector{3}( a,  b*(√2 - 1),  b*(√2 - 1)),
                   SVector{3}( a, -b*(√2 - 1),  b*(√2 - 1)),
                   SVector{3}( a, -b*(√2 - 1), -b*(√2 - 1)),
                   SVector{3}( a,  b*(√2 - 1), -b*(√2 - 1)),
                   SVector{3}( a*(√2 - 1),  b,  b),
                   SVector{3}( a*(√2 - 1), -b,  b),
                   SVector{3}( a*(√2 - 1), -b, -b),
                   SVector{3}( a*(√2 - 1),  b, -b),
                   SVector{3}(-a*(√2 - 1),  b,  b),
                   SVector{3}(-a*(√2 - 1), -b,  b),
                   SVector{3}(-a*(√2 - 1), -b, -b),
                   SVector{3}(-a*(√2 - 1),  b, -b),
                   SVector{3}(-a,  b*(√2 - 1),  b*(√2 - 1)),
                   SVector{3}(-a, -b*(√2 - 1),  b*(√2 - 1)),
                   SVector{3}(-a, -b*(√2 - 1), -b*(√2 - 1)),
                   SVector{3}(-a,  b*(√2 - 1), -b*(√2 - 1)),
                  )
        vertices = Ref(rotate).*tetradecahedron .+ Ref(center)
    else
        vertices = SVector{2}(fA, fB)
    end
end


text:
@testset "Tensor operations" begin
@testset "right propagation of A and B" begin
    nl1 = 2
    nl2 = 3
    na = 4
    nb = 5
    ns = 6

    # Float64 propagation.
    L = reshape(collect(1. : nl1*nl2), (nl1, nl2))
    A = reshape(collect(1. : nl1*ns*na), (nl1, ns, na))
    B = reshape(collect(3. : 2+nl2*ns*nb), (nl2, ns, nb))
    # Do manual contraction first.
    new_L = zeros(ComplexF64, (na, nb))
    for a=1:na, b=1:nb
        for l1=1:nl1, l2=1:nl2, s=1:ns
            new_L[a, b] += L[l1, l2]*A[l1, s, a]*B[l2, s, b]
        end
    end
    @test MPStates.prop_right2(L, A, B) ≈ new_L

    # ComplexF64 propagation.
    L = reshape(complex.(collect(1. : nl1*nl2)), (nl1, nl2))
    A = reshape(complex.(collect(1. : nl1*ns*na), collect(-2. : -3+nl1*ns*na)),
                (nl1, ns, na))
    B = reshape(complex.(collect(3. : 2+nl2*ns*nb), collect(-1. : -2+nl2*ns*nb)),
                (nl2, ns, nb))
    # Do manual contraction first.
    new_L = zeros(ComplexF64, (na, nb))
    for a=1:na, b=1:nb
        for l1=1:nl1, l2=1:nl2, s=1:ns
            new_L[a, b] += L[l1, l2]*A[l1, s, a]*conj(B[l2, s, b])
        end
    end
    @test MPStates.prop_right2(L, A, B) ≈ new_L
end

@testset "right propagation of A, M, and B" begin
    nl1 = 2
    nl2 = 3
    nl3 = 4
    na = 4
    nb = 5
    nc = 3
    ns1 = 6
    ns2 = 2
    L = reshape(complex.(collect(1. : nl1*nl2*nl3)), (nl1, nl2, nl3))
    A = reshape(collect(1. : nl1*ns1*na) + 1im*collect(-2. : -3+nl1*ns1*na),
                (nl1, ns1, na))
    M = reshape(collect(1. : nl2*ns1*ns2*nb) + 1im*collect(0. : -1+nl2*ns1*ns2*nb),
                (nl2, ns1, ns2, nb))
    B = reshape(collect(3. : 2+nl3*ns2*nc) + 1im*collect(-1. : -2+nl3*ns2*nc),
                (nl3, ns2, nc))
    # Do manual contraction first.
    new_L = zeros(ComplexF64, (na, nb, nc))
    for a=1:na, b=1:nb, c=1:nc
        for l1=1:nl1, l2=1:nl2, l3=1:nl3, s1=1:ns1, s2=1:ns2
            new_L[a, b, c] += L[l1, l2, l3]*A[l1, s1, a]*M[l2, s1, s2, b]*conj(B[l3, s2, c])
        end
    end
    @test MPStates.prop_right3(L, A, M, B) ≈ new_L
end

@testset "right propagation of M, and B for DMRG3S" begin
    nl1 = 2
    nl2 = 3
    nl3 = 4
    nr2 = 5
    nr3 = 3
    ns1 = 6
    ns2 = 2
    L = reshape(complex.(collect(1. : nl1*nl2*nl3)), (nl1, nl2, nl3))
    M = reshape(complex.(collect(1. : nl2*ns1*ns2*nr2), collect(0. : -1+nl2*ns1*ns2*nr2)),
                (nl2, ns1, ns2, nr2))
    B = reshape(complex.(collect(3. : 2+nl3*ns2*nr3), collect(-1. : -2+nl3*ns2*nr3)),
                (nl3, ns2, nr3))
    # Do manual contraction first.
    P = zeros(ComplexF64, (nl1, ns1, nr2, nr3))
    for l1=1:nl1, s1=1:ns1, r2=1:nr2, r3=1:nr3
        for l2=1:nl2, l3=1:nl3, s2=1:ns2
            P[l1, s1, r2, r3] += L[l1, l2, l3]*M[l2, s1, s2, r2]*conj(B[l3, s2, r3])
        end
    end
    @test MPStates.prop_right_subexp(L, M, B) ≈ P
end

@testset "right propagation of A, M1, M2, and B" begin
    nl1 = 2
    nl2 = 3
    nl3 = 4
    nl4 = 4
    na = 4
    nb = 5
    nc = 3
    nd = 5
    ns1 = 6
    ns2 = 2
    ns3 = 4
    L = reshape(complex.(collect(1. : nl1*nl2*nl3*nl4)), (nl1, nl2, nl3, nl4))
    A = reshape(collect(1. : nl1*ns1*na) + 1im*collect(-2. : -3+nl1*ns1*na),
                (nl1, ns1, na))
    M1 = reshape(collect(1. : nl2*ns1*ns2*nb) + 1im*collect(0. : -1+nl2*ns1*ns2*nb),
                 (nl2, ns1, ns2, nb))
    M2 = reshape(collect(1. : nl3*ns2*ns3*nc) + 1im*collect(0. : -1+nl3*ns2*ns3*nc),
                 (nl3, ns2, ns3, nc))
    B = reshape(collect(3. : 2+nl4*ns3*nd) + 1im*collect(-1. : -2+nl4*ns3*nd),
                (nl4, ns3, nd))
    # Do manual contraction first.
    new_L = zeros(ComplexF64, (na, nb, nc, nd))
    for a=1:na, b=1:nb, c=1:nc, d=1:nd
        for l1=1:nl1, l2=1:nl2, l3=1:nl3, l4=1:nl4, s1=1:ns1, s2=1:ns2, s3=1:ns3
            new_L[a, b, c, d] += (L[l1, l2, l3, l4]*A[l1, s1, a]*M1[l2, s1, s2, b]
                                  *M2[l3, s2, s3, c]*conj(B[l4, s3, d]))
        end
    end
    @test MPStates.prop_right4(L, A, M1, M2, B) ≈ new_L
end

@testset "left propagation of A and B" begin
    nr1 = 2
    nr2 = 3
    na = 4
    nb = 5
    ns = 6

    # Float64 propagation.
    R = reshape(collect(1. : nr1*nr2), (nr1, nr2))
    A = reshape(collect(1. : nr1*ns*na), (na, ns, nr1))
    B = reshape(collect(3. : 2+nr2*ns*nb), (nb, ns, nr2))
    # Do manual contraction first.
    new_R = zeros(ComplexF64, (na, nb))
    for a=1:na, b=1:nb
        for r1=1:nr1, r2=1:nr2, s=1:ns
            new_R[a, b] += R[r1, r2]*A[a, s, r1]*B[b, s, r2]
        end
    end
    @test MPStates.prop_left2(A, B, R) ≈ new_R

    # ComplexF64 propagation.
    R = reshape(complex.(collect(1. : nr1*nr2)), (nr1, nr2))
    A = reshape(complex.(collect(1. : nr1*ns*na), collect(-2. : -3+nr1*ns*na)),
                (na, ns, nr1))
    B = reshape(complex.(collect(3. : 2+nr2*ns*nb), collect(-1. : -2+nr2*ns*nb)),
                (nb, ns, nr2))
    # Do manual contraction first.
    new_R = zeros(ComplexF64, (na, nb))
    for a=1:na, b=1:nb
        for r1=1:nr1, r2=1:nr2, s=1:ns
            new_R[a, b] += R[r1, r2]*A[a, s, r1]*conj(B[b, s, r2])
        end
    end
    @test MPStates.prop_left2(A, B, R) ≈ new_R
end

@testset "left propagation of A, M, and B" begin
    nr1 = 2
    nr2 = 3
    nr3 = 4
    na = 4
    nb = 5
    nc = 3
    ns1 = 6
    ns2 = 2
    R = reshape(complex.(collect(1. : nr1*nr2*nr3)), (nr1, nr2, nr3))
    A = reshape(collect(1. : nr1*ns1*na) + 1im*collect(-2. : -3+nr1*ns1*na),
                (na, ns1, nr1))
    M = reshape(collect(1. : nr2*ns1*ns2*nb) + 1im*collect(0. : -1+nr2*ns1*ns2*nb),
                (nb, ns1, ns2, nr2))
    B = reshape(collect(3. : 2+nr3*ns2*nc) + 1im*collect(-1. : -2+nr3*ns2*nc),
                (nc, ns2, nr3))
    # Do manual contraction first.
    new_R = zeros(ComplexF64, (na, nb, nc))
    for a=1:na, b=1:nb, c=1:nc
        for r1=1:nr1, r2=1:nr2, r3=1:nr3, s1=1:ns1, s2=1:ns2
            new_R[a, b, c] += R[r1, r2, r3]*A[a, s1, r1]*M[b, s1, s2, r2]*conj(B[c, s2, r3])
        end
    end
    @test MPStates.prop_left3(A, M, B, R) ≈ new_R
end

@testset "left propagation of M, and B for DMRG3S" begin
    nr1 = 2
    nr2 = 3
    nr3 = 4
    nl2 = 5
    nl3 = 3
    ns1 = 6
    ns2 = 2
    R = reshape(complex.(collect(1. : nr1*nr2*nr3)), (nr1, nr2, nr3))
    M = reshape(complex.(collect(1. : nr2*ns1*ns2*nl2), (0. : -1+nr2*ns1*ns2*nl2)),
                (nl2, ns1, ns2, nr2))
    B = reshape(complex.(collect(3. : 2+nr3*ns2*nl3), (-1. : -2+nr3*ns2*nl3)),
                (nl3, ns2, nr3))
    # Do manual contraction first.
    P = zeros(ComplexF64, (ns1, nr1, nl2, nl3))
    for r1=1:nr1, s1=1:ns1, l2=1:nl2, l3=1:nl3
        for r2=1:nr2, r3=1:nr3, s2=1:ns2
            P[s1, r1, l2, l3] += R[r1, r2, r3]*M[l2, s1, s2, r2]*conj(B[l3, s2, r3])
        end
    end
    @test MPStates.prop_left_subexp(M, B, R) ≈ P
end

@testset "left propagation of A, M1, M2, and B" begin
    nr1 = 2
    nr2 = 3
    nr3 = 4
    nr4 = 4
    na = 4
    nb = 5
    nc = 3
    nd = 5
    ns1 = 6
    ns2 = 2
    ns3 = 4
    R = reshape(complex.(collect(1. : nr1*nr2*nr3*nr4)), (nr1, nr2, nr3, nr4))
    A = reshape(collect(1. : nr1*ns1*na) + 1im*collect(-2. : -3+nr1*ns1*na),
                (na, ns1, nr1))
    M1 = reshape(collect(1. : nr2*ns1*ns2*nb) + 1im*collect(0. : -1+nr2*ns1*ns2*nb),
                 (nb, ns1, ns2, nr2))
    M2 = reshape(collect(1. : nr3*ns2*ns3*nc) + 1im*collect(0. : -1+nr3*ns2*ns3*nc),
                 (nc, ns2, ns3, nr3))
    B = reshape(collect(3. : 2+nr4*ns3*nd) + 1im*collect(-1. : -2+nr4*ns3*nd),
                (nd, ns3, nr4))
    # Do manual contraction first.
    new_R = zeros(ComplexF64, (na, nb, nc, nd))
    for a=1:na, b=1:nb, c=1:nc, d=1:nd
        for r1=1:nr1, r2=1:nr2, r3=1:nr3, r4=1:nr4, s1=1:ns1, s2=1:ns2, s3=1:ns3
            new_R[a, b, c, d] += (R[r1, r2, r3, r4]*A[a, s1, r1]*M1[b, s1, s2, r2]
                                  *M2[c, s2, s3, r3]*conj(B[d, s3, r4]))
        end
    end
    @test MPStates.prop_left4(A, M1, M2, B, R) ≈ new_R
end
end # @testset "Tensor operations"


text:
using PyCall
datetime = pyimport("datetime")

import JSON as js
function typing_test()::Int64
    a::Vector{Int64} = [1, 2, 3]
    return a[2]
end

function calendar_test()
    return isoformat(utcnow(dt.datetime))
end

function date_to_json(objDate::dt.datetime)::js
    return Dict(
        "__type__" => "datetime",
        "year" => objDate.year,
        "month" => objDate.month,
        "day" => objDate.day,
        "hour" => objDate.hour,
        "minute" => objDate.minute,
        "second" => objDate.second,
        "microsecond" => objDate.microsecond,
        "tz" => (tzname(objDate.tzinfo, objDate), total_seconds(utcoffset(objDate))),
    )
end

function calendar_json_test()::js
    now = now(dt.datetime, dt.timezone.utc)
    return dumps(js, now, date_to_json)
end

if abspath(PROGRAM_FILE) == @__FILE__
    b = typing_test()
    @assert(b == 2)
    @assert(match(re, "\\d{4}-\\d{2}-\\d{2}T\\d{2}:\\d{2}:\\d{2}.\\d{6}", calendar_test()))
    res = calendar_json_test()
    println(res)
    println("OK")
end


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

mutable struct WinRMListener <: SwaggerModel
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    certificateUrl::Any # spec type: Union{ Nothing, String } # spec name: certificateUrl

    function WinRMListener(;protocol=nothing, certificateUrl=nothing)
        o = new()
        validate_property(WinRMListener, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(WinRMListener, Symbol("certificateUrl"), certificateUrl)
        setfield!(o, Symbol("certificateUrl"), certificateUrl)
        o
    end
end # type WinRMListener

const _property_map_WinRMListener = Dict{Symbol,Symbol}(Symbol("protocol")=>Symbol("protocol"), Symbol("certificateUrl")=>Symbol("certificateUrl"))
const _property_types_WinRMListener = Dict{Symbol,String}(Symbol("protocol")=>"String", Symbol("certificateUrl")=>"String")
Base.propertynames(::Type{ WinRMListener }) = collect(keys(_property_map_WinRMListener))
Swagger.property_type(::Type{ WinRMListener }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_WinRMListener[name]))}
Swagger.field_name(::Type{ WinRMListener }, property_name::Symbol) =  _property_map_WinRMListener[property_name]

const _allowed_WinRMListener_protocol = ["Http", "Https"]

function check_required(o::WinRMListener)
    true
end

function validate_property(::Type{ WinRMListener }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "WinRMListener", :enum, val, _allowed_WinRMListener_protocol)
    end
end


text:
@testset "Convert continuous to discrete type" begin
    for dtype in subtypes(AbstractDiscreteSystem)
        ctype = eval.(Meta.parse.(replace(string(dtype), "Discrete" => "Continuous")))
        @test MathematicalSystems._complementary_type(dtype) == ctype
        @inferred MathematicalSystems._complementary_type(dtype)
    end

    for ctype in subtypes(AbstractContinuousSystem)
        dtype = eval.(Meta.parse.(replace(string(ctype), "Continuous" => "Discrete")))
        @test MathematicalSystems._complementary_type(ctype) == dtype
        @inferred MathematicalSystems._complementary_type(ctype)
    end
end

@testset "Exact discretization of affine continous systems" begin
    algorithm = MathematicalSystems.ExactDiscretization()
    ΔT = 0.1
    A = [0.5 1; 0.0 0.5]
    B = Matrix([0.5 1.0]')
    c = [1.0; 1.0]
    D = [0.3 0.7; -0.5 1.30]
    X = BallInf(zeros(2), 1.0)
    U = BallInf(zeros(1), 2.0)
    W = BallInf(zeros(2), 3.0)
    A_d = exp(A*ΔT)
    M = inv(A)*(A_d - I)
    B_d = M*B
    c_d = M*c
    D_d = M*D
    dict = Dict([:A => A, :B => B, :b => c, :c => c, :D => D,
                 :X => X, :U => U, :W => W])
    dict_discretized = Dict([:A => A_d, :B => B_d, :b => c_d, :c => c_d, :D => D_d,
                 :X => X, :U => U, :W => W])

    # get affine ctypes
    CTYPES = filter(x -> (occursin("Linear", string(x)) || occursin("Affine", string(x))) &&
                         !occursin("Algebraic", string(x)) , subtypes(AbstractContinuousSystem))

    # this test doesn't apply for second order systems
    filter!(x -> x ∉ SECOND_ORDER_CTYPES, CTYPES)

    DTYPES = MathematicalSystems._complementary_type.(CTYPES)
    n_types = length(CTYPES)
    CFIELDS = fieldnames.(CTYPES)
    CValues = [getindex.(Ref(dict), type) for type in CFIELDS]
    DValues = [getindex.(Ref(dict_discretized), type) for type in CFIELDS]
    # Create discretized system manually
    discretized_manually = [DTYPES[i](DValues[i]...) for i=1:n_types]
    for i=1:n_types
        discretized_with_method = discretize(CTYPES[i](CValues[i]...), ΔT, algorithm)
        @test discretized_with_method == discretized_manually[i]
    end
end

@testset "Euler discretization of affine continous systems" begin
    algorithm = MathematicalSystems.EulerDiscretization()
    ΔT = 0.1
    A = [0.5 1; 0.0 0.5]
    B = Matrix([0.5 1.0]')
    c = [1.0; 1.0]
    D = [0.3 0.7; -0.5 1.30]
    X = BallInf(zeros(2), 1.0)
    U = BallInf(zeros(1), 2.0)
    W = BallInf(zeros(2), 3.0)
    A_d = I + ΔT*A
    B_d = ΔT*B
    c_d = ΔT*c
    D_d = ΔT*D
    dict = Dict([:A => A, :B => B, :b => c, :c => c, :D => D,
                 :X => X, :U => U, :W => W])
    dict_discretized = Dict([:A => A_d, :B => B_d, :b => c_d, :c => c_d, :D => D_d,
                 :X => X, :U => U, :W => W])

    # get affine ctypes
    CTYPES = filter(x -> (occursin("Linear", string(x)) || occursin("Affine", string(x))) &&
                         !occursin("Algebraic", string(x)) , subtypes(AbstractContinuousSystem))

    # this test doesn't apply for second order systems
    filter!(x -> x ∉ SECOND_ORDER_CTYPES, CTYPES)

    DTYPES = MathematicalSystems._complementary_type.(CTYPES)
    n_types = length(CTYPES)
    CFIELDS = fieldnames.(CTYPES)
    CValues = [getindex.(Ref(dict), type) for type in CFIELDS]
    DValues = [getindex.(Ref(dict_discretized), type) for type in CFIELDS]
    # Create discretized system manually
    discretized_manually = [DTYPES[i](DValues[i]...) for i=1:n_types]
    for i=1:n_types
        discretized_with_method = discretize(CTYPES[i](CValues[i]...), ΔT, algorithm)
        @test discretized_with_method == discretized_manually[i]
    end
end


text:
@testset "is_raw_ssh_private_key" begin
    content = "- BEGIN END PRIVATE KEY"

    @test CompatHelper.is_raw_ssh_private_key(content)
end

@testset "decode_ssh_private_key" begin
    expected = "- BEGIN END PRIVATE KEY"

    @testset "raw key" begin
        @test CompatHelper.decode_ssh_private_key(expected) == expected
    end

    @testset "base64 encoded" begin
        encoded = Base64.base64encode(expected)

        @test CompatHelper.decode_ssh_private_key(encoded) == expected
    end

    @testset "non-base64 encoded" begin
        encoded = bytes2hex(SHA.sha256(expected))

        @test_throws CompatHelper.UnableToParseSSHKey CompatHelper.decode_ssh_private_key(
            encoded
        )
    end
end


text:
# This file is a part of TextClassification.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

import KCenters: center
export TextCentroidSelection

struct TextCentroidSelection <: AbstractCenterSelection end

function center(::TextCentroidSelection, lst::AbstractVector{<:DVEC})
    u = zero(eltype(lst))

    for v in lst
        TextSearch.add!(u, v)
    end

    TextSearch.normalize!(u)
end


text:
module SloanDigitalSkySurvey

include("PSF.jl")  # PSF must come first because SDSS imports it.
include("SDSS.jl")
include("WCSUtils.jl")

end # module


text:
y = closeness_centrality(g5; normalize=false)
z = closeness_centrality(g5)

@test y == [0.75, 0.6666666666666666, 1.0, 0.0]
@test z == [0.75, 0.4444444444444444, 0.3333333333333333, 0.0]

g = Graph(5)
add_edge!(g,1,2)
z = closeness_centrality(g)
@test z[1] == z[2] == 0.25
@test z[3] == z[4] == z[5] == 0.0


text:
function complementary(v::Array{Float64}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Vector{Float64}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Vector{Int}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Int64, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Array{Int64}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

max_stars_repo_name: kahaaga/SimplexSplitting

max_issues_repo_name: kahaaga/SimplexSplitting

max_forks_repo_name: kahaaga/SimplexSplitting


text:
import DSP, PortAudio

MonoAudioSlice = SampleBuf{T, 1} where T <: AbstractFloat
StereoAudioSlice = SampleBuf{T, 2} where T <: AbstractFloat
AudioSlice = Union{StereoAudioSlice, MonoAudioSlice}

# abstract type AudioRepresentationSlice <: Union{AudioSlice, SpectrumSlice} end

# abstract type AudioStream <: PortAudioStream{T} where T <: AbstractFloat end
# abstract type AudioFeature <: Any end

AudioFeature = Any


text:
using SFML

window = RenderWindow("RenderText example", 800, 600)
set_framerate_limit(window, 60)

event = Event()

text = RenderText()
set_position(text, Vector2f(400.0, 300.0))
set_string(text, "SFML.jl is cool")
set_style(text, TextBold | TextUnderlined | TextItalic | TextStrikeThrough)
set_color(text, SFML.cyan)
set_charactersize(text, 50)
text_size = get_globalbounds(text)
set_origin(text, Vector2f(text_size.width / 2, text_size.height / 2))

println(get_string(text))

while isopen(window)
    while pollevent(window, event)
        if get_type(event) == EventType.CLOSED
            close(window)
        end
    end

    clear(window, SFML.black)
    draw(window, text)
    display(window)
end


text:
using Flux:Dense,Chain
using Flux:destructure,params
using Zygote:Params,Grads
using Zygote:gradient

mutable struct NNValueBaseline <: AbstractValueBaseline
  pars::Params
  pinfo::Vector{Tuple{Int64, Vararg{Int64, N} where N}}
  nn::Union{Chain,Dense}
  re::Any
  theta::Vector
  function NNValueBaseline(nn::Union{Chain,Dense})
    pars = params(nn)
    pinfo = [size(p_) for p_ in pars]
    theta,re = destructure(nn)
    new(pars,pinfo,nn,re,theta)
  end
end

function gradshape(vbl::NNValueBaseline)
  return size(vbl.theta)[1]
end

function callvbl(sys::AbstractSystem,vbl::NNValueBaseline)
  vout::Float64 = vbl.nn(push!(sys.x,sys.t))[1]
  pop!(sys.x)
  return vout
end

"""
Gradient wrt parameters for gaussian value baseline
"""
function gradient!(grad::Vector{Float64},sys::AbstractSystem,vbl::NNValueBaseline)
  # compute gradients of output wrt parametes
  push!(sys.x,sys.t)
  grads::Grads = gradient(() -> vbl.nn(sys.x)[1], vbl.pars)
  pop!(sys.x)
  # store gradients into jacobian array
  bcount::Int64 = 1
  for i in 1:length(vbl.pars)
    # TODO replace with reduce function
    ptup::Tuple = vbl.pinfo[i]
    ecount::Int64 = 1
    for j in 1:length(ptup)
      ecount *= ptup[j]
    end
    ecount += (bcount - 1)
    @views grad[bcount:ecount] .= vec(grads[vbl.pars[i]])
    bcount = ecount + 1
  end
end

"""
Update parameters for the neural network value baseline
"""
function updateparams!(lr::Float64,dtheta::Vector{Float64},vbl::NNValueBaseline)
    # update coefficients for forces
    vbl.theta .+= (lr .* dtheta)
    setparams!(vbl)
end

"""
Set parameters for the neural network model
"""
function setparams!(vbl::NNValueBaseline)
    # update coefficients for forces
    vbl.nn = vbl.re(vbl.theta)
    vbl.pars = params(vbl.nn)
end

function setparams!(theta::Vector{Float64},vbl::NNModel)
    # update coefficients for forces
    vbl.theta .= theta
    setparams!(vbl)
end


text:
using Printf;
using Plots;

function init()
	## 横の長さ
	global lx = 1.0;
	## 縦の長さ
	global ly = 1.0;
	## 横の分割数
	global nx = 11;
	## 縦の分割数
	global ny = 11;
	global dx = lx / (nx - 1);
	global dy = ly / (ny - 1);
	## 係数
	beta_x = 1.0 / dx^2;
	beta_y = 1.0 / dy^2;
	global x = range(0.0, stop=lx, length=nx);
	global y = range(0.0, stop=ly, length=ny);
	global u = zeros(Float64, (ny, nx));
	## 境界条件の適用
	u[:, begin] .= 0.0;
	u[:, end] = 100y;
	u[begin, :] .= 0.0;
	u[end, :] = 100x;
	global u_t = zeros(Float64, (ny, nx));
	for j = 1:ny
		for i = 1:nx
			u_t[j, i] = 100*x[i]*y[j];
		end
	end
	global pre_u = copy(u);
	## SOR法のωの値を決定
	mu = cos(pi/nx) + cos(pi/ny);
	global omega = 2.0/(1.0 + sqrt(1-(mu/2)^2));
	global epsilon = 2e-4;
	global a_ir = beta_x;
	global a_il = beta_x;
	global a_ia = beta_y;
	global a_ib = beta_y;
	global a_0 = -2.0(beta_x + beta_y);
	global rhs = 0.0;
	global iter_max = 100;
	global epsilon = 2e-4;
end

function main()
	init();
	for iter=1:iter_max
		err = 0.0;
		for j = 1:ny-2
			for i = 1:nx-2
				tmp = (rhs - a_ir*u[j+1, i+2] - a_il*u[j+1, i] - a_ia*u[j+2, i+1] - a_ib*u[j, i+1]) / a_0;
				u[j+1, i+1] = u[j+1, i+1] + (tmp - u[j+1, i+1])*omega;;
				err += sqrt((u_t[j+1, i+1] - u[j+1, i+1])^2)/(nx*ny);
			end
		end
		err = sqrt(sum((u - u_t).^2)/(nx*ny))
		#err = sqrt(sum((u - pre_u).^2)/(nx*ny)) global pre_u = copy(u);
		@printf("iter: %d, err: %5e\n", iter, err);
		if err < epsilon
			break;
		end
	end
	contour(u);
	savefig("Laplace.png");
end

main()


text:
"""
    Transducer(o::OnlineStat)

Use an `OnlineStat` as a stateful transducer.

It is implemented as:

```julia
opcompose(Scan(fit!, CopyInit(o)), Map(value))
```

# Examples
```jldoctest
julia> using Transducers
       using OnlineStats: Mean

julia> collect(Transducer(Mean()), 1:4)
4-element Vector{Float64}:
 1.0
 1.5
 2.0
 2.5
```
"""
Transducer(o::OnlineStatsBase.OnlineStat) =
    Map(OnlineStatsBase.value) ∘ Scan(OnlineStatsBase.fit!, CopyInit(o))
# TODO: implement `combine`

"""
    reducingfunction([xf::Transducer,] o::OnlineStat; simd)

Convert an `OnlineStat` to a reducing function.  Returned function can
be used with [`foldl`](@ref), [`foldxt`](@ref), and [`foldxd`](@ref).
Note that input `o` is only used as a "prototype"; i.e., it's not
going to be mutated.

# Examples
```jldoctest
julia> using Transducers
       using OnlineStats: Mean

julia> foldl(reducingfunction(Mean()), Map(x -> x^2), 1:4)
Mean: n=4 | value=7.5

julia> foldl(Mean(), Map(x -> x^2), 1:4)  # equivalent to above
Mean: n=4 | value=7.5

julia> foldxt(Mean(), Map(x -> x^2), 1:4)  # threaded
Mean: n=4 | value=7.5

julia> foldl(Mean(), eduction(x^2 for x in 1:4))  # ditto
Mean: n=4 | value=7.5

julia> foldxt(Mean(), eduction(x^2 for x in 1:4))  # ditto
Mean: n=4 | value=7.5
```

[`foldxd`](@ref) can be used instead of `foldxt`.  However the usual caveats
of code availability for Distributed.jl apply.
"""
reducingfunction(xf::Transducer, stat::OnlineStatsBase.OnlineStat; kwargs...) =
    reducingfunction(xf, reducingfunction(stat); kwargs...)
reducingfunction(stat::OnlineStatsBase.OnlineStat) = OnlineStatReducingFunction(stat)

@inline _reducingfunction(xf, stat::OnlineStatsBase.OnlineStat; kwargs...) =
    _reducingfunction(xf, reducingfunction(stat); kwargs...)

struct OnlineStatReducingFunction{T<:OnlineStatsBase.OnlineStat} <: Function
    init::T
end

InitialValues.hasinitialvalue(::Type{<:OnlineStatReducingFunction}) = true

# `GenericInitialValue{T}` for case when `init=Init(op::OnlineStat)`
# is explicitly passed.
const OSInit{T} =
    Union{GenericInitialValue{T},GenericInitialValue{OnlineStatReducingFunction{T}}}

(rf::OnlineStatReducingFunction{T})(::OSInit{T}, x) where {T<:OnlineStatsBase.OnlineStat} =
    rf(copy(rf.init), x)

function (::OnlineStatReducingFunction)(acc, x)
    OnlineStatsBase.fit!(acc, x)
    return acc
end

combine(::OnlineStatReducingFunction{T}, a::T, b::T) where {T<:OnlineStatsBase.OnlineStat} =
    merge!(a, b)

combine(
    ::OnlineStatReducingFunction{T},
    a::Union{T,OSInit{T}},
    ::OSInit{T},
) where {T<:OnlineStatsBase.OnlineStat} = a

combine(
    ::OnlineStatReducingFunction{T},
    ::OSInit{T},
    b::T,
) where {T<:OnlineStatsBase.OnlineStat} = b

const OSNonZeroNObsError = ArgumentError(
    "An `OnlineStat` with one or more observations cannot be used with " *
    "`foldxt` and `foldxd`.",
)

function validate_reduce_ostat(stat)
    if OnlineStatsBase.nobs(stat) != 0
        throw(OSNonZeroNObsError)
    end
    return stat
end

# Method plumbing:
foldxl(stat::OnlineStatsBase.OnlineStat, xform::Transducer, itr; kwargs...) =
    foldxl(reducingfunction(stat), xform, itr; kwargs...)

foldxl(stat::OnlineStatsBase.OnlineStat, foldable; kwargs...) =
    foldxl(reducingfunction(stat), extract_transducer(foldable)...; kwargs...)

foldxt(stat::OnlineStatsBase.OnlineStat, xform::Transducer, itr; kwargs...) =
    foldxt(reducingfunction(validate_reduce_ostat(stat)), xform, itr; kwargs...)

foldxt(stat::OnlineStatsBase.OnlineStat, foldable; kwargs...) = foldxt(
    reducingfunction(validate_reduce_ostat(stat)),
    extract_transducer(foldable)...;
    kwargs...,
)

foldxd(stat::OnlineStatsBase.OnlineStat, xform::Transducer, itr; kwargs...) =
    foldxd(reducingfunction(validate_reduce_ostat(stat)), xform, itr; kwargs...)

foldxd(stat::OnlineStatsBase.OnlineStat, foldable; kwargs...) = foldxd(
    reducingfunction(validate_reduce_ostat(stat)),
    extract_transducer(foldable)...;
    kwargs...,
)


text:
# This file is a part of UnsafeArrays.jl, licensed under the MIT License (MIT).

using UnsafeArrays
using Test

using Random

@testset "uview" begin
    function rand_array(::Type{T}, Val_N::Val{N}) where {T, N}
        sz_max = (8, 7, 5, 4, 5)
        sz = ntuple(i -> sz_max[i], Val_N)
        rand!(Array{T}(undef, sz...))
    end

    function rand_array(::Type{T}, Val_N::Val{N}) where {T<:Integer, N}
        sz_max = (8, 7, 5, 4, 5)
        sz = ntuple(i -> sz_max[i], Val_N)
        rand!(Array{T}(undef, sz...), 0:99)
    end

    function test_A(test_code::Function, ::Type{T}, Val_N::Val{N}) where {T, N}
        A = rand_array(T, Val_N)
        UnsafeArrays.@gc_preserve A test_code(A)
    end

    @testset "uview" begin
        test_A(Float64, Val(0)) do A
            @test typeof(@inferred(uview(A))) == UnsafeArray{Float64,0}

            @test uview(A) == A
            @test uview(A) == view(A)
        end

        test_A(Float32, Val(1)) do A
            @test typeof(@inferred(uview(A))) == UnsafeArray{Float32,1}
            @test typeof(@inferred(uview(A, :))) == UnsafeArray{Float32,1}
            @test typeof(@inferred(uview(A, 2))) == UnsafeArray{Float32,0}

            @test uview(A) == A
            @test uview(A, :) == view(A, :)
            @test uview(A, 2) == view(A, 2)

            @test_throws BoundsError uview(A, 2:9)
        end

        test_A(ComplexF32, Val(3)) do A
            @test typeof(@inferred(uview(A))) == UnsafeArray{ComplexF32,3}
            @test typeof(@inferred(uview(A, :, 2:4, 3))) == UnsafeArray{ComplexF32,2}
            @test typeof(@inferred(uview(A, :))) == UnsafeArray{ComplexF32,1}
            @test typeof(@inferred(uview(A, 2, 3, 4))) == UnsafeArray{ComplexF32,0}
            # TODO: Type inference fails for some reason:
            # @test typeof(@inferred(uview(A, :, 2:4, :))) <: SubArray
            @test typeof((uview(A, :, 2:4, :))) <: SubArray

            @test uview(A) == A
            @test uview(A, :, 2:4, 3) == view(A, :, 2:4, 3)
            @test uview(A, :) == view(A, :)
            @test uview(A, 2, 3, 4) == view(A, 2, 3, 4)
            @test uview(A, :, 2:4, :) == view(A, :, 2:4, :)

            @test_throws BoundsError uview(A, :, 2:9, 3)

            B = view(A, :, :, 3)

            @test typeof(@inferred(uview(B))) == UnsafeArray{ComplexF32,2}
            @test typeof(@inferred(uview(B, :, 2:4))) == UnsafeArray{ComplexF32,2}

            @test uview(B) == B
            @test uview(B, :, 2:4) == view(B, :, 2:4)

            C = uview(A, :, :, 3)

            @test uview(C) === C
            @test typeof(@inferred(uview(C, :, 2:4))) == UnsafeArray{ComplexF32,2}
            @test uview(C, :, 2:4) == view(B, :, 2:4)
        end

        let A = ["foo", "bar", "baz"]
            UnsafeArrays.@gc_preserve A begin
                @test @inferred(uview(A)) === A
                @test typeof(@inferred(uview(A, 2:3))) <: SubArray
                @test uview(A, 2:3) == view(A, 2:3)
            end
        end
    end

    @testset "uviews" begin
        A = rand(Int32, 8)
        B = rand(ComplexF64, 3, 5, 4)

        uviews(() -> 42) == 42
        uviews(A -> typeof(A), A) == UnsafeArray{Int32, 1}
        uviews((A, B) -> (typeof(A), typeof(B)), A, B) == UnsafeArray{Int32, 1}
    end

    @testset "@uviews" begin
        A = rand(Int32, 8)
        B = rand(ComplexF64, 3, 5, 4)

        @uviews(42) == 42
        @uviews(A, typeof(A)) == UnsafeArray{Int32, 1}
        @uviews(A, B, (typeof(A), typeof(B))) == UnsafeArray{Int32, 1}
    end

    @testset "Base.mightalias" begin
        A = rand(Int32, 10, 15)
        B = rand(Int32, 12, 17)
        C = ["foo", "bar", "a", "b", "c", "d", "e"]

        @uviews A B begin
            @test true == @inferred Base.mightalias(A, A)
            @test false == @inferred Base.mightalias(A, B)

            @test false == @inferred Base.mightalias(view(A, :, 2:5), view(A, :, 6:10))
            @test false == @inferred Base.mightalias(view(A, :, 6:10), view(A, :, 2:5))

            @test true == @inferred Base.mightalias(view(A, :, 2:6), view(A, :, 5:10))
            @test true == @inferred Base.mightalias(view(A, :, 5:10), view(A, :, 2:6))
            @test true == @inferred Base.mightalias(view(A, :, 2:10), view(A, :, 5:6))
            @test true == @inferred Base.mightalias(view(A, :, 5:6), view(A, :, 2:10))

            @test false == @inferred Base.mightalias(view(A, :, 2:5), view(A, :, 6:10))
            @test true == @inferred Base.mightalias(view(A, 2:5, :), view(A, :, 6:10))
            @test true == @inferred Base.mightalias(view(A, :, 2:5), view(A, 6:10, :))
            @test true == @inferred Base.mightalias(view(A, 2:5, :), view(A, 6:10, :))
        end

        @uviews A begin
            @test false == @inferred Base.mightalias(A, B)
            @test false == @inferred Base.mightalias(view(A, 2:5, :), B)
            @test false == @inferred Base.mightalias(A, view(B, :, 2:5))
            @test false == @inferred Base.mightalias(view(A, 2:5, :), view(B, :, 2:5))

            @test false == @inferred Base.mightalias(B, A)
            @test false == @inferred Base.mightalias(B, view(A, 2:5, :))
            @test false == @inferred Base.mightalias(view(B, :, 2:5), A)
            @test false == @inferred Base.mightalias(view(B, :, 2:5), view(A, 2:5, :))
        end

        @uviews A C begin
            @test false == @inferred Base.mightalias(A, C)

            @test false == @inferred Base.mightalias(view(A, 2:5, :), C)
            @test false == @inferred Base.mightalias(A, view(C, 2:5))
            @test false == @inferred Base.mightalias(view(A, 2:5, :), view(C, 2:5))
        end

        @uviews A begin
            @test false == @inferred Base.mightalias(A, C)
            @test false == @inferred Base.mightalias(view(A, 2:5, :), C)
            @test false == @inferred Base.mightalias(A, view(C, 2:5))
            @test false == @inferred Base.mightalias(view(A, 2:5, :), view(C, 2:5))

            @test false == @inferred Base.mightalias(C, A)
            @test false == @inferred Base.mightalias(C, view(A, 2:5, :))
            @test false == @inferred Base.mightalias(view(C, 2:5), A)
            @test false == @inferred Base.mightalias(view(C, 2:5), view(A, 2:5, :))
        end
    end
end


text:
module NVML

using ..CUDA
import ..CUDA: ndevices, getdevice
using Libdl

if Sys.iswindows()
    const libnvml = Libdl.find_library("nvml", [joinpath(ENV["ProgramFiles"],"NVIDIA Corporation","NVSMI")])
else
    const libnvml = Libdl.find_library("libnvidia-ml")
end
isempty(libnvml) && error("NVML cannot be found.")

const API_VERSION = Ref{String}()

function checkresult(result::Cint)
    if result != 0
        p = ccall((:nvmlErrorString,libnvml), Ptr{Cchar}, (Cint,), result)
        throw(unsafe_string(p))
    end
end

function __init__()
    result = ccall((:nvmlInit_v2,libnvml), Cint, ())
    checkresult(result)

    ref = Array{Cchar}(undef, 80)
    result = ccall((:nvmlSystemGetNVMLVersion,libnvml), Cint, (Ptr{Cchar},Cuint), ref, 80)
    checkresult(result)

    API_VERSION[] = unsafe_string(pointer(ref))
    @info "NVML $(API_VERSION[])"
end

include("define.jl")

macro nvml(f, args...)
    f = get(DEFINE, f.value, f.value)
    quote
        result = ccall(($(QuoteNode(f)),libnvml), Cint, $(map(esc,args)...))
        checkresult(result)
    end
end

macro nvml_nocheck(f, args...)
    f = get(DEFINE, f.value, f.value)
    quote
        ccall(($(QuoteNode(f)),libnvml), Cint, $(map(esc,args)...))
    end
end

include("device.jl")

end


text:
"""
    uchain(;encoders, decoders, bridge, connection)

Build a `Chain` with U-Net like architecture. `encoders` and `decoders` are
arrays of encoding/decoding blocks, from top to bottom (see diagram below).
`bridge` is the bottom part of U-Net  architecture.
Each level of the U-Net is connected through a 2-argument callable `connection`.
`connection` could be an array in case the way levels are connected vary from
one level to another.

Notes :
- usually encoder block starts with a 'MaxPool' to downsample image by 2, except
the first level encoder.
- usually decoder block ends with a 'ConvTranspose' to upsample image by 2,
except the first level decoder.

```
+---------+                                                          +---------+
|encoder 1|                                                          |decoder 1|
+---------+                                                          +---------+
     |------------------------------------------------------------------->^
     |                                                                    |
     |   +---------+                                         +---------+  |
     +-->|encoder 2|                                         |decoder 2|--+
         +---------+                                         +---------+
              |-------------------------------------------------->^
              |                                                   |
              |   +---------+                        +---------+  |
              +-->|encoder 3|                        |decoder 3|--+
                  +---------+                        +---------+
                       |--------------------------------->^
                       |                                  |
                       |   +---------+       +---------+  |
                       +-->|encoder 4|       |decoder 4|--+
                           +---------+       +---------+
                                |---------------->^
                                |                 |
                                |    +--------+   |
                                +--->| bridge |---+
                                     +--------+
```
See also [`chcat`](@ref).
"""
function uchain(;encoders, decoders, bridge, connection)
    length(encoders) == length(decoders) || throw(ArgumentError(
        "The number of encoders should be equal to the number of decoders."))

    if isa(connection, AbstractArray)
        length(encoders) == length(connection) || throw(ArgumentError(
            "The number of connections should be equal to the number of encoders/decoders."))
    else
        connection = repeat([connection], length(encoders))
    end

    ite = zip(reverse(encoders[2:end]),
              reverse(decoders[2:end]),
              reverse(connection[1:(end - 1)]))
    l = ubridge(bridge, connection[end])
    for (e, d, c) ∈ ite
        l = uconnect(e, l, d)
        l = SkipConnection(l, c)
    end
    uconnect(encoders[1], l, decoders[1])
end

@inline uconnect(enc, prl, dec) = Chain(enc, prl, dec)

for T1 ∈ [:Chain :AbstractArray], T2 ∈ [:Chain :AbstractArray]
    @eval begin
        @inline uconnect(enc::($T1), prl, dec::($T2)) = Chain(enc..., prl, dec...)
    end
end

for T ∈ [:Chain :AbstractArray]
    @eval begin
        @inline uconnect(enc::($T), prl, dec) = Chain(enc..., prl, dec)
        @inline uconnect(enc, prl, dec::($T)) = Chain(enc, prl, dec...)
    end
end

@inline ubridge(b, c) = SkipConnection(b, c)

@inline ubridge(b::AbstractArray, c) = SkipConnection(Chain(b...), c)


text:
# ode checks
using OrdinaryDiffEq, DiffEqGPU, Test

seed = 100
using Random;Random.seed!(seed)
ra = rand(100)

function f!(du,u,p,t)
 @inbounds begin
     du[1] = 1.01*u[1]
 end
 nothing
end

prob = ODEProblem(f!,[0.5],(0.0,1.0))

function output_func(sol,i)
  last(sol), false
end

function prob_func(prob,i,repeat)
  remake(prob,u0=ra[i]*prob.u0)
end

function reduction(u,batch,I)
  u.+sum(batch),false
end

# no reduction
prob1 = EnsembleProblem(prob,prob_func=prob_func,output_func=output_func)
sim1 = @time solve(prob1,Tsit5(),trajectories=100,batch_size=20)

# reduction and EnsembleThreads()
prob2 = EnsembleProblem(prob,prob_func=prob_func,output_func=output_func,
  reduction=reduction,u_init=Vector{eltype(prob.u0)}([0.0])
  )
sim2 = @time solve(prob2,Tsit5(),trajectories=100,batch_size=20)

# EnsembleCPUArray() and EnsembleGPUArray()
sim3 = @time solve(prob2,Tsit5(),EnsembleCPUArray(),trajectories=100,batch_size=20)
sim4 = @time solve(prob2,Tsit5(),EnsembleGPUArray(),trajectories=100,batch_size=20)

@info  sim2[1]

@test sum(sim1.u) == sim2.u
@test sim2.u ≈ sim3.u
@test sim2.u == sim4.u


text:
export pick_from_pool, cast

const ACTIONS = Any[]
const NAMES = [
    "Roger M",
    "H Roger",
    "Anna",
    "Alex",
    "Bohdan",
    "Isaac",
    "Ejaaz",
    "Wetzel",
    "Estelle",
    "Danny",
    "Lucy",
    "You",
]

"""
    @action <action function>(name)

Register an action function. An action function must take
a `name` as input and return a sentence of type `String`.
"""
macro action(ex)
    def = splitdef(ex)

    quote
        export $(esc(def[:name]))
        $(esc(ex))
        push!($(GlobalRef(PinGo, :ACTIONS)), $(esc(def[:name])))
        $(esc(def[:name]))
    end
end

"""
    does_or_not(name, action; magic_number=0.432)

`name` does `action` or not. By probability of a `magic_number`.
"""
does_or_not(name, action; magic_number=0.432) = string(name, " ", rand() < magic_number ? "not " : "", action)

@action wear_pants(name="Roger") = does_or_not(name, "wearing pants")
@action drink_beer(name="Roger") = does_or_not(name, "drinking beer")
@action squats(name="Roger") = does_or_not(name, "doing squats")
@action pushups(name="Roger") = does_or_not(name, "doing push-ups")
@action cast_spells(name="Roger") = does_or_not(name, "casting spells")
@action learned_so_much(name="Roger") = does_or_not(name * " have", "learned so much")

const BINGO_POOL = readlines(joinpath(@__DIR__, "pool.txt"))

"""
    pick_from_pool()

Randomly pick a phrase from the pool.
"""
pick_from_pool() = rand(BINGO_POOL)

function cast(p=0.8)
    dice = rand()
    if dice < p
        return pick_from_pool()
    else
        rand(ACTIONS)(rand(NAMES))
    end
end


text:
@testset "Common" begin
end


text:
module GeographicData

using CSV
using Dates
using DataFrames
using LinearAlgebra: diagm

export adjacencies
export get_counties, get_county_names
export haversine_distance_matrix

basepath = normpath(@__DIR__, "../../")

########################
######### Data #########
########################

function load_states_google(;states::Array{String,1}=String[], method::Symbol=:cities)::Tuple{Array{Int,2},Array{String,1}}
	@assert method in [:states, :cities]
	if (method == :states) fn = "state_dist_google"
	else fn = "state_cities_dist_google"
	end

	dist_df = DataFrame(CSV.File(joinpath(basepath, "data/geography/$(fn).csv")), copycols=true)

	if (method == :cities) state_abbrev_list = map(s -> String(s)[end-1:end], names(dist_df))
	else state_abbrev_list = names(dist_df)[:]
	end

	if !isempty(states)
		@assert states == sort(states)
		ind = [findfirst(x -> x == s, state_abbrev_list) for s in states]
	else
		ind = sortperm(state_abbrev_list)
	end

	state_abbrev_list = state_abbrev_list[ind]
	dist_matrix = Matrix(dist_df)[ind, ind]

	return dist_matrix, state_abbrev_list
end

function load_states_latlong(;states::Array{String,1}=String[])::Tuple{Array{Float64,2},Array{String,1}}
	state_data = DataFrame(CSV.File(joinpath(basepath, "data/geography/states.csv")), copycols=true)
	if !isempty(states) filter!(row -> row.abbrev in states, state_data) end
	sort!(state_data, :abbrev)
	latlong = hcat(state_data.lat[:], state_data.long[:])
	state_abbrev = state_data.abbrev[:]
	return latlong, state_abbrev
end

function load_counties_latlong(;counties::Array{Int,1}=Int[])::Tuple{Array{Float64,2},Array{Int,1}}
	county_data = DataFrame(CSV.File(joinpath(basepath, "data/geography/counties.csv")), copycols=true)
	if !isempty(counties)
		@assert counties == sort(counties)
		filter!(row -> row.fips in counties, county_data)
	end
	sort!(county_data, :fips)
	latlong = hcat(county_data.lat[:], county_data.long[:])
	county_fips = county_data.fips[:]
	return latlong, county_fips
end

########################
###### Adjacency #######
########################

function adjacencies(locations::Array; level::Symbol=:state, source::Symbol=:google, threshold::Real=4, self_edges::Bool=false)::BitArray{2}
	if source == :fullyconnected
		return fully_connected(length(locations), self_edges=self_edges)
	elseif level == :state && source == :google
		dist_matrix, all_states = load_states_google(states=locations)
		@assert all_states == locations
	elseif level == :state && source == :original
		error("Not yet implemented.")
	elseif level == :state && source == :latlong
		latlong, all_states = load_states_latlong(states=locations)
		@assert all_states == locations
		dist_matrix = haversine_distance_matrix(latlong)
	elseif level == :county
		latlong, all_fips = load_counties_latlong(counties=locations)
		@assert all_fips == locations
		dist_matrix = haversine_distance_matrix(latlong)
	elseif level == :hospital
		dist_matrix = haversine_distance_matrix(locations)
	elseif level == :other
		dist_matrix = haversine_distance_matrix(locations)
	else
		error("Invalid parameters to compute_adjacencies.")
	end

	if self_edges
		adj_matrix = 0 .<= dist_matrix .<= 3600 * threshold
	else
		adj_matrix = 0 .<  dist_matrix .<= 3600 * threshold
	end

	return adj_matrix
end

########################
###### Haversine #######
########################

function haversine_distance(lat1, lon1, lat2, lon2; dist_type=:time, speed_kph=100)
	R = 6371e3

	φ1 = lat1 * π/180
	φ2 = lat2 * π/180
	Δφ = (lat2-lat1) * π/180
	Δλ = (lon2-lon1) * π/180

	a = (sin(Δφ/2) * sin(Δφ/2)) + (cos(φ1) * cos(φ2) * sin(Δλ/2) * sin(Δλ/2))
	c = 2 * atan(sqrt(a), sqrt(1-a))

	dist = R * c

	if dist_type == :time
		dist = dist / 1000 / speed_kph * 3600
	end

	return dist
end

function haversine_distance_matrix(locations::Array{Float64,2}; dist_type=:time, speed_kph=100)
	N = size(locations, 1)
	distancematrix = zeros(Float32, N, N)
	for i in 1:N
		for j in i+1:N
			dist = haversine_distance(locations[i,:]..., locations[j,:]..., dist_type=dist_type, speed_kph=speed_kph)
			distancematrix[i,j] = dist
			distancematrix[j,i] = dist
		end
	end
	return distancematrix
end

########################
######## Other #########
########################

function fully_connected(n::Int; self_edges::Bool=false)
	if (self_edges) return BitArray(ones(Bool, n, n)) end
	adj = BitArray(ones(Bool, n, n) - diagm(ones(Bool, n)))
	return adj
end

function get_counties(states::Array{String,1})
	@assert states == sort(states)
	county_data = DataFrame(CSV.File(joinpath(basepath, "data/geography/counties.csv")), copycols=true)
	filter!(row -> row.state_abbrev in states, county_data)
	sort!(county_data, :state_abbrev)
	return county_data.fips[:]
end

function get_county_names(counties::Array{Int,1})
	@assert counties == sort(counties)
	county_data = DataFrame(CSV.File(joinpath(basepath, "data/geography/counties.csv")), copycols=true)
	filter!(row -> row.fips in counties, county_data)
	sort!(county_data, :fips)
	names = map(c -> c.county * " " * c.state_abbrev, eachrow(county_data))
	return names
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
function setup_rfmip_as(
    ds_lw_in,
    idx_gases,
    exp_no,
    lookup,
    ::Type{FT},
    ::Type{DA},
    ::Type{VMR},
    max_threads,
) where {FT <: AbstractFloat, DA, VMR}
    FTA1D = DA{FT, 1}
    FTA2D = DA{FT, 2}

    deg2rad = FT(π) / FT(180)
    nlay = Int(ds_lw_in.dim["layer"])
    ncol = Int(ds_lw_in.dim["site"])
    nlev = nlay + 1
    ngas = lookup.n_gases
    nbnd = lookup.n_bnd

    lon = DA{FT, 1}(ds_lw_in["lon"][:])
    lat = DA{FT, 1}(ds_lw_in["lat"][:])

    lon = nothing # This example skips latitude dependent gravity computation
    lat = nothing # to be consistent with the FORTRAN RRTMGP test case.

    sfc_emis = DA{FT, 2}(repeat(reshape(Array{FT}(ds_lw_in["surface_emissivity"][:]), 1, :), nbnd, 1)) # all bands use same emissivity
    sfc_alb = DA{FT, 2}(repeat(reshape(Array{FT}(ds_lw_in["surface_albedo"][:]), 1, :), nbnd, 1)) # all bands use same albedo
    #--------------------------------------------------------------
    zenith = Array{FT, 1}(deg2rad .* ds_lw_in["solar_zenith_angle"][:])
    irrad = Array{FT, 1}(ds_lw_in["total_solar_irradiance"][:])
    # block out coluumns with zenith > π/2
    usecol = BitArray(undef, ncol)
    usecol .= 1
    for i in 1:ncol
        if zenith[i] > FT(π) / 2 - 2 * eps(FT)
            irrad[i] = FT(0)
            usecol[i] = 0
        end
    end

    zenith = DA{FT, 1}(zenith)
    irrad = DA{FT, 1}(irrad)
    #--------------------------------------------------------------

    p_lev = ds_lw_in["pres_level"][:]

    lev_ind = p_lev[1, 1] > p_lev[end, 1] ? (1:nlev) : (nlev:-1:1)
    lay_ind = p_lev[1, 1] > p_lev[end, 1] ? (1:nlay) : (nlay:-1:1)

    p_lev[lev_ind[end], :] .= lookup.p_ref_min

    p_lev = DA{FT, 2}(p_lev[lev_ind, :])
    p_lay = DA{FT, 2}(ds_lw_in["pres_layer"][:][lay_ind, :])
    t_lev = DA{FT, 2}(ds_lw_in["temp_level"][:][lev_ind, :, exp_no])
    t_lay = DA{FT, 2}(ds_lw_in["temp_layer"][:][lay_ind, :, exp_no])

    t_sfc = DA{FT, 1}(ds_lw_in["surface_temperature"][:, exp_no])
    col_dry = DA{FT, 2}(undef, nlay, ncol)

    # Reading volume mixing ratios

    vmr_h2o = FTA2D(ds_lw_in["water_vapor"][:][lay_ind, :, exp_no]) # vmr of H2O and O3
    vmr_o3 = FTA2D(ds_lw_in["ozone"][:][lay_ind, :, exp_no])       # vary with height

    vmrat = zeros(FT, ngas)

    vmrat[idx_gases["co2"]] =
        FT(ds_lw_in["carbon_dioxide_GM"][exp_no]) * parse(FT, ds_lw_in["carbon_dioxide_GM"].attrib["units"])

    vmrat[idx_gases["n2o"]] =
        FT(ds_lw_in["nitrous_oxide_GM"][exp_no]) * parse(FT, ds_lw_in["nitrous_oxide_GM"].attrib["units"])

    vmrat[idx_gases["co"]] =
        FT(ds_lw_in["carbon_monoxide_GM"][exp_no]) * parse(FT, ds_lw_in["carbon_monoxide_GM"].attrib["units"])

    vmrat[idx_gases["ch4"]] = FT(ds_lw_in["methane_GM"][exp_no]) * parse(FT, ds_lw_in["methane_GM"].attrib["units"])

    vmrat[idx_gases["o2"]] = FT(ds_lw_in["oxygen_GM"][exp_no]) * parse(FT, ds_lw_in["oxygen_GM"].attrib["units"])

    vmrat[idx_gases["n2"]] = FT(ds_lw_in["nitrogen_GM"][exp_no]) * parse(FT, ds_lw_in["nitrogen_GM"].attrib["units"])

    vmrat[idx_gases["ccl4"]] =
        FT(ds_lw_in["carbon_tetrachloride_GM"][exp_no]) * parse(FT, ds_lw_in["carbon_tetrachloride_GM"].attrib["units"])

    vmrat[idx_gases["cfc11"]] = FT(ds_lw_in["cfc11_GM"][exp_no]) * parse(FT, ds_lw_in["cfc11_GM"].attrib["units"])

    vmrat[idx_gases["cfc12"]] = FT(ds_lw_in["cfc12_GM"][exp_no]) * parse(FT, ds_lw_in["cfc12_GM"].attrib["units"])

    vmrat[idx_gases["cfc22"]] = FT(ds_lw_in["hcfc22_GM"][exp_no]) * parse(FT, ds_lw_in["hcfc22_GM"].attrib["units"])

    vmrat[idx_gases["hfc143a"]] = FT(ds_lw_in["hfc143a_GM"][exp_no]) * parse(FT, ds_lw_in["hfc143a_GM"].attrib["units"])

    vmrat[idx_gases["hfc125"]] = FT(ds_lw_in["hfc125_GM"][exp_no]) * parse(FT, ds_lw_in["hfc125_GM"].attrib["units"])

    vmrat[idx_gases["hfc23"]] = FT(ds_lw_in["hfc23_GM"][exp_no]) * parse(FT, ds_lw_in["hfc23_GM"].attrib["units"])

    vmrat[idx_gases["hfc32"]] = FT(ds_lw_in["hfc32_GM"][exp_no]) * parse(FT, ds_lw_in["hfc32_GM"].attrib["units"])

    vmrat[idx_gases["hfc134a"]] = FT(ds_lw_in["hfc134a_GM"][exp_no]) * parse(FT, ds_lw_in["hfc134a_GM"].attrib["units"])

    vmrat[idx_gases["cf4"]] = FT(ds_lw_in["cf4_GM"][exp_no]) * parse(FT, ds_lw_in["hfc23_GM"].attrib["units"])

    #    vmr[idx_gases["no2"]] = FT(ds_lw_in["no2_GM"][exp_no]) *                # missing from input file
    #                                         parse(FT, ds_lw_in["hfc32_GM"].attrib["units"])

    # This example skips latitude dependent gravity compution to be consistent with the
    # FORTRAN RRTMGP test case.
    compute_col_dry!(p_lev, col_dry, param_set, vmr_h2o, lat) # the example skips lat based gravity calculation

    vmr = VMR(vmr_h2o, vmr_o3, FTA1D(vmrat))
    cld_r_eff_liq = nothing
    cld_r_eff_ice = nothing
    cld_path_liq = nothing
    cld_path_ice = nothing
    cld_mask = nothing
    ice_rgh = 1
    #------------------
    return (
        AtmosphericState{
            FT,
            typeof(t_sfc),
            typeof(lat),
            typeof(p_lev),
            typeof(cld_r_eff_liq),
            typeof(cld_mask),
            typeof(vmr),
            Int,
        }(
            lon,
            lat,
            p_lay,
            p_lev,
            t_lay,
            t_lev,
            t_sfc,
            col_dry,
            vmr,
            cld_r_eff_liq,
            cld_r_eff_ice,
            cld_path_liq,
            cld_path_ice,
            cld_mask,
            ice_rgh,
            nlay,
            ncol,
            ngas,
        ),
        sfc_emis,
        sfc_alb,
        zenith,
        irrad,
        usecol,
    )

    #=
    return (
        ClearAtmosphericState{
            FT,
            DA{FT,1},
            typeof(lat),
            DA{FT,2},
            typeof(vmr),
            Int,
        }(
            lon,
            lat,
            p_lay,
            p_lev,
            t_lay,
            t_lev,
            t_sfc,
            col_dry,
            vmr,
            nlay,
            ncol,
            ngas,
        ),
        sfc_emis,
        sfc_alb,
        zenith,
        irrad,
    )
    =#

end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
using ImageInpainting
using ColorTypes
using TestImages
using ReferenceTests
using Test

# list of tests
testfiles = [
  "criminisi.jl"
]

@testset "ImageInpainting.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
abstract type AbstractInputEncoding end
abstract type AbstractEncodingData end

struct RandomMapping{I,T} <: AbstractInputEncoding
    permutations::I
    expansion_size::T
end

"""
    RandomMapping(permutations, expansion_size)
    RandomMapping(permutations; expansion_size=40)
    RandomMapping(;permutations=8, expansion_size=40)

Random mapping of the input data directly in the reservoir. The ```expansion_size```
determines the dimension of the single reservoir, and ```permutations``` determines the
number of total reservoirs that will be connected, each with a different mapping.
The detail of this implementation can be found in [1].

[1] Nichele, Stefano, and Andreas Molund. “Deep reservoir computing using cellular
automata.” arXiv preprint arXiv:1703.02806 (2017).
"""
function RandomMapping(;permutations=8, expansion_size=40)
    RandomMapping(permutations, expansion_size)
end

function RandomMapping(permutations; expansion_size=40)
    RandomMapping(permutations, expansion_size)
end

struct RandomMaps{T,E,G,M,S} <: AbstractEncodingData
    permutations::T
    expansion_size::E
    generations::G
    maps::M
    states_size::S
    ca_size::S
end

function create_encoding(rm::RandomMapping, input_data, generations)
    maps = init_maps(size(input_data, 1), rm.permutations, rm.expansion_size)
    states_size = generations*rm.expansion_size*rm.permutations
    ca_size = rm.expansion_size*rm.permutations
    return RandomMaps(
        rm.permutations, rm.expansion_size, generations, maps, states_size, ca_size
    )
end

function reca_create_states(rm::RandomMaps, automata, input_data)
    train_time = size(input_data, 2)
    states = zeros(rm.states_size, train_time)
    init_ca = zeros(rm.expansion_size*rm.permutations)

    for i in 1:train_time
        init_ca = encoding(rm, input_data[:,i], init_ca)
        ca = CellularAutomaton(automata, init_ca, rm.generations+1)
        ca_states = ca.evolution[2:end ,:]
        states[:,i] = reshape(transpose(ca_states), rm.states_size)
        init_ca = ca.evolution[end, :]
    end

    return states
end

function encoding(rm::RandomMaps, input_vector, tot_encoded_vector)
    input_size = size(input_vector,1)
    #single_encoded_size = Int(size(tot_encoded_vector, 1)/permutations)
    new_tot_enc_vec = copy(tot_encoded_vector)

    for i in 1:rm.permutations
       new_tot_enc_vec[(i-1)*rm.expansion_size+1:i*rm.expansion_size]  = single_encoding(
           input_vector,new_tot_enc_vec[(i-1)*rm.expansion_size+1:i*rm.expansion_size], rm.maps[i,:]
        )
    end

    return new_tot_enc_vec
end

#function obtain_res_size(rm::RandomMapping, generations)
#    generations*rm.expansion_size*rm.permutations
#end

function single_encoding(input_vector, encoded_vector, map)
    new_enc_vec = copy(encoded_vector)

    for i in 1:size(input_vector, 1)
        new_enc_vec[map[i]] = input_vector[i]
    end

    return new_enc_vec
end

function init_maps(input_size, permutations, mapped_vector_size)

    maps = Array{Int}(undef, permutations, input_size)
    #tot_size = input_size*permutations

    for i in 1:permutations
        maps[i,:] = mapping(input_size, mapped_vector_size)
    end

    return maps
end

function mapping(input_size, mapped_vector_size)
    return sample(1:mapped_vector_size, input_size, replace=false)
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork