text:
using ARCTR
using Test
using NLPModels
using JuMP

# first test the LDLt function
include("../src/Utilities/testLDLt.jl")

# test all solvers with the well known Woods test function
include("woods.jl")
nlp = MathProgNLPModel(woods(10), name="woods");

solver = ALL_solvers[25]
(x, f, gNorm, iter, optimal, tired, status) = solver(nlp, verbose=false);
    @printf("%-15s  %8d  %9.2e  %7.1e  %5d  %5d  %6d  %s\n",
            nlp.meta.name, nlp.meta.nvar, f, gNorm,
            nlp.counters.neval_obj, nlp.counters.neval_grad,
            nlp.counters.neval_hprod, status)    #stats = run_solver(solver, model, verbose=false)
    #@test (all([stats...] .>= 0))
    reset!(nlp)
solver = ALL_solvers[26]
(x, f, gNorm, iter, optimal, tired, status) = solver(nlp, verbose=false);
    @printf("%-15s  %8d  %9.2e  %7.1e  %5d  %5d  %6d  %s\n",
            nlp.meta.name, nlp.meta.nvar, f, gNorm,
            nlp.counters.neval_obj, nlp.counters.neval_grad,
            nlp.counters.neval_hprod, status)    #stats = run_solver(solver, model, verbose=false)

#reset!(nlp);


text:
using ArnoldiMethod
import ContinuumArrays.QuasiArrays: AbstractQuasiMatrix
using LazyArrays

using ClassicalOrthogonalPolynomials

function orb_label((n,ℓ))
    ℓs = "spdfg"[ℓ+1]
    "$n$(ℓs)"
end

"""
    hydredwfn(n,ℓ[,Z=1])

Returns a function of `r` (in atomic units, i.e. Bohr radii) that
evaluates the _reduced_ hydrogenic radial orbital corresponding to the
quantum numbers `n` and `ℓ`, and for a nuclear charge of `Z`.
"""
function hydredwfn(n,ℓ,Z=1)
    n > 0 || throw(ArgumentError("Invalid n = $n"))
    0 ≤ ℓ && ℓ<n || throw(ArgumentError("Invalid ℓ = $ℓ"))

    q = 2Z/n
    N = √(q^3*factorial(n-ℓ-1)/(2n*factorial(n+ℓ)))

    r -> N*r*(q*r)^ℓ*laguerrel(n-ℓ-1,2ℓ+1,q*r)*exp(-q*r/2)
end

function get_orbitals(R::B, ℓ::Int, Z, nev::Int,
                      mode::Symbol) where {T,B<:AbstractQuasiMatrix{T}}
    if mode == :arnoldi
        println("Finding eigenstates for Z = $Z, ℓ = $ℓ")
        Tℓ = CoulombIntegrals.get_double_laplacian(R,ℓ,T)
        Tℓ ./= -2
        r = axes(R,1)
        H = Tℓ + R' * QuasiDiagonal(-Z ./ r) * R

        schur,history = partialschur(H, nev=nev, tol=sqrt(eps()), which=SR())
        println(history)
        println(diag(schur.R))

        [(n+ℓ,ℓ) => normalize!(applied(*, R, (schur.Q[:,n]*sign(schur.Q[1,n])))) for n = 1:nev]
    elseif mode == :symbolic
        r = axes(R,1)
        [(n,ℓ) => applied(*, R, R \ hydredwfn(n,ℓ,Z).(r))
         for n ∈ ℓ .+ (1:nev)]
    else
        throw(ArgumentError("Unknown mode $(mode)"))
    end
end

ext: jl


text:
# ---
# title: 1092. Shortest Common Supersequence
# id: problem1092
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/shortest-common-supersequence/description/>
# hidden: true
# ---
#
# Given two strings `str1` and `str2`, return the shortest string that has both
# `str1` and `str2` as subsequences.  If multiple answers exist, you may return
# any of them.
#
# _(A string S is a subsequence of string T if deleting some number of
# characters from T (possibly 0, and the characters are chosen _anywhere_ from
# T) results in the string S.)_
#
#
#
# **Example 1:**
#
#
#
#     Input: str1 = "abac", str2 = "cab"
#     Output: "cabac"
#     Explanation:
#     str1 = "abac" is a subsequence of "cabac" because we can delete the first "c".
#     str2 = "cab" is a subsequence of "cabac" because we can delete the last "ac".
#     The answer provided is the shortest such string that satisfies these properties.
#
#
#
#
# **Note:**
#
#   1. `1 <= str1.length, str2.length <= 1000`
#   2. `str1` and `str2` consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Distributions, Parameters, Random
import Distributions: pdf,logpdf,rand
export LBA,pdf,logpdf,rand

#  See discussions at https://discourse.julialang.org/t/dynamichmc-reached-maximum-number-of-iterations/24721

############################################
#                              Model functions
############################################

mutable struct LBA{T1,T2,T3,T4} <: ContinuousUnivariateDistribution
    ν::T1
    A::T2
    k::T3
    τ::T4
    σ::Float64
end

Base.broadcastable(x::LBA)=Ref(x)

LBA(;τ,A,k,ν,σ=1.0) = LBA(ν,A,k,τ,σ)

function selectWinner(dt)
    if any(x->x >0,dt)
        mi,mv = 0,Inf
        for (i,t) in enumerate(dt)
            if (t > 0) && (t < mv)
                mi = i
                mv = t
            end
        end
    else
        return 1,-1.0
    end
    return mi,mv
end

function sampleDriftRates(ν,σ)
    noPositive=true
    v = similar(ν)
    while noPositive
        v = [rand(Normal(d,σ)) for d in ν]
        any(x->x>0,v) ? noPositive=false : nothing
    end
    return v
end

function rand(d::LBA)
    @unpack τ,A,k,ν,σ = d
    b=A+k
    N = length(ν)
    v = sampleDriftRates(ν,σ)
    a = rand(Uniform(0,A),N)
    dt = @. (b-a)/v
    choice,mn = selectWinner(dt)
    rt = τ .+ mn
    return choice,rt
end

function rand(d::LBA,N::Int)
    choice = fill(0,N)
    rt = fill(0.0,N)
    for i in 1:N
        choice[i],rt[i]=rand(d)
    end
    return (choice=choice,rt=rt)
end

logpdf(d::LBA,choice,rt) = log(pdf(d,choice,rt))

function logpdf(d::LBA,data::T) where {T<:NamedTuple}
    return sum(logpdf.(d,data...))
end

function logpdf(dist::LBA,data::Array{<:Tuple,1})
    LL = 0.0
    for d in data
        LL += logpdf(dist,d...)
    end
    return LL
end

function pdf(d::LBA,c,rt)
    @unpack τ,A,k,ν,σ = d
    b=A+k; den = 1.0
    rt < τ ? (return 1e-10) : nothing
    for (i,v) in enumerate(ν)
        if c == i
            den *= dens(d,v,rt)
        else
            den *= (1-cummulative(d,v,rt))
        end
    end
    pneg = pnegative(d)
    den = den/(1-pneg)
    den = max(den,1e-10)
    isnan(den) ? (return 0.0) : (return den)
end

function dens(d::LBA,v,rt)
    @unpack τ,A,k,ν,σ = d
    dt = rt-τ; b=A+k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    dens = (1/A)*(-v*cdf(Normal(0,1),n1) + σ*pdf(Normal(0,1),n1) +
        v*cdf(Normal(0,1),n2) - σ*pdf(Normal(0,1),n2))
    return dens
end

function cummulative(d::LBA,v,rt)
    @unpack τ,A,k,ν,σ = d
    dt = rt-τ; b=A+k
    n1 = (b-A-dt*v)/(dt*σ)
    n2 = (b-dt*v)/(dt*σ)
    cm = 1 + ((b-A-dt*v)/A)*cdf(Normal(0,1),n1) -
        ((b-dt*v)/A)*cdf(Normal(0,1),n2) + ((dt*σ)/A)*pdf(Normal(0,1),n1) -
        ((dt*σ)/A)*pdf(Normal(0,1),n2)
    return cm
end

function pnegative(d::LBA)
    @unpack ν,σ=d
    p=1.0
    for v in ν
        p*= cdf(Normal(0,1),-v/σ)
    end
    return p
end

function simulateLBA(;Nd=10,v=[1.0,1.5,2.0],A=.8,k=.2,tau=.4,kwargs...)
    return (rand(LBA(ν=v,A=A,k=k,τ=tau),Nd)...,N=Nd,Nc=length(v))
end


text:
# As described in Algorith, 7.3.4 in [CGTBOOK]
struct CGT <: TRSPSolver
end


text:
"""
    warp(img, tform, [indices], [degree = Linear()], [fill = NaN]) -> imgw

Transform the coordinates of `img`, returning a new `imgw`
satisfying `imgw[I] = img[tform(I)]`. This approach is known as
backward mode warping. The transformation `tform` must accept a
`SVector` as input. A useful package to create a wide variety of
such transformations is
[CoordinateTransformations.jl](https://github.com/FugroRoames/CoordinateTransformations.jl).

# Reconstruction scheme

During warping, values for `img` must be reconstructed at
arbitrary locations `tform(I)` which do not lie on to the lattice
of pixels. How this reconstruction is done depends on the type of
`img` and the optional parameter `degree`.

When `img` is a plain array, then on-grid b-spline interpolation
will be used. It is possible to configure what degree of b-spline
to use with the parameter `degree`. For example one can use
`degree = Linear()` for linear interpolation, `degree =
Constant()` for nearest neighbor interpolation, or `degree =
Quadratic(Flat())` for quadratic interpolation.

In the case `tform(I)` maps to indices outside the original
`img`, those locations are set to a value `fill` (which defaults
to `NaN` if the element type supports it, and `0` otherwise). The
parameter `fill` also accepts extrapolation schemes, such as
`Flat()`, `Periodic()` or `Reflect()`.

For more control over the reconstruction scheme --- and how
beyond-the-edge points are handled --- pass `img` as an
`AbstractInterpolation` or `AbstractExtrapolation` from
[Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl).

# The meaning of the coordinates

The output array `imgw` has indices that would result from
applying `inv(tform)` to the indices of `img`. This can be very
handy for keeping track of how pixels in `imgw` line up with
pixels in `img`.

If you just want a plain array, you can "strip" the custom
indices with `parent(imgw)`.

# Examples: a 2d rotation (see JuliaImages documentation for pictures)

```
julia> using Images, CoordinateTransformations, TestImages, OffsetArrays

julia> img = testimage("lighthouse");

julia> axes(img)
(Base.OneTo(512),Base.OneTo(768))

# Rotate around the center of `img`
julia> tfm = recenter(RotMatrix(-pi/4), center(img))
AffineMap([0.707107 0.707107; -0.707107 0.707107], [-196.755,293.99])

julia> imgw = warp(img, tfm);

julia> axes(imgw)
(-196:709,-68:837)

# Alternatively, specify the origin in the image itself
julia> img0 = OffsetArray(img, -30:481, -384:383);  # origin near top of image

julia> rot = LinearMap(RotMatrix(-pi/4))
LinearMap([0.707107 -0.707107; 0.707107 0.707107])

julia> imgw = warp(img0, rot);

julia> axes(imgw)
(-293:612,-293:611)

julia> imgr = parent(imgw);

julia> axes(imgr)
(Base.OneTo(906),Base.OneTo(905))
```
"""
function warp(img::AbstractExtrapolation{T}, tform, inds::Tuple = autorange(img, inv(tform))) where T
    out = similar(Array{T}, inds)
    warp!(out, img, try_static(tform, img))
end

function warp!(out, img::AbstractExtrapolation, tform)
    @inbounds for I in CartesianIndices(axes(out))
        out[I] = _getindex(img, tform(SVector(I.I)))
    end
    out
end

function warp(img::AbstractArray, tform, inds::Tuple, args...)
    etp = box_extrapolation(img, args...)
    warp(etp, try_static(tform, img), inds)
end

function warp(img::AbstractArray, tform, args...)
    etp = box_extrapolation(img, args...)
    warp(etp, try_static(tform, img))
end

"""
    imrotate(img, θ, [indices], [degree = Linear()], [fill = NaN]) -> imgr

Rotate image `img` by `θ`∈[0,2π) in a clockwise direction around its center point. To rotate the image counterclockwise, specify a negative value for angle.

By default, rotated image `imgr` will not be cropped. Bilinear interpolation will be used and values outside the image are filled with `NaN` if possible, otherwise with `0`.

# Examples
```julia
julia> img = testimage("cameraman")

# rotate with bilinear interpolation but without cropping
julia> imrotate(img, π/4)

# rotate with bilinear interpolation and with cropping
julia> imrotate(img, π/4)[axes(img)...]

# rotate with nearest interpolation but without cropping
julia> imrotate(img, π/4, Constant())
```

See also [`warp`](@ref), [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl).
"""
function imrotate(img::AbstractArray{T}, θ::Real, args...) where T
    θ = floor(mod(θ,2pi)*typemax(Int16))/typemax(Int16) # periodic discretezation
    tform = recenter(RotMatrix{2}(θ), center(img))
    warp(img, tform, args...)
end


text:
@testset "Sorts" begin
    sorts = [
        bubble_sort!
        bucket_sort!
        counting_sort!
        exchange_sort!
        insertion_sort!
        merge_sort!
        quick_sort!
        selection_sort!
    ]

    for f in sorts
        x = [3,5,1,4,2]
        f(x)
        @test x == [1,2,3,4,5]
    end
end

max_stars_repo_name: frankschmitt/Julia

max_issues_repo_name: frankschmitt/Julia

max_forks_repo_name: frankschmitt/Julia


text:
#-------- Inverses of DPR1

function inv{T}(A::SymDPR1{T},i::Integer,tols::Vector{Float64})

    # COMPUTES: inverse of a shifted SymDPR1 matrix A=diagm(A.D)+A.r*A.u*A.u',
    # inv(A-A.D[i]*I) which is a SymArrow.
    # Uses higher precision to compute top of the arrow element accurately, if
    # needed.
    # tols=[tolb,tolz] are tolerances, usually [1e3, 10*n]
    # [0.0,0.0] forces DoubleDouble, [1e50,1e50] would never use it
    # RETURNS:  SymArrow(D,z,b,i), Kb, Kz, Qout
    # Kb - condition Kb, Kz - condition Kz, Qout = 1 / 0 - double was / was not used

    n=length(A.D)
    D=Array(T,n-1)
    z=Array(T,n-1)
    wz=one(T)/A.u[i]
    shift=A.D[i]

    for k=1:i-1
        D[k]=one(T)/(A.D[k]-shift)
        z[k]=-A.u[k]*D[k]*wz
    end
    for k=i+1:n
        D[k-1]=one(T)/(A.D[k]-shift)
        z[k-1]=-A.u[k]*D[k-1]*wz
    end

    # compute the sum in a plain loop
    P=zero(T)
    Q=zero(T)
    Qout=0

    for k=1:i-1
        D[k]>0.0 ? P+=A.u[k]^2*D[k] : Q+=A.u[k]^2*D[k]
    end
    for k=i:n-1
        D[k]>0.0 ? P+=A.u[k+1]^2*D[k] : Q+=A.u[k+1]^2*D[k]
    end
    A.r>0 ? P=P+one(T)/A.r : Q=Q+one(T)/A.r

    Kb=(P-Q)/abs(P+Q)

    # compute Kz

    Kz=(sumabs(A.u)-abs(A.u[i]))*abs(wz)

    # Kz=maxabs(A.z)*abs(wz)

    if Kb<tols[1] ||  Kz<tols[2]
        b=(P+Q)*wz*wz
    else  # recompute in Double
        Qout=1
        shiftd=map(Double,A.D[i])
        Dd=[Double{Float64}[Double(A.D[k])-shiftd for k=1:i-1];
            Double{Float64}[Double(A.D[k])-shiftd for
                            k=i+1:length(A.D)]]
        wzd=Double(A.u[i])

        Pd,Qd=map(Double,(0.0,0.0))

        for k=1:i-1
            Dd[k].hi>0.0 ? Pd+=Double(A.u[k])^2/Dd[k] :
            Qd+=Double(A.u[k])^2/Dd[k]
        end

        for k=i+1:n
            Dd[k-1].hi>0.0 ? Pd+=Double(A.u[k])^2/Dd[k-1] :
            Qd+=Double(A.u[k])^2/Dd[k-1]
            # @show P,Q
        end

        A.r > 0 ?   Pd+=Double(1.0)/Double(A.r)  : Qd+=Double(1.0)/Double(A.r)

        bd=(Pd+Qd)/(wzd*wzd)
        b=bd.hi+bd.lo
    end

    # return this
    SymArrow(D,z,b,i),Kb,Kz,Qout

end # inv

function inv{T}(A::SymDPR1{T}, shift::Float64, tolr::Float64)

    # COMPUTES: inverse of the shifted SymDPR1 A = diagm(A.D)+A.r*A.u*A.u',
    # inv(A-shift*I) = D + rho*u*u', shift!=A.D[i], which is again a SymDPR1
    # uses DoubleDouble to compute A.r accurately, if needed.
    # tolr is tolerance, usually 1e3,  0.0 forces Double, 1e50 would never use it
    # RETURNS: SymDPR1(D,u,rho), Krho, Qout
    # Krho - condition Krho, Qout = 1 / 0 - Double was / was not used

    n=length(A.D)
    D=Array(T,n)
    u=Array(T,n)

    for k=1:n
        D[k]=one(T)/(A.D[k]-shift)
        u[k]=A.u[k]*D[k]
    end

    # compute gamma and Kgamma
    #--- compute the sum in a plain loop

    P=zero(T)
    Q=zero(T)
    Qout=0

    for k=1:n
        D[k]>0.0 ? P=P+A.u[k]^2*D[k] : Q=Q+A.u[k]^2*D[k]
    end

    A.r>0 ? P=P+one(T)/A.r : Q=Q+one(T)/A.r

    # Condition of rho
    Krho=(P-Q)/abs(P+Q)

    if Krho < tolr
        rho=-one(T)/(P+Q)

    else  # recompute in Double
        Qout=1
        Pd,Qd=map(Double,(0.0,0.0))
        shiftd=Double(shift)

        for k=1:n
            D[k]>0.0 ? Pd=Pd+Double(A.u[k])^2/(Double(A.D[k])-shiftd) : Qd=Qd+Double(A.u[k])^2/(Double(A.D[k])-shiftd)
        end

        A.r > 0 ?   Pd+=Double(1.0)/Double(A.r)  : Qd+=Double(1.0)/Double(A.r)
        r=Double(1.0)/(Pd+Qd)
        rho=-(r.hi+r.lo)

    end

    # returns the following
    SymDPR1(D,u,rho), Krho, Qout

end # inv

function inv{T}(A::SymDPR1{T}, shift::Double)

    # COMPUTES: inverse of the shifted SymDPR1 A = diagm(A.D)+A.r*A.u*A.u',
    # inv(A-shift*I) = D + rho*u*u', shift!=A.D[i], which is again a SymDPR1
    # here shift is Double so it uses Double to compute everything
    # RETURNS: SymDPR1(D,u,rho), Qout
    # Qout = 1 on exit meaning Double was used

    n=length(A.D)
    D=Array(Double,n)
    u=Array(Double,n)

    for k=1:n
        D[k]=Double(A.D[k])-shift
    end

    u=map(Double,A.u)

    oned=Double(1.0)
    zerod=Double(0.0)

    for k=1:n
        u[k]=u[k]/D[k]
        D[k]=oned/D[k]
    end

    # compute rho
    # compute the sum in a plain loop

    P,Q=zerod,zerod
    Qout=1

    for k=1:n
        D[k].hi > 0.0 ? P=P+Double(A.u[k])*u[k] : Q=Q+Double(A.u[k])*u[k]
    end

    A.r > 0 ?   P+=Double(1.0)/Double(A.r)  : Q+=Double(1.0)/Double(A.r)

    r=oned/(P+Q)
    rho=-(r.hi+r.lo)

    # returns the following
    # SymDPR1(T[x.hi+x.lo for x=D],T[x.hi+x.lo for x=u],rho), Qout

    D1=Array(T,n)
    u1=Array(T,n)

    for k=1:n
        D1[k]=D[k].hi+D[k].lo
        u1[k]=u[k].hi+u[k].lo
    end

    SymDPR1(D1,u1,rho), Qout

end # inv

function  eig( A::SymDPR1,k::Integer,tols::Vector{Float64})

    # COMPUTES: k-th eigenpair of an ordered irreducible SymDPR1
    # A = diagm(A.D)+A.r*A.u*A.u', A.r > 0
    # tols=[tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3]
    # RETURNS: lambda, v, Sind, Kb, Kz, Knu, Krho, Qout
    # lambda - k-th eigenvalue in descending order
    # v - lambda's normalized eigenvector
    # Kb, Kz, Knu, Krho - condition numbers
    # Qout = 1 / 0 - Double was / was not used

    # Set the dimension
    n = length(A.D)

    # Set all conditions initially to zero
    Kb,Kz,Knu,Krho=0.0,0.0,0.0,0.0
    Qout=0
    v=zeros(n)
    # Kz is former kappa_nu

    # Determine the shift sigma, the shift index i, and whether lambda
    # is on the left or the right side of the nearest pole
    # Exterior eigenvalues (k = 1 or k = n):

    if k == 1
        sigma,i,side = A.D[1],1,'R'
    else
        # Interior eigenvalues (k in (2,...n-1) ):
        Dtemp = A.D-A.D[k]
        middle = Dtemp[k-1]/2.0
        Fmiddle = 1.0+A.r*sum(A.u.^2./(Dtemp-middle))
        sigma,i,side = Fmiddle > 0.0 ? (A.D[k],k,'R') : (A.D[k-1],k-1,'L')
    end

    # Compute the inverse of the shifted matrix, A_i^(-1), Kb and Kz

    Ainv,Kb,Kz,Qout = inv(A,i,tols[1:2])

    # Compute the eigenvalue of the inverse shifted matrix

    nu = bisect( Ainv,side )

    #  nu=fastzero([invD1; 0; invD2], [w1;wz;w2], b, side); # implement later
    #  [nu-nu1]/nu, [nu-nueva]/nu, pause  # and compare

    if abs(nu)==Inf
        # this is nonstandard
        # Deflation in dpr1eig (nu=Inf)
        v[i]=1.0
        lambda=sigma
    else
        # standard case, full computation
        # nu1 is the F- or 1-norm of the inverse of the shifted matrix
        # nu10=maximum([sum(abs(Ainv.z))+abs(Ainv.a), maximum(abs(Ainv.D)+abs(Ainv.z))])
        nu1=0.0
        for k=1:n-1
            nu1=max(nu1,abs(Ainv.D[k])+abs(Ainv.z[k]))
        end
        nu1=max(sumabs(Ainv.z)+abs(Ainv.a), nu1)

        Knu=nu1/abs(nu)
        if Knu<tols[3]
            # Accuracy is fine, compute the eigenvector
            mu = 1.0/nu
            # v=[ A.u./((A.D-sigma)-mu)]
            for l=1:n
                v[l]= A.u[l]/((A.D[l]-sigma)-mu)
            end
            #        for k=1:n-1
            #	   v[k] = A.z[k]/((A.D[k]-sigma)-mu)
            #        end
            #        v[n]=-one
            lambda, v = mu+sigma, v/norm(v)
        else
            # Remedies according to Remark 3 - we shift between original
            # eigenvalues and compute DPR1 matrix
            # 1/nu1+sigma, 1/nu+sigma
            println("Remedy 3 ")
            nu = side=='R' ? abs(nu) : -abs(nu)
            nu1=-sign(nu)*nu1
            sigma1=(nu1+nu)/(2.0*nu*nu1)+sigma

            if findfirst(A.D-sigma1,0.0)>0 # we came back with a pole
                # recompute sigmav more accurately according with dekker
                sigmav=(Double(nu1)+Double(nu))/(Double(2.0)*Double(nu)*Double(nu1))+Double(sigma)
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv,Qout1=inv(A,sigmav) # Ainv is Float64
                nu1=bisect(Ainv,side)
                mu1 = 1.0/nu1
                D0=map(A.D,Double)-sigmav
                D1=D0.hi+D0.lo
                if findfirst(D1-mu1,0.0)>0
                    ind=find(D1-mu1==0.0);
                    v=zeros(n)
                    v[ind]=1.0;
                else
                    v=[ A.u./(D1-mu1)]
                end
                # Shift the eigenvalue back in Double
                lam = Double(1.0)/Double(nu1)+sigmav
                # Return this
                lambda,v = lam.hi+lam.lo, v/norm(v)
            else
                # Compute the inverse of the shifted arrowhead (DPR1)
                Ainv, Krho,Qout1=inv(A,sigma1,tols[4]) # Ainv is Float64
                # Compute the eigenvalue by bisect for DPR1
	        # Note: instead of bisection could use dlaed4 (need a wrapper) but
	        # it is not faster. There norm(u)==1
                nu1= bisect(Ainv,side)
                mu1=1.0/nu1
                # standard v
                v=A.u./((A.D-sigma1)-mu1)
                # Return this - shift the eigenvalue back and normalize the vector
                lambda, v = mu1+sigma1, v/norm(v)
            end
            Qout==1 && (Qout=Qout+2)
        end

        # Remedy according to Remark 1 - we recompute the the eigenvalue
        # near zero from the inverse of the original matrix (a DPR1 matrix).
        if (abs(A.D[i])+abs(1.0/nu))/abs(lambda)>tols[5]

            if (k==1 && A.D[1]<0.0 || side=='L' && sign(A.D[i])+sign(A.D[i+1])==0 ||
                side=='R' && sign(A.D[i])+sign(A.D[i-1])==0)
                println("Remedy 1 ")
                # Compute the inverse of the original arrowhead (DPR1)
                Ainv,Krho,Qout1 = inv(A,0.0,tols[4]) # Ainv is Float64
                Qout==1 && (Qout=Qout+4)
                if abs(Ainv.r)==Inf
                    lambda=0.0
                else
                    # Here we do not need bisection. We compute the Rayleigh
                    # quotient by using already computed vectors which is
                    # componentwise accurate
                    nu1=sum(v.^2.*Ainv.D)+Ainv.r*sum(v.*Ainv.u)^2;
                    lambda=1.0/nu1
                end
            end

        end
    end

    # Return this
    lambda,v,i,Kb,Kz,Knu,Krho,Qout

end # eig (k)

function eig(A::SymDPR1, tols::Vector{Float64})

    # COMPUTES: all eigenvalues and eigenvectors of a real symmetric SymDPR1
    # A = diagm(A.D)+A.r*A.u*A.u'
    # tols = [tolb,tolz,tolnu,tolrho,tollambda] = [1e3,10.0*n,1e3,1e3,1e3] or similar
    # RETURNS: U, E, Sind, Kb, Kz, Knu, Krho, Qout
    # U = eigenvectors, E = eigenvalues in decreasing order
    # Sind[k] - shift index i for the k-th eigenvalue
    # Kb, Kz, Knu, Krho [k] - respective conditions for the k-th eigenvalue
    # Qout[k] = 1 / 0 - Double was / was not used when computing k-th eigenvalue

    n=length(A.D)
    n0=n

    # Checking if A.r > 0
    signr =  A.r > 0.0 ? 1.0 : -1.0

    # Ordering the matrix
    D=signr*A.D
    is=sortperm(D,rev=true)
    D=D[is]
    z=A.u[is]
    rho=signr*A.r

    U=eye(n,n)
    E=zeros(n)
    Kb=zeros(n); Kz=zeros(n); Knu=zeros(n); Krho=zeros(n)
    Qout=zeros(Int,n); Sind=zeros(Int,n)

    # Quick return for 1x1, this is trivial for SymArrow, not so trivial here :)
    if n==1
        U=1;
        if (D==0)&&((rho==0)|(z==0))
            E=0
        else
            E=A.D[1]+A.r*A.u[1]^2
            # Higher accuracy if needed
            KD=(abs(A.D[1])+abs(A.r)*A.u[1]^2)/abs(E)
            if KD>tols[1]
                Ed=Double(A.D[1])+Double(A.r)*Double(A.u[1])^2
                E=Ed.hi+Ed.lo
            end
        end
        return U, [E], Sind,Kb,Kz,Knu,Krho,Qout
    end

    #  test for deflation in z
    z0=find(z.==0)
    zx=find(z.!=0)

    if isempty(zx)  # nothing to do
        E=A.D
        isE=sortperm(E,rev=true)
        E=E[isE]
        U=U[:,isE]
        return U,E,Sind,Kb,Kz,Knu,Krho,Qout
    end

    if !isempty(z0)
        E[z0]=D[z0]
        D=D[zx]
        z=z[zx]
        if !isempty(z)
            # return U,E,Sind,Kb,Kz,Knu,Krho,Qout
            n=length(z)
        end
    end

    #  Test for deflation in D
    g=D[1:n-1]-D[2:n]
    # Can play with inexact deflation
    # g0=find(abs(g)<eps)
    # gx=find(abs(g)>=eps)
    # Only exact deflation !!
    g0=find(g.==0.0)
    gx=find(g.!=0.0)
    if !isempty(g0)
        # Deflation
        Dgx=D[gx]; zgx=z[gx]
        lg0=length(g0)
        R=Array(Tuple{Givens{Float64},Float64},lg0)
        for l=lg0:-1:1
            R[l]=givens(z[g0[l]],z[g0[l]+1],zx[g0[l]],zx[g0[l]+1])
            z[g0[l]]=R[l][2]; z[g0[l]+1]=0.0
            # A_mul_Bc!(U,R) # multiply R'*U later
            E[zx[g0[l]+1]]=D[g0[l]+1]
        end
        # remains
        gx=[0;gx]+1
        nn=length(gx)

        zxx=zx[gx]
        for k=1:nn
            E[zxx[k]],U[zxx,zxx[k]],Sind[zxx[k]],Kb[zxx[k]],Kz[zxx[k]],Knu[zxx[k]],Krho[zxx[k]],Qout[zxx[k]]=
            eig(SymDPR1(D[gx],z[gx],rho),k,tols)
        end

        for l=1:lg0
            U=R[l][1]'*U
        end

    else

        # No deflation in D
        for k=1:n
            E[zx[k]],U[zx,zx[k]],Sind[zx[k]],Kb[zx[k]],Kz[zx[k]],Knu[zx[k]],Krho[zx[k]],Qout[zx[k]]=
            eig(SymDPR1(D,z,rho),k,tols)
        end
    end

    # back premutation of vectors
    isi=sortperm(is)

    # change the sign if A.r was negative
    # must sort E once more
    E=signr*E
    es=sortperm(E,rev=true)
    E=E[es]

    U=U[isi,es]

    # Return this
    U,E,Sind[es],Kb[es],Kz[es],Knu[es],Krho[es],Qout[es]

end # eig (all)


text:
module error

import builtin
import io
import struct

class Token {
  new(field text,
      field line,
      field col,
      field fn,
      field line_text) {}
  method to_s() {
    concat('Token("', text, '":',
           line, ',',
           col, ',',
           fn, ',',
           line_text, ')')
  }
}

class StackLine {
  field mod, caller, ip, token
  new(block) {
    mod = block.$module
    caller = block.$caller
    ip = block.$ip
    token = Token(builtin.token__(mod, ip))
  }
  method to_s() {
    tos = ''
    if caller {
      if caller is MethodInstance {
        tos.extend(caller.$method.method_path())
      } else if caller is Method {
        tos.extend(caller.method_path())
      } else {
        tos.extend(concat(mod.name,'.', caller.name))
      }
    } else {
      tos.extend('\t').extend(mod.name)
    }
    tos.extend(concat('(', token.fn, ':', token.line, ')'))
  }
}

class Error {
  field stack_lines
  new(msg) {
    self.msg = if (msg is String) then msg else if (~msg) 'None' else concat(msg)
    stack_lines = []
    ; Start at 1 since 0-th block will be this loop.
    for i=1, i<$thread.$saved_blocks.len, i=i+1 {
      block = $thread.$saved_blocks[i]
      stack_lines.append(StackLine(block))
    }
  }
  method to_s() {
    head = stack_lines[0]
    tok = head.token
    tos = concat('Error in ', tok.fn,
                 ' at line ', tok.line, ' col ', tok.col, ': ',
                 self.msg, '\n')
    tos.extend(concat(tok.line, ':', tok.line_text))
    if tok.line_text and ~str(tok.line_text).ends_with('\n') {
      tos.extend('\n')
    }
    for i=0, i<tok.col + 2, i=i+1 {
      tos.extend(' ')
    }
    for i=0, i<tok.text.len, i=i+1 {
      tos.extend('^')
    }
    tos.extend('\n')
    for i=0, i<stack_lines.len, i=i+1 {
      sl = stack_lines[i]
      tos.extend(sl.to_s()).extend('\n')
    }
    tos
  }
}

def display_error(e) {
  io.fprintln(io.ERROR, e)
}

def raise_error(e) {
  io.fprintln(io.ERROR, e)
  exit
}


text:
using LinearAlgebra
using SparseArrays
using Random
using ProximalOperators
using Test

@testset "IndAffine" begin

# Full matrix

m, n = 10, 30
A = randn(m, n)
b = randn(m)
f = IndAffine(A, b)
x = randn(n)

predicates_test(f)

@test ProximalOperators.is_smooth(f) == false
@test ProximalOperators.is_quadratic(f) == false
@test ProximalOperators.is_generalized_quadratic(f) == true
@test ProximalOperators.is_set(f) == true

call_test(f, x)
y, fy = prox_test(f, x)

@test f(y) == 0.0

# Sparse matrix

m, n = 10, 30
A = sprandn(m, n, 0.5)
b = randn(m)
f = IndAffine(A, b)
x = randn(n)

call_test(f, x)
y, fy = prox_test(f, x)

@test f(y) == 0.0

# Iterative version

m, n = 200, 500
A = sprandn(m, n, 0.5)
b = randn(m)
f = IndAffine(A, b; iterative=true)
x = randn(n)

call_test(f, x)
y, fy = prox_test(f, x)

@test f(y) == 0.0

end


text:
using SymArrays
using Documenter

DocMeta.setdocmeta!(SymArrays, :DocTestSetup, :(using SymArrays); recursive=true)

makedocs(;
    modules=[SymArrays],
    authors="Johannes Feist <johannes.feist@gmail.com> and contributors",
    repo="https://github.com/jfeist/SymArrays.jl/blob/{commit}{path}#{line}",
    sitename="SymArrays.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jfeist.github.io/SymArrays.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jfeist/SymArrays.jl",
    devbranch = "main",
)


text:
abstract type Game{H,K} end

infokeytype(::Game{H,K}) where {H,K} = K
histtype(::Game{H,K}) where {H,K} = H

"""
    `initialhist(game::Game)`

Return initial history with which to start the game
"""
function initialhist end

"""
    `isterminal(game::Game, h)`

Returns boolean - whether or not current history is terminal \n
i.e h ∈ Z
"""
function isterminal end

"""
    `utility(game::Game, i::Int, h)`

Returns utility of some history h for some player i
"""
function utility end

"""
    `player(game::Game{H,K}, h::H)`

Returns integer id corresponding to which player's turn it is at history h
0 - Chance Player
1 - Player 1
2 - Player 2
\n
If converting to IIE to Matrix Game need to implement:
    `player(game::Game{H,K}, k::K)`
"""
function player end

"""
    `chance_action(game::Game, h)`

Return randomly sampled action from chance player at a given history
"""
function chance_action end

"""
    chance_actions(game::Game, h)

Return all chance actions available for chance player at history h
"""
function chance_actions end

"""
    `next_hist(game::Game, h, a)`

Given some history and action return the next history
`h′ = next_hist(game, h, a)`
"""
function next_hist end

"""
    `infokey(game::Game, h)`

Returns unique identifier corresponding to some information set \n
`infokey(game, h1) == infokey(game, h2)` ⟺ h1 and h2 belong to the same info set \n
(key must be immutable as it's being stored as a key in a dictionary)
"""
function infokey end

"""
    `actions(game::Game, h)`

Returns all actions available at some history
"""
function actions end

@inline other_player(i) = 3-i

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
include("gs.jl")

using .GS

nos = [No(1, 0, 0; rx=true, ry=true),
       No(2, 1, 0),
       No(3, 2, 0; fy=-10),
       No(4, 0, 1; rx=true, ry=true),
       No(5, 1, 1),
       No(6, 2, 1)]

material = Material(1, 7e7)
vol = 0.0003
Lt = 12.65685424949238
area = 1 #vol / Lt

els = [1 2
       1 4
       1 5
       2 3
       2 4
       2 5
       2 6
       3 5
       3 6
       4 5
       5 6]

elementos = [Elemento(i, nos[els[i, 1]], nos[els[i, 2]], area, material) for i=1:size(els)[1]]

estrutura = Estrutura(nos, elementos)

otm = Otimizador(vol, estrutura)
otimizar!(otm)
plotar_estrutura_otimizada(otm)
#plotar_estrutura_inicial(estrutura)

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


text:
import Base: BitSet, convert, reverse, print_array, sum, promote_rule
import LinearAlgebra: normalize, BLAS.dotu

const IntCol = Union{AbstractVector{Int},BitSet,Set{Int},NTuple{N,Int} where N}

convert(::Type{BitSet}, r::IntCol) = BitSet(r)

BitSet(elems::Vararg{Int}) = BitSet(elems)

reverse(n::Number) = n

normalize(z::Complex) = z == 0 ? one(z)/abs(one(z)) : z/abs(z)

normalize(z::Real) = one(z)

trace(A::AbstractMatrix) = tr(A)

sum(A::AbstractArray,i::Integer) = Base._sum(A,i)

parseb2(s::AbstractString) = Base.tryparse_internal(Int,s,firstindex(s),lastindex(s),2,true)

const MatrixWrapper{T,S} = Union{Adjoint{T,S}, Transpose{T,S}}
const AllAbstractSparseMatrix{T} = Union{
    AbstractSparseMatrix{T},
    MatrixWrapper{T,<:AbstractSparseMatrix{T}},
    MatrixWrapper{T,<:MatrixWrapper{T,<:AbstractSparseMatrix{T}}}}

promote_rule(::Type{<:AllAbstractSparseMatrix{T1}},::Type{<:AllAbstractSparseMatrix{T2}}) where {T1,T2} = SparseMatrixCSC{promote_type(T1,T2),Int}

promote_rule(::Type{<:AllAbstractSparseMatrix{T1}},::Type{<:AbstractMatrix{T2}}) where {T1,T2} = Matrix{promote_type(T1,T2)}
promote_rule(::Type{<:AbstractMatrix{T1}},::Type{<:AllAbstractSparseMatrix{T2}}) where {T1,T2} = Matrix{promote_type(T1,T2)}

promote_rule(::Type{<:AllAbstractSparseMatrix{T1}},::Type{Matrix{T2}}) where {T1,T2} = Matrix{promote_type(T1,T2)}
promote_rule(::Type{Matrix{T1}},::Type{<:AllAbstractSparseMatrix{T2}}) where {T1,T2} = Matrix{promote_type(T1,T2)}

promote_rule(::Type{Matrix{T1}},::Type{<:MatrixWrapper{T2,<:Matrix{T2}}}) where {T1,T2} = Matrix{promote_type(T1,T2)}
promote_rule(::Type{<:MatrixWrapper{T1,<:Matrix{T1}}},::Type{Matrix{T2}}) where {T1,T2} = Matrix{promote_type(T1,T2)}

promote_rule(::Type{SparseVector{T1,S1}},::Type{SparseVector{T2,S2}}) where {T1,T2,S1,S2} = SparseVector{promote_type(T1,T2),promote_type(S1,S2)}

dotu(x::AbstractVector{T1},y::AbstractVector{T2}) where {T1<:Real,T2<:Real} = dot(x,y)


text:
"""

Automatically generated from SageMath calculations

Fergus Baker - 9th Nov 2021
             - 10th Feb 2022: updated to include Jacobian method
             - 1st Apr 2022: fix sign error in mu
"""
module EddingtonFinkelsteinCoords

using ..ComputedGeodesicEquations

@inline function geodesic_eq(u, v, M)
    @fastmath ComputedGeodesicEquations.@let_unpack u v begin
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        out1 =
            2 * (
                M * r^3 * v_phi^2 * sin_theta^2 + M * r^3 * v_theta^2 - M^2 * v_t^2 -
                (M^2 + M * r) * v_r^2 - (2 * M^2 + M * r) * v_r * v_t
            ) / r^3
        out2 =
            -(
                (2 * M * r^3 - r^4) * v_phi^2 * sin_theta^2 - 4 * M^2 * v_r * v_t -
                (2 * M^2 + M * r) * v_r^2 - (2 * M^2 - M * r) * v_t^2 +
                (2 * M * r^3 - r^4) * v_theta^2
            ) / r^3
        out3 = (r * v_phi^2 * cos_theta * sin_theta - 2 * v_r * v_theta) / r
        out4 =
            -2 * (r * v_phi * v_theta * cos_theta + v_phi * v_r * sin_theta) /
            (r * sin_theta)

        (out1, out2, out3, out4)
    end
end

@inline function constrain(μ, u, v, M)
    @fastmath ComputedGeodesicEquations.@let_unpack u v begin
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        -(
            2 * M * v_r - sqrt(
                -(2 * M * r^3 - r^4) * v_phi^2 * sin_theta^2 - 2 * M * μ^2 * r +
                μ^2 * r^2 +
                r^2 * v_r^2 - (2 * M * r^3 - r^4) * v_theta^2,
            )
        ) / (2 * M - r)
    end
end

@inline function jacobian(u, M)
    @fastmath let t = u[1], r = u[2], theta = u[3], phi = u[4]
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        comp1 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 0 0
            0 0 0 0
        ]
        comp2 = ComputedGeodesicEquations.@SMatrix [
            -2*M/r^2 -2*M/r^2 0 0
            -2*M/r^2 -2*M/r^2 0 0
            0 0 2*r 0
            0 0 0 2*r*sin_theta^2
        ]
        comp3 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 0 0
            0 0 0 2*r^2*cos_theta*sin_theta
        ]
        comp4 = ComputedGeodesicEquations.@SMatrix [
            0 0 0 0
            0 0 0 0
            0 0 0 0
            0 0 0 0
        ]
        (comp1, comp2, comp3, comp4)
    end
end

@inline function metric(u, M)
    @fastmath let t = u[1], r = u[2], theta = u[3], phi = u[4]
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        ComputedGeodesicEquations.@SMatrix [
            2*M/r-1 2*M/r 0 0
            2*M/r 2*M/r+1 0 0
            0 0 r^2 0
            0 0 0 r^2*sin_theta^2
        ]
    end
end

@inline function inverse_metric(u, M)
    @fastmath let t = u[1], r = u[2], theta = u[3], phi = u[4]
        cos_theta = cos(theta)
        sin_theta = sin(theta)

        ComputedGeodesicEquations.@SMatrix [
            -(2 * M + r)/r 2*M/r 0 0
            2*M/r -(2 * M - r)/r 0 0
            0 0 r^(-2) 0
            0 0 0 1/(r^2*sin_theta^2)
        ]
    end
end

end # module

@with_kw struct EddingtonFinkelstein{T} <: AbstractMetricParams{T}
    @deftype T
    M = 1.0
end

@with_kw struct EddingtonFinkelsteinJac{T} <: AbstractMetricParams{T}
    @deftype T
    M = 1.0
end

geodesic_eq(m::EddingtonFinkelstein{T}, u, v) where {T} =
    EddingtonFinkelsteinCoords.geodesic_eq(u, v, m.M)
geodesic_eq(m::EddingtonFinkelsteinJac{T}, u, v) where {T} = jac_geodesic_eq(m, u, v)

constrain(m::EddingtonFinkelstein{T}, u, v; μ::T = 0.0) where {T} =
    EddingtonFinkelsteinCoords.constrain(μ, u, v, m.M)
constrain(m::EddingtonFinkelsteinJac{T}, u, v; μ::T = 0.0) where {T} =
    EddingtonFinkelsteinCoords.constrain(μ, u, v, m.M)

# specialisations
metric(m::EddingtonFinkelstein{T}, u) where {T} = EddingtonFinkelsteinCoords.metric(u, m.M)
inverse_metric(m::EddingtonFinkelstein{T}, u) where {T} =
    EddingtonFinkelsteinCoords.inverse_metric(u, m.M)
jacobian(m::EddingtonFinkelstein{T}, u) where {T} =
    EddingtonFinkelsteinCoords.jacobian(u, m.M)
metric(m::EddingtonFinkelsteinJac{T}, u) where {T} =
    EddingtonFinkelsteinCoords.metric(u, m.M)
inverse_metric(m::EddingtonFinkelsteinJac{T}, u) where {T} =
    EddingtonFinkelsteinCoords.inverse_metric(u, m.M)
jacobian(m::EddingtonFinkelsteinJac{T}, u) where {T} =
    EddingtonFinkelsteinCoords.jacobian(u, m.M)

export EddingtonFinkelsteinCoords, EddingtonFinkelstein, EddingtonFinkelsteinJac

# additional specializations
inner_radius(m::EddingtonFinkelstein{T}) where {T} = 2 * m.M
inner_radius(m::EddingtonFinkelsteinJac{T}) where {T} = 2 * m.M


text:
## plot.jl : functions for easy visualization of Gaussian random fields

# 1D
plot(grf::GaussianRandomField{G, C}; kwargs...) where {G, C <: AbstractCovarianceFunction{1}} = plot(grf.pts[1], sample(grf); kwargs...)

plot!(grf::GaussianRandomField{G, C}; kwargs...) where {G, C <: AbstractCovarianceFunction{1}} = plot!(grf.pts[1], sample(grf); kwargs...)

# 2D
for f in [:surface, :contour, :contourf, :heatmap]
    expr = quote
        function $f(grf::GaussianRandomField{G, C}; kwargs...) where {G, C <: AbstractCovarianceFunction{2}}
            x, y = grf.pts
            $f(x, y, sample(grf); kwargs...)
        end
    end
    eval(expr)
end

# Plot eigenvalues
plot_eigenvalues_sub(eigenvalues) = plot(1:length(eigenvalues), eigenvalues, xaxis=:log, yaxis=:log, xlabel="eigenvalue", ylabel="magnitude")

plot_eigenvalues(grf::GaussianRandomField{<:KarhunenLoeve}) = plot_eigenvalues_sub(grf.data.eigenval.^2)

plot_eigenvalues(grf::GaussianRandomField{<:CirculantEmbedding}) = begin
    data = grf.data[1]
    eigenvalues = sort(vec(data), rev=true).^2 .* length(data)
    plot_eigenvalues_sub(eigenvalues)
end

plot_eigenvalues(grf::GaussianRandomField{<:KarhunenLoeve,<:SeparableCovarianceFunction}) = begin
    order, data = grf.data
    eigenvalues = Array{Float64}(undef, length(order))
    @inbounds for (i, o) in enumerate(order)
        eigenvalues[i] = prod(data[j].eigenval[o[j]] for j in 1:length(data))^2
    end
    plot_eigenvalues_sub(eigenvalues)
end

# Plot eigenfunctions
plot_eigenfunction_1d(x, eigenfunction; kwargs...) = plot(x, eigenfunction; kwargs...)

plot_eigenfunction_2d(x, y, eigenfunction; kwargs...) = contourf(x, y, eigenfunction; kwargs...)

plot_eigenfunction(grf::GaussianRandomField{<:KarhunenLoeve, C}, n; kwargs...) where C <: AbstractCovarianceFunction{1} = begin
    x = grf.pts[1]
    plot_eigenfunction_1d(x, view(grf.data.eigenfunc, :, n); kwargs...)
end

plot_eigenfunction(grf::GaussianRandomField{<:KarhunenLoeve, C}, n) where C <: SeparableCovarianceFunction{1} = begin
    x = grf.pts[1]
    order, data = grf.data
    ordern = order[n]
    eigenfunction = view(data[1].eigenfunc, :, ordern[1])
    plot_eigenfunction_1d(x, eigenfunction)
end

plot_eigenfunction(grf::GaussianRandomField{<:KarhunenLoeve, C}, n; kwargs...) where C <: AbstractCovarianceFunction{2} = begin
    x, y = grf.pts
    plot_eigenfunction_2d(x, y, view(grf.data.eigenfunc, :, n); kwargs...)
end

plot_eigenfunction(grf::GaussianRandomField{<:KarhunenLoeve, C}, n) where C <: SeparableCovarianceFunction{2} = begin
    x, y = grf.pts
    order, data = grf.data
    ordern = order[n]
    eigenfunction = kron((view(data[j].eigenfunc, :, ordern[j]) for j = 1:length(grf.cov.cov))...)
    plot_eigenfunction_2d(x, y, eigenfunction)
end

# Plot covariance matrix
function plot_covariance_matrix(c::CovarianceFunction{d},
                                pts::Vararg{AbstractRange, d}) where d
    heatmap(apply(c, pts, pts))
end


text:
"""
Function to obtain the output current time series of a Dynamic Inverter model out of the DAE Solution. It receives the simulation inputs,
the dynamic device and bus voltage. It is dispatched for device type to compute the specific current.

"""
function compute_output_current(
    sim::Simulation,
    dynamic_device::G,
    V_R::Vector{Float64},
    V_I::Vector{Float64},
) where {G <: PSY.DynamicInverter}

    #Obtain Data
    sys = get_system(sim)

    #Get machine
    filt = PSY.get_filter(dynamic_device)
    Sbase = PSY.get_base_power(sys)
    basepower = PSY.get_base_power(dynamic_device)
    base_power_ratio = basepower / Sbase
    return _filter_current(
        filt,
        PSY.get_name(dynamic_device),
        V_R,
        V_I,
        base_power_ratio,
        sim,
    )
end

"""
Function to obtain the output current time series of a LCL Filter model out of the DAE Solution. It is dispatched via the Filter type.

"""
function _filter_current(
    ::PSY.LCLFilter,
    name::String,
    ::Vector{Float64},
    ::Vector{Float64},
    base_power_ratio::Float64,
    sim::Simulation,
)
    ir_filter = post_proc_state_series(sim, (name, :ir_filter))
    ii_filter = post_proc_state_series(sim, (name, :ii_filter))

    return base_power_ratio * ir_filter, base_power_ratio * ii_filter
end


text:
mutable struct TaggedIntrusiveSList{T, Tag}
    tail::Union{T, Nothing}

    TaggedIntrusiveSList{T, Tag}() where {T, Tag} = new(nothing)
end

const IntrusiveSList{T} = TaggedIntrusiveSList{T, :default}

@inline function checkbounds(list::TaggedIntrusiveSList)
    isempty(list) && throw(BoundsError(list))
end

Base.eltype(::Type{TaggedIntrusiveSList{T, Tag}}) where {T, Tag} = T

@inline Base.isempty(list::TaggedIntrusiveSList) = isnothing(list.tail)
@inline Base.empty!(list::TaggedIntrusiveSList) = list.tail = nothing

# TODO: this is duplicate of list code
function Base.length(list::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    c = 0
    for n in list
        c += 1
    end
    c
end

function Base.show(io::IO, list::TaggedIntrusiveSList)
    first = true
    print(io, "[")
    for x in list
        !first && print(io, "→")
        first = false
        print(io, x)
    end
    print(io, "]")
end

@inline function Base.iterate(list::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    isempty(list) && return nothing
    tail = list.tail::T
    head = getnext(tail, Val{Tag}())::T
    (head, head)
end

@inline function Base.iterate(list::TaggedIntrusiveSList{T, Tag}, state::T) where {T, Tag}
    if state === list.tail
        nothing
    else
        x = getnext(state, Val{Tag}())::T
        (x, x)
    end
end

@inline function Base.first(list::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    @boundscheck checkbounds(list)
    getnext(list.tail::T, Val{Tag}())::T
end

@inline function Base.last(list::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    @boundscheck checkbounds(list)
    list.tail::T
end

function Base.push!(list::TaggedIntrusiveSList{T, Tag}, node::T) where {T, Tag}
    if isempty(list)
        setnext!(node, node, Val{Tag}())
        list.tail = node
    else
        tail = list.tail::T
        head = getnext(tail, Val{Tag}())::T
        setnext!(tail, node, Val{Tag}())
        setnext!(node, head, Val{Tag}())
        list.tail = node
    end
    list
end

function Base.pop!(list::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    @boundscheck checkbounds(list)
    tail = list.tail::T
    head = getnext(tail, Val{Tag}())::T
    if head === tail
        list.tail = nothing
    else
        prev = circslist_prev(tail, head, Val{Tag}())
        setnext!(prev, head, Val{Tag}())
        list.tail = prev
    end
    tail
end

function Base.pushfirst!(list::TaggedIntrusiveSList{T, Tag}, node::T) where {T, Tag}
    if isempty(list)
        setnext!(node, node, Val{Tag}())
        list.tail = node
    else
        tail = list.tail::T
        head = getnext(tail, Val{Tag}())::T
        setnext!(tail, node, Val{Tag}())
        setnext!(node, head, Val{Tag}())
    end
    list
end

function Base.popfirst!(list::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    @boundscheck checkbounds(list)
    tail = list.tail::T
    head = getnext(tail, Val{Tag}())::T
    if head === tail
        list.tail = nothing
    else
        setnext!(tail, getnext(head, Val{Tag}())::T, Val{Tag}())
    end
    head
end

function deleteafter!(list::TaggedIntrusiveSList{T, Tag}, node::T) where {T, Tag}
    next = getnext(node, Val{Tag}())::T
    setnext!(node, getnext(next, Val{Tag}())::T, Val{Tag}())
    if next === list.tail
        if next === node
            list.tail = nothing
        else
            list.tail = node
        end
    end
    list
end

function insertafter!(list::TaggedIntrusiveSList{T, Tag}, after::T, node::T) where {T, Tag}
    next = getnext(after, Val{Tag}())::T
    setnext!(after, node, Val{Tag}())
    setnext!(node, next, Val{Tag}())
    if after === list.tail
        list.tail = node
    end
    list
end

function Base.delete!(list::TaggedIntrusiveSList{T, Tag}, node::T) where {T, Tag}
    @boundscheck checkbounds(list)
    prev = circslist_prev(node, list.tail::T, Val{Tag}())
    deleteafter!(list, prev)
    list
end

function Base.append!(list::TaggedIntrusiveSList{T, Tag}, list2::TaggedIntrusiveSList{T, Tag}) where {T, Tag}
    isempty(list2) && return list
    if !isempty(list)
        tail = list.tail::T
        tail2 = list2.tail::T
        head = getnext(tail, Val{Tag}())::T
        head2 = getnext(tail2, Val{Tag}())::T
        setnext!(tail, head2, Val{Tag}())
        setnext!(tail2, head, Val{Tag}())
    end
    list.tail = list2.tail::T
    list2.tail = nothing
    list
end


text:
# Calculate to total volume of regions

using Plots
include("../mustache-tables.jl")
include("load-one.jl")

# regions = [11]

# # assume independence: gives super small error
# samples = 200
# inds = 1:100
# totvol = zeros(length(inds))
# for out in [out1, out2]
#     for o in out
#         !(o isa Dict) && continue
#         vs = o[:mc_vols_above]
#         for i=inds
#             totvol[i] += vs[i]./1e9 #km3
#         end
#     end
# end
# Plots.histogram(totvol)

# # do not assume independence
# totvol = 0.0;
# totvol_std = 0.0;
# for out in [out1, out2]
#     for o in out
#         !(o isa Dict) && continue
#         # @show o[:vol_above_km3][:mean], std(o[:mc_vols_above])/1e9
#         totvol += o[:vol_above_km3][:mean]
#         totvol_std += std(o[:mc_vols_above])/1e9
#     end
# end
# println("Region $region total vol $(round(Int, totvol)) ± $(round(Int, totvol_std)) km3")

using Measurements, DataStructures

"Values from the G2TI paper"
G2TI = OrderedDict(3 => 28.33±7.35,
                   4 => 8.61±2.23,
                   7 => 7.47±1.94,
                   11 =>0.13±0.03,
#                   14 =>2.87±0.74)  # note, this is for both subregions
                   14 => 2.19±0.56) # this is for subregion 2 only (error is just adjusted from whole)

"""
Total volume in 10^3km^3

The errors (stdev) are just added.  I.e. no independence is assumed
For glaciers which errored it uses the Volume-Area volume instead (if
use_VA=true, default).

Returns:
- total volume with error
- the error-glaciers' V-A volume (should be small)
- total number of glaciers
- number of error glaciers
- total area
- area of error glaciers
"""
function get_vols(o1s, o2s, gnr1, gnr2; use_VA=true)
    out = Dict()
    va_out = Dict()
    nerrs = Dict()
    aerrs = Dict()
    ntot = Dict()
    atot = Dict()
    for reg in keys(o1s)
        glnrs, areas = BM.read_attrib_file(reg)[[1, 5]]
        area = Dict(g=>a/1e6 for (g,a) in zip(glnrs, areas))

        totvol = 0.0;
        totvol_std = 0.0;
        va_vol = 0.0;
        va_vol_std = 0.0;
        nerrs[reg] = 0
        aerrs[reg] = 0.0
        ntot[reg] = length(o1s[reg]) + length(o2s[reg])
        atot[reg] = sum(values(area))
        if reg==14
            # only subregion 2
            atot[reg] = 22871
        end

        for (oo,nr) in zip([o1s[reg], o2s[reg]], [gnr1[reg], gnr2[reg]])
            for (o,n) in zip(oo,nr)
                if !(o isa Dict)
                    a = area[n]
                    if use_VA
                        v = BM.volume_area_mean_h(a*1e6)*a*1e6 /1e9/1000
                        stdv = 0.2*v
                        totvol += v
                        totvol_std += stdv
                        va_vol += v
                        va_vol_std += stdv
                    end
                    nerrs[reg] += 1
                    aerrs[reg] += a
                else
                    # @show o[:vol_above_km3][:mean], std(o[:mc_vols_above])/1e9
                    totvol += o[:vol_above_km3][:mean]/1000
                    totvol_std += std(o[:mc_vols_above])/1e9/1000
                end
            end
        end
        out[reg] = totvol ± totvol_std
        va_out[reg] = va_vol ± va_vol_std
    end
    return out, va_out, ntot, nerrs, atot, aerrs
end

# https://github.com/JuliaPlots/Plots.jl/issues/143
# https://github.com/JuliaPlots/StatsPlots.jl/issues/32

"""
Compares the regional G2TI estimates with ours.
"""
function plot_compare_G2TI(vols; kws...)
    ks = [k for k in keys(G2TI) if haskey(vols, k)]
    vs = hcat([[vols[k].val, G2TI[k].val] for k in ks]...)
    errs = hcat([[vols[k].err, G2TI[k].err] for k in ks]...)
    ind11 = findfirst(ks.==11)
    # p1= StatPlots.groupedbar(map(VAWTools.int2str2, repeat(ks,inner=2)), vs', err=errs', marker=stroke(2, :black),
    #                      xlabel="RGI region", ylabel="V (10^3 km^3)", label="") #group=repeat(["BITE", "G2TI"], inner=length(ks)))

    inds1 = [findfirst(ks.==r) for r in [3,4,7,14] if findfirst(ks.==r)>0]
    inds2 = [findfirst(ks.==r) for r in [11] if findfirst(ks.==r)>0]
    p1= StatPlots.groupedbar(map(VAWTools.int2str2, repeat(ks[inds1],inner=2)), vs[:,inds1]', err=errs[:,inds1]',
                             markercolor=:black, markerstrokecolor=:black, markerstrokewidth=2,
                             xlabel="RGI region", ylabel="Volume (10³ km³)", label="", bar_width=0.7); #group=repeat(["BITE", "G2TI"], inner=length(ks)))

    p0 = StatPlots.groupedbar!([1,2], NaN*vs[:,1:2]', err=NaN*errs[:,1:2]', group=["BITE", "G2TI"]);

    p2 = StatPlots.groupedbar(map(VAWTools.int2str2, repeat(ks[inds2],inner=2)), vs[:,inds2]', err=errs[:,inds2]',
                              markercolor=:black, markerstrokecolor=:black, markerstrokewidth=2,
                              xlabel="", ylabel="", bar_width=0.7, label=""); #, group=repeat(["BITE", "G2TI"], inner=length(ks[inds2])))

    plot(p1, p2, p0, layout=@layout [a{0.8w} c{0.15w} b{0.05w}]; kws...)
#        plot(p1, p2, layout=@layout [a{0.8w} c{0.15w}])
end

function plot_compare_G2TI_2(vols, volsdhdt, atot; kws...)
    ks = [k for k in keys(G2TI) if haskey(vols, k)]
    kks = ["$r ($k)" for (r,k) in zip(["A. Canada N.", "A. Canada S.", "Svalbard", "C. Europe", "Karakoram"], keys(G2TI)) if haskey(vols, k)]
    vs = hcat([[vols[k].val, volsdhdt[k].val, G2TI[k].val] for k in ks]...)
    errs = hcat([[vols[k].err, volsdhdt[k].err, G2TI[k].err] for k in ks]...)
    ind11 = findfirst(ks.==11)

    hmean = 1e6*hcat([[vols[k].val/atot[k] ± vols[k].err/atot[k], volsdhdt[k].val/atot[k] ± volsdhdt[k].err/atot[k], G2TI[k].val/atot[k] ± G2TI[k].err/atot[k]] for k in ks]...)
#    hmean_errs = hcat([[vols[k].err/atot[k], volsdhdt[k].err/atot[k], G2TI[k].err] for k in ks]...)

    # p1= StatPlots.groupedbar(map(VAWTools.int2str2, repeat(ks,inner=2)), vs', err=errs', marker=stroke(2, :black),
    #                      xlabel="RGI region", ylabel="V (10^3 km^3)", label="") #group=repeat(["BITE", "G2TI"], inner=length(ks)))

    inds1 = [findfirst(ks.==r) for r in [3,4,7,14] if findfirst(ks.==r)>0]
    inds2 = [findfirst(ks.==r) for r in [11] if findfirst(ks.==r)>0]
    p1= StatPlots.groupedbar(repeat(kks[inds1],inner=3), vs[:,inds1]', err=errs[:,inds1]',
                             markercolor=:black, markerstrokecolor=:black, markerstrokewidth=1.5,
                             xlabel="RGI region", ylabel="Volume (10³ km³)", label="", bar_width=0.7); #group=repeat(["BITE", "G2TI"], inner=length(ks)))

    off = 0.22
    font = text("").font
    font.rotation = pi/2
    font.pointsize = 10
    for (i,ind) in enumerate(inds1)
        for j=1:3
            isnan(hmean[j,ind]) && continue
            txt = text("$(round(Int, hmean[j,ind]))±$(round(Int, hmean[j,ind].err))m", font)
            annotate!((i-0.5) + off*(j-2), 15, txt)
        end
    end

    #p0 = StatPlots.groupedbar!([1,2,3], NaN*vs[:,1:3]', err=NaN*errs[:,1:3]', group=["BITE", "BITE-dhdt", "G2TI"]);
    p0 = StatPlots.groupedbar([1,2,3], NaN*vs[:,1:3]', err=NaN*errs[:,1:3]', group=["BITE", "BITE-dhdt", "G2TI"], frame=:none);

    p2 = StatPlots.groupedbar(repeat(kks[inds2],inner=3), vs[:,inds2]', err=errs[:,inds2]',
                              markercolor=:black, markerstrokecolor=:black, markerstrokewidth=1.5,
                              xlabel="", ylabel="", bar_width=0.7, label=""); #, group=repeat(["BITE", "G2TI"], inner=length(ks[inds2])))

    for (i,ind) in enumerate(inds2)
        for j=1:3
            isnan(hmean[j,ind]) && continue
            txt = text("$(round(Int, hmean[j,ind]))±$(round(Int, hmean[j,ind].err))m", font)
            annotate!((i-0.5) + off*(j-2), 0.07, txt)
        end
    end

#    plot(p1, p2, layout=@layout [a{0.8w} c{0.15w}]; kws...)
    plot(p1, p0, p2, layout=@layout [a{0.8w} b{0.001w} c{0.15w}]; kws...)
#        plot(p1, p2, layout=@layout [a{0.8w} c{0.15w}])
end

""
function split_camel_case(cc)
    out = []
    for c in cc
        if isupper(c)
            push!(out, ' ')
        end
        push!(out, c)
    end
    out = "$(out...)"
    strip(out)
end

"Total area of regions"
areas = Dict(3=>104920,
             4=>40860,
             7=>33932,
             11=>2091,
             14=>33561)

function table_compare_G2TI(vols)
    ks = [k for k in keys(G2TI) if haskey(vols, k)]

    # names = [:Region, :Volume, :Error, :Volume_G2TI, :Error_G2TI]
    # vals = hcat([Any[k, vols[k].val, vols[k].err, G2TI[k].val, G2TI[k].err] for k in ks]...)'
    # df = DataFrame(vals, names)

    header = ["RGI Region", "N", "A (km²)", "V (10³ km³)", "\$\\bar\{\\mathrm\{h\}\}\$ (m)", "G2TI V (10³ km³)", "G2TI \$\\bar\{\\mathrm\{h\}\}\$ (m)"]
    names = Symbol.(["RGI Region", "N", "A (km²)", "V (10³ km³)", "h", "G2TI V (10³ km³)", "G2TI h (m)"])
    Vals = [Any[replace(replace(split_camel_case(string(BM.reg2dir[k])), "North", "N"), "South", "S") * " ($k)",
                length(get_glaciers(k)[3]),
                areas[k],
                map(x->signif(x.val, 3)±signif(x.err, 3), vols[k]),
#                signif.(vols[k], 3),
                round.(vols[k]*1e6./areas[k]),
                signif.(G2TI[k], 3),
                round.(G2TI[k]*1e6./areas[k])
                ] for k in ks]
    vals = permutedims(hcat(vals...), [2,1])
    return DataFrame(vals, names), header
end

"""
Plot the distribution of parameters
"""
function plot_paras(out1; kws...)
    exp_dist = hcat([o[:thetas_expect] for o in out1 if o isa Dict]...)
    mean_exp = calculate_prior_update(out1, false)[3]
    println(mean_exp)
    pls = []
    for (n,inds) in out1[1][:theta0][:names]
        for i in inds
            yl = i==1 ? BM.getrgi_region(out1[1][:rgi]) : ""
            push!(pls, Plots.histogram(exp_dist[i,:], xlabel=n, label="", ylabel=yl))
            Plots.scatter!([mean_exp[i]], [0], label="")
        end
    end
    Plots.plot(pls...; layout=grid(1,length(mean_exp)), kws...)
end

function plot_paras(o1s::Dict)
    pls = []
    for (k,o) in o1s
        if !isempty(o)
            push!(pls, plot_paras(o))
        end
    end
    Plots.plot(pls..., layout=grid(length(pls),1))
end

## Takes TIME
##       3,4      4,7,11      14
# dirbase = "output/rasters"
# dirs = ["mar-23", "apr-02", "apr-03"]
# o1s, o2s = Dict(), Dict()
# gnr1s, gnr2s = Dict(), Dict()
# for dir in dirs
#     @show dir
#     o1, o2, gnr1, gnr2 = load_all(dir, dirbase)
#     merge!(o1s, o1)
#     merge!(o2s, o2)
#     merge!(gnr1s, gnr1)
#     merge!(gnr2s, gnr2)
# end

# nosigma, withsigma = nothing, nothing
withsigma = load_all("production-without-dhdt", "output/rasters-fit_sigma");#, regions=[11])
withsigma_dhdt_ = load_all("production-with-dhdt", "output/rasters-fit_sigma");#, regions=[11])

# o1s, o2s, gnr1s, gnr2s = [nosigma, withsigma, withsigma_dhdt_][2];

volsdhdt, va_vols, ntot, nerrs, atot, aerrs = get_vols(withsigma_dhdt_...);
vols, va_vols, ntot, nerrs, atot, aerrs = get_vols(withsigma...);

volsdhdt[3] = volsdhdt[7]*NaN
volsdhdt[4] = volsdhdt[7]*NaN
if !haskey(volsdhdt, 14)
    volsdhdt[14] = volsdhdt[7]*NaN
end
# df_vols, header = table_compare_G2TI(vols)
# print(df_to_table(df_vols, header=header, fmt="lrrcrcr"))

display(plot_compare_G2TI_2(vols,volsdhdt,atot))

#display(plot_compare_G2TI(vols))
#savefig("paper1/figs/g2ti_comp.png")
savefig("figs/g2ti_comp.pdf")

error()

## get expectation of region 14
em = []
ss = []
for r in [3,4,7,11]
    if r in [3,4]
        o1s = withsigma[1]
    else
        o1s = withsigma_dhdt_[1]
    end
    @show r
    @show o1s[r][1][:thetas_expect]
    @show o1s[r][1][:theta0][:names]
    if !haskey(o1s, r)
        println("Not using region $r")
        continue
    end
    _,expect_median,_,stdev = calculate_prior_update(o1s[r], false)
    push!(em, expect_median)
    push!(ss, stdev)
end
println(median(hcat(em...),2)[:])
println(mean(hcat(em...),2)[:])
println(std(hcat(ss...),2)[:])


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

using DelimitedFiles, Random, Sockets

mktempdir() do dir

tasks = []

# Create test file
filename = joinpath(dir, "file.txt")
text = "C1,C2\n1,2\na,b\n"

# List of IO producers
l = Vector{Tuple{AbstractString,Function}}()

# File
io = (text) -> begin
    write(filename, text)
    Base.Filesystem.open(filename, Base.Filesystem.JL_O_RDONLY)
end
s = io(text)
@test isa(s, IO)
@test isa(s, Base.Filesystem.File)
close(s)
push!(l, ("File", io))

# IOStream
io = (text) -> begin
    write(filename, text)
    open(filename)
end
s = io(text)
@test isa(s, IO)
@test isa(s, IOStream)
close(s)
push!(l, ("IOStream", io))

# IOBuffer
io = (text)->IOBuffer(text)
s = io(text)
@test isa(s, IO)
@test isa(s, IOBuffer)
close(s)
push!(l, ("IOBuffer", io))

function run_test_server(srv, text)
    push!(tasks, @async begin
        try
            sock = accept(srv)
            try
                write(sock, text)
            catch e
                if !(isa(e, Base.UVError) && e.code == Base.UV_EPIPE)
                    if !(isa(e, Base.UVError) && e.code == Base.UV_ECONNRESET)
                        rethrow(e)
                    end
                end
            finally
                close(sock)
            end
        finally
            close(srv)
        end
    end)
    yield()
end

# TCPSocket
io = (text) -> begin
    port, srv = listenany(rand(2000:4000))
    run_test_server(srv, text)
    connect(port)
end
s = io(text)
@test isa(s, IO)
@test isa(s, TCPSocket)
close(s)
push!(l, ("TCPSocket", io))

# PipeEndpoint
io = (text) -> begin
    a = "\\\\.\\pipe\\uv-test-$(randstring(6))"
    b = joinpath(dir, "socket-$(randstring(6))")
    socketname = Sys.iswindows() ? a : b
    srv = listen(socketname)
    run_test_server(srv, text)
    connect(socketname)
end
s = io(text)
@test isa(s, IO)
@test isa(s, Base.PipeEndpoint)
close(s)
push!(l, ("PipeEndpoint", io))

#FIXME See https://github.com/JuliaLang/julia/issues/14747
#      Reading from open(::Command) seems to deadlock on Linux/Travis
#=
if !Sys.iswindows()

# Windows type command not working?
# See "could not spawn `type 'C:\Users\appveyor\AppData\Local\Temp\1\jul3516.tmp\file.txt'`"
#https://ci.appveyor.com/project/StefanKarpinski/julia/build/1.0.12733/job/hpwjs4hmf03vs5ag#L1244

# Pipe
io = (text) -> begin
    write(filename, text)
    open(`$(Sys.iswindows() ? "type" : "cat") $filename`)[1]
#    Was open(`echo -n $text`)[1]
#    See https://github.com/JuliaLang/julia/issues/14747
end
s = io(text)
@test isa(s, IO)
@test isa(s, Pipe)
close(s)
push!(l, ("Pipe", io))

end
=#

open_streams = []
function cleanup()
    for s_ in open_streams
        try close(s_) end
    end
    empty!(open_streams)
    for tsk in tasks
        Base._wait(tsk)
    end
    empty!(tasks)
end

verbose = false

for (name, f) in l
    local f
    local function io(text=text)
        local s = f(text)
        push!(open_streams, s)
        return s
    end

    verbose && println("$name readuntil...")
    for (t, s, m, kept) in [
            ("a", "ab", "a", "a"),
            ("b", "ab", "b", "b"),
            ("α", "αγ", "α", "α"),
            ("ab", "abc", "ab", "ab"),
            ("bc", "abc", "bc", "bc"),
            ("αβ", "αβγ", "αβ", "αβ"),
            ("aaabc", "ab", "aa", "aaab"),
            ("aaabc", "ac", "aaabc", "aaabc"),
            ("aaabc", "aab", "a", "aaab"),
            ("aaabc", "aac", "aaabc", "aaabc"),
            ("αααβγ", "αβ", "αα", "αααβ"),
            ("αααβγ", "ααβ", "α", "αααβ"),
            ("αααβγ", "αγ", "αααβγ", "αααβγ"),
            ("barbarbarians", "barbarian", "bar", "barbarbarian"),
            ("abcaabcaabcxl", "abcaabcx", "abca", "abcaabcaabcx"),
            ("abbaabbaabbabbaax", "abbaabbabbaax", "abba", "abbaabbaabbabbaax"),
            ("abbaabbabbaabbaabbabbaax", "abbaabbabbaax", "abbaabbabba", "abbaabbabbaabbaabbabbaax"),
           ]
        local t, s, m, kept
        @test readuntil(io(t), s) == m
        @test readuntil(io(t), s, keep=true) == kept
        @test readuntil(io(t), SubString(s, firstindex(s))) == m
        @test readuntil(io(t), SubString(s, firstindex(s)), keep=true) == kept
        @test readuntil(io(t), GenericString(s)) == m
        @test readuntil(io(t), GenericString(s), keep=true) == kept
        @test readuntil(io(t), unsafe_wrap(Vector{UInt8},s)) == unsafe_wrap(Vector{UInt8},m)
        @test readuntil(io(t), unsafe_wrap(Vector{UInt8},s), keep=true) == unsafe_wrap(Vector{UInt8},kept)
        @test readuntil(io(t), collect(s)::Vector{Char}) == Vector{Char}(m)
        @test readuntil(io(t), collect(s)::Vector{Char}, keep=true) == Vector{Char}(kept)
    end
    cleanup()

    write(filename, text)

    verbose && println("$name read...")
    @test read(io(), UInt8) == read(IOBuffer(text), UInt8)
    @test read(io(), UInt8) == read(filename, UInt8)
    @test read(io(), Int) == read(IOBuffer(text), Int)
    @test read(io(), Int) == read(filename,Int)
    s1 = io()
    s2 = IOBuffer(text)
    @test read!(s1, Vector{UInt32}(uninitialized, 2)) == read!(s2, Vector{UInt32}(uninitialized, 2))
    @test !eof(s1)
    @test read!(s1, Vector{UInt8}(uninitialized, 5)) == read!(s2, Vector{UInt8}(uninitialized, 5))
    @test !eof(s1)
    @test read!(s1, Vector{UInt8}(uninitialized, 1)) == read!(s2, Vector{UInt8}(uninitialized, 1))
    @test eof(s1)
    @test_throws EOFError read(s1, UInt8)
    @test eof(s1)
    close(s1)
    close(s2)

    verbose && println("$name eof...")
    n = length(text) - 1
    @test read!(io(), Vector{UInt8}(uninitialized, n)) ==
          read!(IOBuffer(text), Vector{UInt8}(uninitialized, n))
    @test (s = io(); read!(s, Vector{UInt8}(uninitialized, n)); !eof(s))
    n = length(text)
    @test read!(io(), Vector{UInt8}(uninitialized, n)) ==
          read!(IOBuffer(text), Vector{UInt8}(uninitialized, n))
    @test (s = io(); read!(s, Vector{UInt8}(uninitialized, n)); eof(s))
    n = length(text) + 1
    @test_throws EOFError read!(io(), Vector{UInt8}(uninitialized, n))
    @test_throws EOFError read!(io(), Vector{UInt8}(uninitialized, n))

    old_text = text
    cleanup()

    for text_ in [
        old_text,
        String(Char['A' + i % 52 for i in 1:(div(Base.SZ_UNBUFFERED_IO,2))]),
        String(Char['A' + i % 52 for i in 1:(    Base.SZ_UNBUFFERED_IO -1)]),
        String(Char['A' + i % 52 for i in 1:(    Base.SZ_UNBUFFERED_IO   )]),
        String(Char['A' + i % 52 for i in 1:(    Base.SZ_UNBUFFERED_IO +1)])
    ]
        text = text_
        write(filename, text)

        verbose && println("$name read(io, String)...")
        @test read(io(), String) == text

        @test read(io(), String) == read(filename, String)

        verbose && println("$name read...")
        @test read(io()) == unsafe_wrap(Vector{UInt8},text)

        @test read(io()) == read(filename)

        cleanup()

        verbose && println("$name readbytes!...")
        l = length(text)
        for n = [1, 2, l-2, l-1, l, l+1, l+2]
            a1 = Vector{UInt8}(uninitialized, n)
            a2 = Vector{UInt8}(uninitialized, n)
            s1 = io()
            s2 = IOBuffer(text)
            n1 = readbytes!(s1, a1)
            n2 = readbytes!(s2, a2)
            @test n1 == n2
            @test length(a1) == length(a2)
            @test a1[1:n1] == a2[1:n2]
            @test n <= length(text) || eof(s1)
            @test n <= length(text) || eof(s2)

            cleanup()
        end

        verbose && println("$name read!...")
        l = length(text)
        for n = [1, 2, l-2, l-1, l]
            @test read!(io(), Vector{UInt8}(uninitialized, n)) ==
                  read!(IOBuffer(text), Vector{UInt8}(uninitialized, n))
            @test read!(io(), Vector{UInt8}(uninitialized, n)) ==
                  read!(filename, Vector{UInt8}(uninitialized, n))

            cleanup()
        end
        @test_throws EOFError read!(io(), Vector{UInt8}(uninitialized, length(text)+1))

        verbose && println("$name readuntil...")
        for keep in [false, true]
            @test readuntil(io(), '\n', keep=keep) == readuntil(IOBuffer(text),'\n', keep=keep)
            @test readuntil(io(), '\n', keep=keep) == readuntil(filename,'\n', keep=keep)
            @test readuntil(io(), "\n", keep=keep) == readuntil(IOBuffer(text),"\n", keep=keep)
            @test readuntil(io(), "\n", keep=keep) == readuntil(filename,"\n", keep=keep)
            @test readuntil(io(), ',', keep=keep)  == readuntil(IOBuffer(text),',', keep=keep)
            @test readuntil(io(), ',', keep=keep)  == readuntil(filename,',', keep=keep)
        end

        cleanup()

        verbose && println("$name readline...")
        @test readline(io(), keep=true) == readline(IOBuffer(text), keep=true)
        @test readline(io(), keep=true) == readline(filename, keep=true)

        verbose && println("$name readlines...")
        @test readlines(io(), keep=true) == readlines(IOBuffer(text), keep=true)
        @test readlines(io(), keep=true) == readlines(filename, keep=true)
        @test readlines(io()) == readlines(IOBuffer(text))
        @test readlines(io()) == readlines(filename)
        @test collect(eachline(io(), keep=true)) == collect(eachline(IOBuffer(text), keep=true))
        @test collect(eachline(io(), keep=true)) == collect(eachline(filename, keep=true))
        @test collect(eachline(io())) == collect(eachline(IOBuffer(text)))
        @test collect(@inferred(eachline(io()))) == collect(@inferred(eachline(filename))) #20351

        cleanup()

        verbose && println("$name countlines...")
        @test countlines(io()) == countlines(IOBuffer(text))

        verbose && println("$name readdlm...")
        @test readdlm(io(), ',') == readdlm(IOBuffer(text), ',')
        @test readdlm(io(), ',') == readdlm(filename, ',')

        cleanup()
    end

    text = old_text
    write(filename, text)

    if !(typeof(io()) in [Base.PipeEndpoint, Pipe, TCPSocket])
        verbose && println("$name position...")
        @test (s = io(); read!(s, Vector{UInt8}(uninitialized, 4)); position(s))  == 4

        verbose && println("$name seek...")
        for n = 0:length(text)-1
            @test readlines(seek(io(), n)) == readlines(seek(IOBuffer(text), n))
            cleanup()
        end
        verbose && println("$name skip...")
        for n = 0:length(text)-1
            @test readlines(seek(io(), n)) == readlines(seek(IOBuffer(text), n))
            @test readlines(skip(io(), n)) == readlines(skip(IOBuffer(text), n))
            cleanup()
        end
        verbose && println("$name seekend...")
        @test read(seekend(io()), String) == ""
    end

    verbose && println("$name write(::IOStream, ...)")
    to = open("$filename.to", "w")
    write(to, io())
    close(to)
    @test read("$filename.to", String) == text

    verbose && println("$name write(filename, ...)")
    write("$filename.to", io())
    @test read("$filename.to", String) == text

    verbose && println("$name write(::IOBuffer, ...)")
    to = IOBuffer(Vector{UInt8}(codeunits(text)), read=false, write=true)
    write(to, io())
    @test String(take!(to)) == text

    cleanup()
end

function test_read_nbyte()
    fn = tempname()
    # Write one byte. One byte read should work once
    # but 2-byte read should throw EOFError.
    open(fn, "w+") do f
        write(f, 0x55)
        flush(f)
        seek(f, 0)
        @test read(f, UInt8) == 0x55
        @test_throws EOFError read(f, UInt8)
        seek(f, 0)
        @test_throws EOFError read(f, UInt16)
    end
    # Write 2 more bytes. Now 2-byte read should work once
    # but 4-byte read should fail with EOFError.
    open(fn, "a+") do f
        write(f, 0x4444)
        flush(f)
        seek(f, 0)
        @test read(f, UInt16) == 0x4455
        @test_throws EOFError read(f, UInt16)
        seek(f, 0)
        @test_throws EOFError read(f, UInt32)
    end
    # Write 4 more bytes. Now 4-byte read should work once
    # but 8-byte read should fail with EOFError.
    open(fn, "a+") do f
        write(f, 0x33333333)
        flush(f)
        seek(f, 0)
        @test read(f, UInt32) == 0x33444455
        @test_throws EOFError read(f, UInt32)
        seek(f, 0)
        @test_throws EOFError read(f, UInt64)
    end
    # Writing one more byte should allow an 8-byte
    # read to proceed.
    open(fn, "a+") do f
        write(f, 0x22)
        flush(f)
        seek(f, 0)
        @test read(f, UInt64) == 0x2233333333444455
    end
    rm(fn)
end
test_read_nbyte()

# devnull
@test !isreadable(devnull)
@test iswritable(devnull)
@test isopen(devnull)
@test write(devnull, 0xff) === 1
@test write(devnull, Int32(1234)) === 4
@test_throws EOFError read(devnull, UInt8)
@test close(devnull) === nothing
@test flush(devnull) === nothing
@test eof(devnull)
@test print(devnull, "go to /dev/null") === nothing

let s = "qwerty"
    @test read(IOBuffer(s)) == codeunits(s)
    @test read(IOBuffer(s), 10) == codeunits(s)
    @test read(IOBuffer(s), 1) == codeunits(s)[1:1]

    # Test growing output array
    x = UInt8[]
    n = readbytes!(IOBuffer(s), x, 10)
    @test x == codeunits(s)
    @test n == length(x)
end

# Filesystem.File
f = joinpath(dir, "test.txt")
open(io->write(io, "123"), f, "w")
f1 = open(f)
f2 = Base.Filesystem.open(f, Base.Filesystem.JL_O_RDONLY)
@test read(f1, UInt8) == read(f2, UInt8) == UInt8('1')
@test read(f1, UInt8) == read(f2, UInt8) == UInt8('2')
@test read(f1, UInt8) == read(f2, UInt8) == UInt8('3')
@test_throws EOFError read(f1, UInt8)
@test_throws EOFError read(f2, UInt8)
close(f1)
close(f2)

a = UInt8[0,0,0]
f1 = open(f)
f2 = Base.Filesystem.open(f, Base.Filesystem.JL_O_RDONLY)
@test read!(f1, a) == read!(f2, a) == UInt8['1','2','3']
@test_throws EOFError read!(f1, a)
@test_throws EOFError read!(f2, a)
close(f1)
close(f2)

a = UInt8[0,0,0,0]
f1 = open(f)
f2 = Base.Filesystem.open(f, Base.Filesystem.JL_O_RDONLY)
@test_throws EOFError read!(f1, a)
@test_throws EOFError read!(f2, a)
close(f1)
close(f2)
rm(f)

io = Base.Filesystem.open(f, Base.Filesystem.JL_O_WRONLY | Base.Filesystem.JL_O_CREAT | Base.Filesystem.JL_O_EXCL, 0o000)
@test write(io, "abc") == 3
close(io)
if !Sys.iswindows() && get(ENV, "USER", "") != "root" && get(ENV, "HOME", "") != "/root"
    # msvcrt _wchmod documentation states that all files are readable,
    # so we don't test that it correctly set the umask on windows
    @test_throws SystemError open(f)
    @test_throws Base.UVError Base.Filesystem.open(f, Base.Filesystem.JL_O_RDONLY)
else
    Sys.iswindows() || @warn "File permissions tests skipped due to running tests as root (not recommended)"
    close(open(f))
end
chmod(f, 0o400)
f1 = open(f)
f2 = Base.Filesystem.open(f, Base.Filesystem.JL_O_RDONLY)
for i = 1:2
    @test !eof(f1)
    @test !eof(f2)
    @test position(f1) == 0
    @test position(f2) == 0
    @test read(f1, String) == read(f2, String) == "abc"
    @test read(f1, String) == read(f2, String) == ""
    @test position(f1) == 3
    @test position(f2) == 3
    @test eof(f1)
    @test eof(f2)
    @test seekstart(f1) == f1
    @test seekstart(f2) == f2
end
@test seekend(f1) == f1
@test seekend(f2) == f2
@test eof(f1)
@test eof(f2)
@test skip(f1, -2) == f1
@test skip(f2, -2) == f2
@test position(f1) == 1
@test position(f2) == 1
@test_throws SystemError skip(f1, -2)
@test_throws SystemError skip(f2, -2)
@test position(f1) == 1
@test position(f2) == 1
@test skip(f1, 300) == f1
@test skip(f2, 300) == f2
@test position(f1) == 301
@test position(f2) == 301
@test eof(f1)
@test eof(f2)
@test_throws ArgumentError write(f1, '*')
@test_throws Base.UVError write(f2, '*')
close(f1)
close(f2)
@test eof(f1)
@test_throws Base.UVError eof(f2)
if get(ENV, "USER", "") != "root" && get(ENV, "HOME", "") != "/root"
    @test_throws SystemError open(f, "r+")
    @test_throws Base.UVError Base.Filesystem.open(f, Base.Filesystem.JL_O_RDWR)
else
    @warn "File permissions tests skipped due to running tests as root (not recommended)"
end
chmod(f, 0o600)
f1 = open(f, "r+")
f2 = Base.Filesystem.open(f, Base.Filesystem.JL_O_RDWR)
@test skip(f1, 10) == f1
@test skip(f2, 10) == f2
@test eof(f1)
@test eof(f2)
@test write(f1, '*') == 1
@test flush(f1) === nothing
@test !eof(f2)
@test skip(f2, 1) == f2
@test write(f2, '*') == 1
@test !eof(f1)
@test seekstart(f1) == f1
@test seekstart(f2) == f2
@test read(f1, String) == read(f2, String) == "abc\0\0\0\0\0\0\0**"
close(f1)
close(f2)
rm(f)

end # mktempdir() do dir

@testset "countlines" begin
    @test countlines(IOBuffer("")) == 0
    @test countlines(IOBuffer("\n")) == 1
    @test countlines(IOBuffer("\n"), eol = '\r') == 1
    @test countlines(IOBuffer("\r\r\n\r"), eol = '\r') == 3
    @test countlines(IOBuffer("\n\n\n\n\n\n\n\n\n\n")) == 10
    @test countlines(IOBuffer("\n \n \n \n \n \n \n \n \n \n")) == 10
    @test countlines(IOBuffer("\r\n \r\n \r\n \r\n \r\n")) == 5
    @test countlines(IOBuffer("foo\nbar")) == length(readlines(IOBuffer("foo\nbar"))) == 2
    file = tempname()
    write(file,"Spiffy header\nspectacular first row\neven better 2nd row\nalmost done\n")
    @test countlines(file) == 4
    @test countlines(file, eol = '\r') == 1
    @test countlines(file, eol = '\n') == 4
    rm(file)
end

let p = Pipe()
    Base.link_pipe!(p, reader_supports_async=true, writer_supports_async=true)
    t = @schedule read(p)
    @sync begin
        @async write(p, zeros(UInt16, 660_000))
        for i = 1:typemax(UInt16)
            @async write(p, UInt16(i))
        end
        @async close(p.in)
    end
    s = reinterpret(UInt16, fetch(t))
    @test length(s) == 660_000 + typemax(UInt16)
    @test s[(end - typemax(UInt16)):end] == UInt16.(0:typemax(UInt16))
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
module OrgMode

export Org, @org_str,                 # types/org.jl, here
    org, term, tableofcontents, html, # render/*
    parsetree                         # analysis/parsetree.jl

include("types/org.jl")
include("cache/cache.jl")

include("parse/interpret.jl")

macro org_str(content::String)
    parse(Org, content)
end

include("types/documentation.jl")
include("parse/documentation.jl")

include("render/utils.jl")
include("render/formatting.jl")
include("render/org.jl")
include("render/term.jl")
include("render/html.jl")

include("analysis/parsetree.jl")
include("analysis/diff.jl")

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    AbstractGraph,
    UndirectedGraph, DirectedGraph, HyperGraph,
    OrderedGraph, OrderedDiGraph, OrderedHyperGraph,
    AbstractNode, AbstractEdge, UndirectedEdge, DirectedEdge,
    hascachestore, hascache, getcache, setcache!, clearcache!, @cachefirst,
    clone,
    neighbors, outneighbors, inneighbors,
    findedgekey, findalledgekeys, getedge, hasedge, nodeattr, edgeattr,
    adjacencies, successors, predecessors,
    incidences, outincidences, in_incidences,
    nodeset, edgeset, edgesiter, neighborsiter, nodeattrs, edgeattrs,
    nodecount, edgecount, neighborcount, degree, indegree, outdegree,
    addnode!, addedge!, setnodeattr!, setedgeattr!,
    unlinknodes, unlinkedges,
    nodeattrtype, edgeattrtype

# Graph types

abstract type AbstractGraph end
abstract type UndirectedGraph <: AbstractGraph end
abstract type DirectedGraph <: AbstractGraph end
abstract type HyperGraph <: AbstractGraph end
abstract type OrderedGraph <: UndirectedGraph end
abstract type OrderedDiGraph <: DirectedGraph end
abstract type OrderedHyperGraph <: HyperGraph end

# Node and edge attributes

abstract type AbstractNode end
abstract type AbstractEdge end
abstract type UndirectedEdge <: AbstractEdge end # TODO: unnecessary?
abstract type DirectedEdge <: AbstractEdge end # TODO: unnecessary?

"""
    hascachestore(graph::AbstractGraph) -> Bool

Return whether the graph has calculated property cache store.
"""
hascachestore(graph::AbstractGraph) = isdefined(graph, :cache)

"""
    hascache(graph::AbstractGraph, key) -> Bool

Return whether the graph has calculated properties with the given key or not.
"""
hascache(graph::AbstractGraph, key) = haskey(graph.cache, key)

"""
    getcache(graph::AbstractGraph, key) -> Any

Get cached calculated properties.
"""
getcache(graph::AbstractGraph, key) = graph.cache[key]

"""
    setcache!(graph::AbstractGraph, key; kwargs...)

Set calculated property caches.
"""
function setcache!(graph::AbstractGraph, key; kwargs...)
    graph.cache[key] = getfield(Graph, key)(graph; kwargs...)
end

"""
    clearcache!(graph::AbstractGraph)

Clear calculated property caches.

Calling `clearcache!` is recommended when the graph nodes/edges are added, removed or reindexed. You can `deepcopy` or [`Graph.clone`](@ref) the graph before destructive operation instead.
"""
clearcache!(graph::AbstractGraph) = empty!(graph.cache)

"""
    @cachefirst expression

A macro for molecular descriptor array cache mechanism.
"""
macro cachefirst(ex)
    func = ex.args[1].args[1]
    dummy = gensym()
    ex.args[1].args[1] = dummy
    return esc(quote
        $ex
        Core.@__doc__ function $func(graph; kwargs...)
            if hascachestore(graph) && isempty(kwargs)
                # Return cache
                symf = nameof($func)
                if hascache(graph, symf)
                    return getcache(graph, symf)
                end
            end
            return $dummy(graph; kwargs...)
        end
    end)
end

"""
    clone(graph::AbstractGraph) -> AbstractGraph

Return deep copy of the graph.
"""
clone = deepcopy  # TODO: use deepcopy

# Lookup

"""
    neighbors(graph, i) -> Dict

Return the mapping of incident edges and adjacent nodes of node `i`.
If the graph is directed graph, both outneighbors and inneighbors are mapped.
"""
neighbors(graph::UndirectedGraph, i) = graph.neighbormap[i]
neighbors(graph::DirectedGraph, i) = merge(outneighbors(graph, i), inneighbors(graph, i))

neighbors(graph::HyperGraph, i) = Dict(inc => graph.edges[inc] for inc in graph.incidences[i])

"""
    outneighbors(graph::DirectedGraph, i) -> Dict

Return the mapping of successor node keys and out edge keys connected to
the given node.
"""
outneighbors(graph::DirectedGraph, i) = graph.outneighbormap[i]

"""
    inneighbors(graph::DirectedGraph, i) -> Dict

Return the mapping of predecessor node keys and in edge keys connected to
the given node.
"""
inneighbors(graph::DirectedGraph, i) = graph.inneighbormap[i]

"""
    findedgekey(graph::UndirectedGraph, u, v) -> Any
    findedgekey(graph::DirectedGraph, u, v) -> Any
    findedgekey(graph::HyperGraph, u, v) -> Any

Return the first edge found that connects nodes u and v.
"""
function findedgekey(graph::UndirectedGraph, u, v)
    for (inc, adj) in neighbors(graph, u)
        adj == v && return inc
    end
end

function findedgekey(graph::DirectedGraph, source, target)
    for (inc, adj) in outneighbors(graph, source)
        adj == target && return inc
    end
end

function findedgekey(graph::HyperGraph, u, v)
    for (inc, adj) in neighbors(graph, u)
        v in adj && return inc
    end
end

"""
    findalledgekeys(graph::UndirectedGraph, u, v) -> Any
    findalledgekeys(graph::DirectedGraph, u, v) -> Any
    findalledgekeys(graph::HyperGraph, u, v) -> Any

Return all edges that connect nodes u and v.
"""
function findalledgekeys(graph::UndirectedGraph, u, v)
    hasmultiedge(graph) || return [findedgekey(graph, u, v)]
    keys = Int[]
    for (inc, adj) in neighbors(graph, u)
        adj == v && push!(keys, inc)
    end
    return keys
end

function findalledgekeys(graph::DirectedGraph, source, target)
    hasmultiedge(graph) || return [findedgekey(graph, source, target)]
    keys = Int[]
    for (inc, adj) in outneighbors(graph, source)
        adj == target && push!(keys, inc)
    end
    return keys
end

function findalledgekeys(graph::HyperGraph, u, v)
    hasmultiedge(graph) || return [findedgekey(graph, u, v)]
    keys = Int[]
    for (inc, adj) in neighbors(graph, u)
        v in adj && push!(keys, inc)
    end
    return keys
end

"""
    getedge(graph::AbstractGraph, i) -> Any

Return an edge.
"""
getedge(graph::AbstractGraph, i) = graph.edges[i]

"""
    hasedge(graph, u, v) -> Bool

Return whether the given two nodes are connected by at least one edge.
"""
hasedge(graph, u, v) = findedgekey(graph, u, v) !== nothing

"""
    nodeattr(graph::AbstractGraph, i) -> AbstractNode

Return the attribute object of node `i`.
"""
nodeattr(graph::AbstractGraph, i) = graph.nodeattrs[i]

"""
    edgeattr(graph::AbstractGraph, i) -> AbstractEdge

Return the attribute object of edge `i`.
"""
edgeattr(graph::AbstractGraph, i) = graph.edgeattrs[i]

"""
    edgeattr(graph::AbstractGraph, u, v) -> Union{AbstractEdge,Nothing}

Return the attribute object of an edge that connects `u` and `v`. If not found,
return nothing.
"""
function edgeattr(graph::AbstractGraph, u, v)
    k = findedgekey(graph, u, v)
    return k === nothing ? nothing : edgeattr(graph, k)
end

# Collection

adjacencies(g::AbstractGraph, i) = Set(values(neighbors(g, i)))
successors(g::DirectedGraph, i) = Set(values(outneighbors(g, i)))
predecessors(g::DirectedGraph, i) = Set(values(inneighbors(g, i)))

incidences(g::AbstractGraph, i) = Set(keys(neighbors(g, i)))
outincidences(g::DirectedGraph, i) = Set(keys(outneighbors(g, i)))
in_incidences(g::DirectedGraph, i) = Set(keys(inneighbors(g, i)))

"""
    nodeset(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) -> Set{Int}

Return the set of node keys.
"""
nodeset(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) = Set(1:nodecount(graph))

"""
    edgeset(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) -> Set{Int}

Return the set of edge keys.
"""
edgeset(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) = Set(1:edgecount(graph))

edgesiter(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) = graph.edges

neighborsiter(graph::OrderedGraph) = graph.neighbormap

"""
    nodeattrs(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) -> Vector{AbstractNode}

Return graph node attributes.
"""
nodeattrs(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) = graph.nodeattrs

"""
    edgeattrs(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) -> Vector{AbstractEdge}

Return graph edge attributes.
"""
edgeattrs(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}) = graph.edgeattrs

# Counting

"""
    nodecount(graph::AbstractGraph) -> Int

Return the number of graph nodes.
"""
nodecount(graph::UndirectedGraph) = length(graph.neighbormap)
nodecount(graph::DirectedGraph) = length(graph.outneighbormap)
nodecount(graph::HyperGraph) = length(graph.incidences)

"""
    edgecount(graph::AbstractGraph) -> Int

Return the number of graph edges.
"""
edgecount(graph::AbstractGraph) = length(graph.edges)

"""
    neighborcount(graph::AbstractGraph, n) -> Int
    degree(graph::AbstractGraph, n) -> Int

Return the number of adjacent nodes of the node 'n'.
"""
neighborcount(graph::AbstractGraph, n) = length(neighbors(graph, n))
degree = neighborcount

"""
    outdegree(graph::DirectedGraph, n) -> Int

Return the number of outneighbors of the node 'n'.
"""
outdegree(graph::DirectedGraph, n) = length(outneighbors(graph, n))

"""
    indegree(graph::DirectedGraph, n) -> Int

Return the number of inneighbors of the node 'n'.
"""
indegree(graph::DirectedGraph, n) = length(inneighbors(graph, n))

# Editing

"""
    addnode!(graph) -> Int
    addnode!(graph, attr) -> Int

Add new node and return the node index. If the node attribute type is required,
specify the node attribute object by `node` keyword.
"""
function addnode!(graph::UndirectedGraph, attr)
    push!(graph.neighbormap, Dict())
    push!(graph.nodeattrs, attr)
    return length(graph.neighbormap)
end

function addnode!(graph::UndirectedGraph)
    isdefined(graph, :nodeattrs) && throw(ErrorException("nodeattr required"))
    push!(graph.neighbormap, Dict())
    return length(graph.neighbormap)
end

function addnode!(graph::DirectedGraph, attr)
    push!(graph.outneighbormap, Dict())
    push!(graph.inneighbormap, Dict())
    push!(graph.nodeattrs, attr)
    return length(graph.outneighbormap)
end

function addnode!(graph::DirectedGraph)
    isdefined(graph, :nodeattrs) && throw(ErrorException("nodeattr required"))
    push!(graph.outneighbormap, Dict())
    push!(graph.inneighbormap, Dict())
    return length(graph.outneighbormap)
end

function addnode!(graph::HyperGraph, attr)
    push!(graph.incidences, Set())
    push!(graph.nodeattrs, attr)
    return length(graph.incidences)
end

function addnode!(graph::HyperGraph)
    isdefined(graph, :nodeattrs) && throw(ErrorException("nodeattr required"))
    push!(graph.incidences, Set())
    return length(graph.incidences)
end

"""
    addedge!(graph, u, v) -> Int
    addedge!(graph, u, v, attr) -> Int

Add new edge and return the edge index. If the edge attribute type is required,
specify the edge attribute object by `attr` keyword.
"""
function _addedge!(graph::UndirectedGraph, u, v)
    push!(graph.edges, (u, v))
    i = edgecount(graph)
    graph.neighbormap[u][i] = v
    graph.neighbormap[v][i] = u
    return i
end

function addedge!(graph::UndirectedGraph, u, v, attr)
    push!(graph.edgeattrs, attr)
    return _addedge!(graph, u, v)
end

function addedge!(graph::UndirectedGraph, u, v)
    isdefined(graph, :edgeattrs) && throw(ErrorException("edgeattr required"))
    return _addedge!(graph, u, v)
end

function _addedge!(graph::DirectedGraph, s, t)
    push!(graph.edges, (s, t))
    i = edgecount(graph)
    graph.outneighbormap[s][i] = t
    graph.inneighbormap[t][i] = s
    return i
end

function addedge!(graph::DirectedGraph, s, t, attr)
    push!(graph.edgeattrs, attr)
    return _addedge!(graph, s, t)
end

function addedge!(graph::DirectedGraph, s, t)
    isdefined(graph, :edgeattrs) && throw(ErrorException("edgeattr required"))
    return _addedge!(graph, s, t)
end

function _addedge!(graph::HyperGraph, edge)
    push!(graph.edges, edge)
    i = edgecount(graph)
    for n in edge
        push!(graph.incidences[n], i)
    end
    return i
end

function addedge!(graph::HyperGraph, edge, attr)
    push!(graph.edgeattrs, attr)
    return _addedge!(graph, edge)
end

function addedge!(graph::HyperGraph, edge)
    isdefined(graph, :edgeattrs) && throw(ErrorException("edgeattr required"))
    return _addedge!(graph, edge)
end

"""
    setnodeattr!(graph::AbstractGraph, i, attr::AbstractNode)

Update the node attribute.
"""
function setnodeattr!(graph::AbstractGraph, i, attr)
    graph.nodeattrs[i] = attr
    return
end

"""
    setedgeattr!(graph::AbstractGraph, i, attr::AbstractNode)

Update the edge attribute.
"""
function setedgeattr!(graph::AbstractGraph, i, attr)
    graph.edgeattrs[i] = attr
    return
end

"""
    unlinknodes(graph, nodes)

Delete given nodes and its incident edges from the graph.
"""
function unlinknodes end

"""
    unlinkedges(graph, edges)

Delete given edges from the graph.
"""
function unlinkedges end

# Types

"""
    nodeattrtype(graph) -> Type

Return node attribute type.
"""
nodeattrtype(mol::AbstractGraph) = eltype(mol.nodeattrs)

"""
    edgeattrtype(graph) -> Type

Return edge attribute type.
"""
edgeattrtype(mol::AbstractGraph) = eltype(mol.edgeattrs)


text:
using Game2048
using Test

@testset "Game2048.jl" begin
    @test true == true
end

board = initboard()
rand_dir = rand((left, right, up, down))
@time move!(board, rand_dir); board


text:
using CoordinateTransformations, TestImages, ImageCore, StaticArrays, OffsetArrays, Interpolations, LinearAlgebra
using Test, ReferenceTests

refambs = detect_ambiguities(CoordinateTransformations, Base, Core)
using ImageTransformations
ambs = detect_ambiguities(ImageTransformations, CoordinateTransformations, Base, Core)
@test isempty(setdiff(ambs, refambs))

tests = [
    "autorange.jl",
    "resizing.jl",
    "interpolations.jl",
    "warp.jl",
]

@testset "ImageTransformations" begin
for t in tests
    @testset "$t" begin
        include(t)
    end
end
end

nothing