text:
## Elliptical copulas

"""
    gausscopulagen(t::Int, Σ::Matrix{Float64} = [1. 0.5; 0.5 1.])

Returns: t x n matrix of t realisations of multivariate data generated
using gaussian copula with Σ - correlation matrix. If the symmetric covariance
matrix is imputed, it will be converted into a correlation matrix automatically.

```jldoctest

julia> Random.seed!(43);

julia> gausscopulagen(10)
10×2 Array{Float64,2}:
 0.589188  0.815308
 0.708285  0.924962
 0.747341  0.156994
 0.227634  0.183116
 0.227575  0.957376
 0.271558  0.364803
 0.445691  0.52792
 0.585362  0.23135
 0.498593  0.48266
 0.190283  0.594451
```
"""
function gausscopulagen(t::Int, Σ::Matrix{Float64} = [1. 0.5; 0.5 1.])
  z = transpose(rand(MvNormal(Σ),t))
  for i in 1:size(Σ, 1)
    d = Normal(0, sqrt.(Σ[i,i]))
    z[:,i] = cdf.(d, z[:,i])
  end
  Array(z)
end

"""
  tstudentcopulagen(t::Int, Σ::Matrix{Float64} = [1. 0.5; 0.5 1.], ν::Int=10)

Generates data using t-student Copula given Σ - correlation matrix, ν - degrees of freedom.
If the symmetric covariance matrix is imputed, it will be converted into a
correlation matrix automatically.

```jldoctest
julia> Random.seed!(43);

julia> tstudentcopulagen(10)
10×2 Array{Float64,2}:
 0.658199  0.937148
 0.718244  0.92602
 0.809521  0.0980325
 0.263068  0.222589
 0.187187  0.971109
 0.245373  0.346428
 0.452336  0.524498
 0.57113   0.272525
 0.498443  0.48082
 0.113788  0.633349
 ```
"""
function tstudentcopulagen(t::Int, Σ::Matrix{Float64} = [1. 0.5; 0.5 1.], ν::Int=10)
  z = transpose(rand(MvNormal(Σ),t))
  U = rand(Chisq(ν), size(z, 1))
  for i in 1:size(Σ, 1)
    x = z[:,i].*sqrt.(ν./U)./sqrt(Σ[i,i])
    z[:,i] = cdf.(TDist(ν), x)
  end
  Array(z)
end

### Frechet familly

"""
  frechetcopulagen(t::Int, n::Int, α::Union{Int, Float64})

Returns t realisation of n variate data generated from one parameter frechet multidimentional copula,
a combination of maximal copla with  weight α and independent copula with  weight 1-α

```jldoctest
julia> Random.seed!(43);

julia> frechetcopulagen(10, 2, 0.5)
10×2 Array{Float64,2}:
 0.180975  0.661781
 0.775377  0.775377
 0.888934  0.125437
 0.924876  0.924876
 0.408278  0.408278
 0.912603  0.740184
 0.828727  0.00463791
 0.400537  0.0288987
 0.429437  0.429437
 0.955881  0.851275
```
"""
function frechetcopulagen(t::Int, n::Int, α::Union{Int, Float64})
  0 <= α <= 1 || throw(DomainError("generaton not supported for α ∉ [0,1]"))
  u = rand(t, n)
  for j in 1:t
    if (α >= rand())
      for i in 1:n
        u[j,i] = u[j, 1]
      end
    end
  end
  u
end

"""
  frechetcopulagen(t::Int, n::Int, α::Union{Int, Float64}, β::Union{Int, Float64})

Two parameters Frechet copula C = α C_{max} + β C_{min} + (1- α - β) C_{⟂}, supported
only for n == 2

```jldoctest
julia> Random.seed!(43);

julia> frechetcopulagen(10, 2, 0.4, 0.2)
10×2 Array{Float64,2}:
 0.180975  0.661781
 0.775377  0.775377
 0.888934  0.125437
 0.924876  0.924876
 0.408278  0.591722
 0.912603  0.740184
 0.828727  0.171273
 0.400537  0.0288987
 0.429437  0.429437
 0.955881  0.851275
```
"""
function frechetcopulagen(t::Int, n::Int, α::Union{Int, Float64}, β::Union{Int, Float64})
  n == 2 || throw(AssertionError("two parameters Frechet copula supported only for n = 2"))
  0 <= α+β <= 1 || throw(DomainError("α+β must be in range [0,1]"))
  u = rand(t,2)
  for j in 1:t
    v = rand()
    if (α >= v)
      u[j,2] = u[j, 1]
    elseif (α < v <= α+β)
      u[j,2] = 1-u[j, 1]
    end
  end
  u
end

### Marshall - Olkin familly

"""
  marshallolkincopulagen(t::Int, λ::Vector{Float64})

Returns: t x n Matrix{Float}, t realisations of n-variate data generated from Marshall-Olkin
copula with parameter vector λ of non-negative elements λₛ.
Number of marginals is n = ceil(Int, log(2, length(λ)-1)).
Parameters are ordered as follow:
λ = [λ₁, λ₂, ..., λₙ, λ₁₂, λ₁₃, ..., λ₁ₙ, λ₂₃, ..., λₙ₋₁ₙ, λ₁₂₃, ..., λ₁₂...ₙ]
If reversed = true, returns data from reversed Marshall-Olkin copula.

```jldoctest

julia> Random.seed!(43)

julia> marshallolkincopulagen(10, [0.2, 1.2, 1.6])
10×2 Array{Float64,2}:
 0.99636   0.994344
 0.167268  0.0619408
 0.977418  0.965093
 0.495167  0.0247053
 0.410336  0.250159
 0.778989  0.678064
 0.50927   0.350059
 0.925875  0.887095
 0.353646  0.219006
 0.782477  0.686799
 ```
"""
function marshallolkincopulagen(t::Int, λ::Vector{Float64} = rand(7); reverse::Bool = false)
  minimum(λ) >= 0 || throw(AssertionError("all parameters must by >= 0 "))
  n = floor(Int, log(2, length(λ)+1))
  U = mocopula(rand(t,2^n-1), n, λ)
  reverse ? 1 .-U : U
end

"""
  mocopula(u::Matrix{Float64}, n::Int, λ::Vector{Float64})

  Returns: t x n Matrix{Float}, t realisations of n-variate data generated from Marshall-Olkin
  copula with parameter vector λ of non-negative elements λₛ, given [0,1]ᵗˣˡ ∋ u, where
  l = 2ⁿ-1

```jldoctest

  julia> mocopula([0.2 0.3 0.4; 0.3 0.4 0.6; 0.4 0.5 0.7], 2, [1., 1.5, 2.])
  3×2 Array{Float64,2}:
   0.252982  0.201189
   0.464758  0.409039
   0.585662  0.5357

```
"""
function mocopula(u::Matrix{Float64}, n::Int, λ::Vector{Float64})
  s = collect(combinations(1:n))
  t,l = size(u)
  U = zeros(t, n)
    for j in 1:t
      for i in 1:n
        inds = findall([i in s[k] for k in 1:l])
        x = minimum([-log(u[j,k])./(λ[k]) for k in inds])
        Λ = sum(λ[inds])
        U[j,i] = exp.(-Λ*x)
      end
    end
    U
end

# transforms univariate distributions
"""
  convertmarg!(X::Matrix, d::UnionAll, p::Vector{Vector})

Takes matrix X of realisations of size(X,2) = n dimensional random variable, with
uniform marginals numbered by i, and convert those marginals to common distribution
d with parameters p[i].
If `testunif = true` each marginal is tested for uniformity.

```jldoctest
julia> Random.seed!(43);

julia> x = rand(10,2);

julia> convertmarg!(x, Normal, [[0, 1],[0, 10]])

julia> x
10×2 Array{Float64,2}:
 -0.911655    4.17328
  0.756673  -14.4472
  1.22088   -11.4823
  1.43866   -13.1053
 -0.231978  -11.2415
  1.35696     6.43914
  0.949145  -26.0172
 -0.251957  -18.9723
 -0.177808    0.54172
  1.70477    10.4192
```
"""
VecVec = Union{Vector{Vector{Int64}}, Vector{Vector{Float64}}}

function convertmarg!(U::Matrix{T}, d::UnionAll, p::VecVec = [fill([0,1], size(U, 2))...];
                                                testunif::Bool = true) where T <: AbstractFloat
  for i = 1:size(U, 2)
    if testunif
      pvalue(ExactOneSampleKSTest(U[:,i],Uniform(0,1)))>0.0001 || throw(AssertionError("$i marg. not unif."))
    end
    @inbounds U[:,i] = quantile.(d(p[i]...), U[:,i])
  end
end

  # generates covariance matrix

"""
    cormatgen(n::Int, ρ::Float64 = 0.5, ordered::Bool = false, altersing::Bool = true)

Returns symmetric correlation matrix Σ of size n x n, with reference correlation 0 < ρ < 1.
If ordered = false, Σ elements varies arround ρ, i.e. σᵢⱼ ≈ ρ+δ else they drop
as indices differences rise, i.e. σᵢⱼ ≳ σᵢₖ as |i-j| < |i-k|.
If altersing = true, some σ are positive and some negative, else ∀ᵢⱼ σᵢⱼ ≥ 0.

```jldoctest
julia> Random.seed!(43);

julia> julia> cormatgen(2)
2×2 Array{Float64,2}:
 1.0       0.660768
 0.660768  1.0
```
"""
function cormatgen(n::Int = 20)
  a = rand(n,n)
  b = a*a'
  c = b./maximum(b)
  c .- Matrix(Diagonal(c)) .+ Matrix(1.0I, size(c)...)

end

function cormatgen_rand(n::Int = 20)
  a = rand(n,n)
  b = a*a'
  diagb = Matrix(Diagonal(1 ./sqrt.(LinearAlgebra.diag(b))))
  b = diagb*b*diagb
  (b+b')/2
end

function cormatgen_constant(n::Int, α::Float64)
  @assert 0 <= α <= 1 "α should satisfy 0 <= α <= 1"
  α .*ones(n, n) .+(1-α) .*Matrix(1.0I, n,n)
end

function random_unit_vector(dim::Int)
  result = rand(Normal(), dim, 1)
  result /= norm(result)
end

function cormatgen_constant_noised(n::Int, α::Float64; ϵ::Float64 = (1 .-α)/2.)
  @assert 0 <= ϵ <= 1-α "ϵ must satisfy 0 <= ϵ <= 1-α"
  result = cormatgen_constant(n, α)
  u = hcat([random_unit_vector(n) for i=1:n]...)
  result += ϵ .*(u'*u)
  result - ϵ .*Matrix(1.0I, size(result)...)
end

function cormatgen_two_constant(n::Int, α::Float64, β::Float64)
  @assert α > β "First argument must be greater"
  result = fill(β, (n,n))
  result[1:div(n,2),1:div(n,2)] = fill(α, (div(n,2),div(n,2)))
  result += Matrix(1.0I, size(result)...) - Matrix(Diagonal(result))
  result
end

function cormatgen_two_constant_noised(n::Int, α::Float64, β::Float64; ϵ::Float64= (1-α)/2)
  @assert ϵ <= 1-α
  result = cormatgen_two_constant(n, α, β)
  u = hcat([random_unit_vector(n) for i=1:n]...)
  result += ϵ .*(u'*u)
  result - ϵ .*Matrix(1.0I, size(result)...)
end

function cormatgen_toeplitz(n::Int, ρ::Float64)
  @assert 0 <= ρ <= 1 "ρ needs to satisfy 0 <= ρ <= 1"
  [ρ^(abs(i-j)) for i=0:n-1, j=0:n-1]
end

function cormatgen_toeplitz_noised(n::Int, ρ::Float64; ϵ=(1-ρ)/(1+ρ)/2)
  @assert 0 <= ϵ <= (1-ρ)/(1+ρ) "ϵ must satisfy 0 <= ϵ <= (1-ρ)/(1+ρ)"

  result = cormatgen_toeplitz(n, ρ)
  u = hcat([random_unit_vector(n) for i=1:n]...)
  result += ϵ .*(u'*u)
  result - ϵ .*Matrix(1.0I, size(result)...)
end


text:
function MLJBase.fit(model::MClassifier, verbosity::Int, X, Y)
    return fit(model, verbosity, X, Y)
end

function MLJBase.predict(classifier::MClassifier, fitresult, instance)
    return update_predict(classifier, fitresult, instance)
end

ext: jl


text:
include("markovDecision.jl")
include("simulate.jl")

board1 = [0,0,0,0,0,2,0,0,1,1,2,2,0,0,0]
board2 = [0,1,2,1,2,2,1,0,0,0,0,0,0,0,0]
board3 = ones(Int64, 15)

boards = [board1, board2, board3]
c = [false, true]

for board in boards
    println("--- Board : ", board, " ---")
    for circular in c
        if !circular
            println("Non-circular")
        else
            println("Circular")
        end
        (expec, dice) = markovdecision(board, circular)
        println("\tExpected cost = ", round(expec[1],3))
        allsim = Vector{Int64}()
        for i=1:1e5
        	push!(allsim, simulate(board, dice, circular))
        end
        println("\tSimulated cost = ", round(mean(allsim),3))
    end
end

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
using StaticArrays

const tetrahedron1an = [
    SVector{4}([0.2500000000000000  0.2500000000000000  0.2500000000000000  0.2500000000000000])
]

const tetrahedron1wn = [
  1.0,
]

#x = [0.13819660112501053]
#w = [ 0.25]

const tetrahedron4an = [
    SVector{4}([0.138196601125011 0.138196601125011 0.138196601125011 0.585410196624968]),
    SVector{4}([0.138196601125011 0.138196601125011 0.585410196624968 0.138196601125011]),
    SVector{4}([0.138196601125011 0.585410196624968 0.138196601125011 0.138196601125011]),
    SVector{4}([0.585410196624968 0.138196601125011 0.138196601125011 0.138196601125011])
]

const tetrahedron4wn = [
  0.25,
  0.25,
  0.25,
  0.25
]

#x = [0.07394411115584681, 0.40623507869707387]
#w = [0.047721831932829764, 0.13485211204478015]

const tetrahedron10an = [
  SVector{4}([0.0739441111558468 0.0739441111558468 0.0739441111558468 0.778167666532460]),
  SVector{4}([0.0739441111558468 0.0739441111558468 0.778167666532460 0.0739441111558468]),
  SVector{4}([0.0739441111558468 0.778167666532460 0.0739441111558468 0.0739441111558468]),
  SVector{4}([0.778167666532460 0.0739441111558468 0.0739441111558468 0.0739441111558468]),
  SVector{4}([0.406235078697074 0.406235078697074 0.0937649213029261 0.0937649213029261] ),
  SVector{4}([0.406235078697074 0.0937649213029261 0.406235078697074 0.0937649213029261] ),
  SVector{4}([0.406235078697074 0.0937649213029261 0.0937649213029261 0.406235078697074] ),
  SVector{4}([0.0937649213029261 0.406235078697074 0.406235078697074 0.0937649213029261] ),
  SVector{4}([0.0937649213029261 0.406235078697074 0.0937649213029261 0.406235078697074] ),
  SVector{4}([0.0937649213029261 0.0937649213029261 0.406235078697074 0.406235078697074] ),
]

const tetrahedron10wn = [
  0.047721831932829764,
  0.047721831932829764,
  0.047721831932829764,
  0.047721831932829764,
  0.13485211204478015,
  0.13485211204478015,
  0.13485211204478015,
  0.13485211204478015,
  0.13485211204478015,
  0.13485211204478015
]

#x = [0.09273339107792229, 0.31088554482739006, 0.0455085138515946, 0.4565216399641138]
#w = [0.07348669573167155, 0.1126851232260692, 0.02127606034741975]
const tetrahedron20an = [
  SVector{4}([0.0927333910779223 0.0927333910779223 0.0927333910779223 0.721799826766233]),
  SVector{4}([0.0927333910779223 0.0927333910779223 0.721799826766233 0.0927333910779223]),
  SVector{4}([0.0927333910779223 0.721799826766233 0.0927333910779223 0.0927333910779223]),
  SVector{4}([0.721799826766233 0.0927333910779223 0.0927333910779223 0.0927333910779223]),
  SVector{4}([0.310885544827390 0.310885544827390 0.310885544827390 0.0673433655178298]  ),
  SVector{4}([0.310885544827390 0.310885544827390 0.0673433655178298 0.310885544827390]  ),
  SVector{4}([0.310885544827390 0.0673433655178298 0.310885544827390 0.310885544827390]  ),
  SVector{4}([0.0673433655178298 0.310885544827390 0.310885544827390 0.310885544827390]  ),
  SVector{4}([0.0455085138515946 0.0455085138515946 0.456521639964114 0.452461332332697] ),
  SVector{4}([0.0455085138515946 0.0455085138515946 0.452461332332697 0.456521639964114] ),
  SVector{4}([0.0455085138515946 0.456521639964114 0.0455085138515946 0.452461332332697] ),
  SVector{4}([0.0455085138515946 0.456521639964114 0.452461332332697 0.0455085138515946] ),
  SVector{4}([0.0455085138515946 0.452461332332697 0.0455085138515946 0.456521639964114] ),
  SVector{4}([0.0455085138515946 0.452461332332697 0.456521639964114 0.0455085138515946] ),
  SVector{4}([0.456521639964114 0.0455085138515946 0.0455085138515946 0.452461332332697] ),
  SVector{4}([0.456521639964114 0.0455085138515946 0.452461332332697 0.0455085138515946] ),
  SVector{4}([0.456521639964114 0.452461332332697 0.0455085138515946 0.0455085138515946] ),
  SVector{4}([0.452461332332697 0.0455085138515946 0.0455085138515946 0.456521639964114] ),
  SVector{4}([0.452461332332697 0.0455085138515946 0.456521639964114 0.0455085138515946] ),
  SVector{4}([0.452461332332697 0.456521639964114 0.0455085138515946 0.0455085138515946] )
]
const tetrahedron20wn = [
  0.07348669573167155,
  0.07348669573167155,
  0.07348669573167155,
  0.07348669573167155,
  0.1126851232260692,
  0.1126851232260692,
  0.1126851232260692,
  0.1126851232260692,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975,
  0.02127606034741975
]

#x = [0.31570114977769714, 0.4495101774013808, 0.021265472540388547, 0.8108302411014925, 0.18883383102391316, 0.0471607003604522]
#w = [0.09548528946195194, 0.04232958121262396, 0.03189692783318026, 0.008110770829544959, 0.037207130727827595]
const tetrahedron35an = [
  SVector{4}([0.2500000000000000 0.2500000000000000 0.2500000000000000 0.250000000000000]),
  SVector{4}([0.315701149777697 0.315701149777697 0.315701149777697 0.0528965506669086] ),
  SVector{4}([0.315701149777697 0.315701149777697 0.0528965506669086 0.315701149777697] ),
  SVector{4}([0.315701149777697 0.0528965506669086 0.315701149777697 0.315701149777697] ),
  SVector{4}([0.0528965506669086 0.315701149777697 0.315701149777697 0.315701149777697] ),
  SVector{4}([0.449510177401381 0.449510177401381 0.0504898225986192 0.0504898225986192]),
  SVector{4}([0.449510177401381 0.0504898225986192 0.449510177401381 0.0504898225986192]),
  SVector{4}([0.449510177401381 0.0504898225986192 0.0504898225986192 0.449510177401381]),
  SVector{4}([0.0504898225986192 0.449510177401381 0.449510177401381 0.0504898225986192]),
  SVector{4}([0.0504898225986192 0.449510177401381 0.0504898225986192 0.449510177401381]),
  SVector{4}([0.0504898225986192 0.0504898225986192 0.449510177401381 0.449510177401381]),
  SVector{4}([0.0212654725403885 0.0212654725403885 0.810830241101492 0.146638813817730]),
  SVector{4}([0.0212654725403885 0.0212654725403885 0.146638813817730 0.810830241101492]),
  SVector{4}([0.0212654725403885 0.810830241101492 0.0212654725403885 0.146638813817730]),
  SVector{4}([0.0212654725403885 0.810830241101492 0.146638813817730 0.0212654725403885]),
  SVector{4}([0.0212654725403885 0.146638813817730 0.0212654725403885 0.810830241101492]),
  SVector{4}([0.0212654725403885 0.146638813817730 0.810830241101492 0.0212654725403885]),
  SVector{4}([0.810830241101492 0.0212654725403885 0.0212654725403885 0.146638813817730]),
  SVector{4}([0.810830241101492 0.0212654725403885 0.146638813817730 0.0212654725403885]),
  SVector{4}([0.810830241101492 0.146638813817730 0.0212654725403885 0.0212654725403885]),
  SVector{4}([0.146638813817730 0.0212654725403885 0.0212654725403885 0.810830241101492]),
  SVector{4}([0.146638813817730 0.0212654725403885 0.810830241101492 0.0212654725403885]),
  SVector{4}([0.146638813817730 0.810830241101492 0.0212654725403885 0.0212654725403885]),
  SVector{4}([0.188833831023913 0.188833831023913 0.0471607003604522 0.575171637591722] ),
  SVector{4}([0.188833831023913 0.188833831023913 0.575171637591722 0.0471607003604522] ),
  SVector{4}([0.188833831023913 0.0471607003604522 0.188833831023913 0.575171637591722] ),
  SVector{4}([0.188833831023913 0.0471607003604522 0.575171637591722 0.188833831023913] ),
  SVector{4}([0.188833831023913 0.575171637591722 0.188833831023913 0.0471607003604522] ),
  SVector{4}([0.188833831023913 0.575171637591722 0.0471607003604522 0.188833831023913] ),
  SVector{4}([0.0471607003604522 0.188833831023913 0.188833831023913 0.575171637591722] ),
  SVector{4}([0.0471607003604522 0.188833831023913 0.575171637591722 0.188833831023913] ),
  SVector{4}([0.0471607003604522 0.575171637591722 0.188833831023913 0.188833831023913] ),
  SVector{4}([0.575171637591722 0.188833831023913 0.188833831023913 0.0471607003604522] ),
  SVector{4}([0.575171637591722 0.188833831023913 0.0471607003604522 0.188833831023913] ),
  SVector{4}([0.575171637591722 0.0471607003604522 0.188833831023913 0.188833831023913] ),
]

const tetrahedron35wn = [
  0.09548528946195194,
  0.04232958121262396,
  0.04232958121262396,
  0.04232958121262396,
  0.04232958121262396,
  0.03189692783318026,
  0.03189692783318026,
  0.03189692783318026,
  0.03189692783318026,
  0.03189692783318026,
  0.03189692783318026,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.008110770829544959,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595,
  0.037207130727827595
]

#x = [0.01007724952363584, 0.08345441021380152, 0.03680559344436588, 0.1920920117640443, 0.0505176627051748, 0.45182539755827644, 0.22978489045560146, 0.027391908147254534, 0.22283977125262663, 0.12426032219090977]
#w = [ 0.0012530310990666355, 0.012033660585619839, 0.012884550846146703, 0.013643172516067632, 0.026693411503923556, 0.02568330123896662]

const tetrahedron56an = [
  SVector{4}([0.0100772495236358 0.0100772495236358 0.0100772495236358 0.969768251429093]),
  SVector{4}([0.0100772495236358 0.0100772495236358 0.969768251429093 0.0100772495236358]),
  SVector{4}([0.0100772495236358 0.969768251429093 0.0100772495236358 0.0100772495236358]),
  SVector{4}([0.969768251429093 0.0100772495236358 0.0100772495236358 0.0100772495236358]),
  SVector{4}([0.0834544102138015 0.0834544102138015 0.0834544102138015 0.749636769358595]),
  SVector{4}([0.0834544102138015 0.0834544102138015 0.749636769358595 0.0834544102138015]),
  SVector{4}([0.0834544102138015 0.749636769358595 0.0834544102138015 0.0834544102138015]),
  SVector{4}([0.749636769358595 0.0834544102138015 0.0834544102138015 0.0834544102138015]),
  SVector{4}([0.0368055934443659 0.0368055934443659 0.192092011764044 0.734296801347224] ),
  SVector{4}([0.0368055934443659 0.0368055934443659 0.734296801347224 0.192092011764044] ),
  SVector{4}([0.0368055934443659 0.192092011764044 0.0368055934443659 0.734296801347224] ),
  SVector{4}([0.0368055934443659 0.192092011764044 0.734296801347224 0.0368055934443659] ),
  SVector{4}([0.0368055934443659 0.734296801347224 0.0368055934443659 0.192092011764044] ),
  SVector{4}([0.0368055934443659 0.734296801347224 0.192092011764044 0.0368055934443659] ),
  SVector{4}([0.192092011764044 0.0368055934443659 0.0368055934443659 0.734296801347224] ),
  SVector{4}([0.192092011764044 0.0368055934443659 0.734296801347224 0.0368055934443659] ),
  SVector{4}([0.192092011764044 0.734296801347224 0.0368055934443659 0.0368055934443659] ),
  SVector{4}([0.734296801347224 0.0368055934443659 0.0368055934443659 0.192092011764044] ),
  SVector{4}([0.734296801347224 0.0368055934443659 0.192092011764044 0.0368055934443659] ),
  SVector{4}([0.734296801347224 0.192092011764044 0.0368055934443659 0.0368055934443659] ),
  SVector{4}([0.0505176627051748 0.0505176627051748 0.451825397558276 0.447139277031374] ),
  SVector{4}([0.0505176627051748 0.0505176627051748 0.447139277031374 0.451825397558276] ),
  SVector{4}([0.0505176627051748 0.451825397558276 0.0505176627051748 0.447139277031374] ),
  SVector{4}([0.0505176627051748 0.451825397558276 0.447139277031374 0.0505176627051748] ),
  SVector{4}([0.0505176627051748 0.447139277031374 0.0505176627051748 0.451825397558276] ),
  SVector{4}([0.0505176627051748 0.447139277031374 0.451825397558276 0.0505176627051748] ),
  SVector{4}([0.451825397558276 0.0505176627051748 0.0505176627051748 0.447139277031374] ),
  SVector{4}([0.451825397558276 0.0505176627051748 0.447139277031374 0.0505176627051748] ),
  SVector{4}([0.451825397558276 0.447139277031374 0.0505176627051748 0.0505176627051748] ),
  SVector{4}([0.447139277031374 0.0505176627051748 0.0505176627051748 0.451825397558276] ),
  SVector{4}([0.447139277031374 0.0505176627051748 0.451825397558276 0.0505176627051748] ),
  SVector{4}([0.447139277031374 0.451825397558276 0.0505176627051748 0.0505176627051748] ),
  SVector{4}([0.229784890455601 0.229784890455601 0.0273919081472545 0.513038310941543]  ),
  SVector{4}([0.229784890455601 0.229784890455601 0.513038310941543 0.0273919081472545]  ),
  SVector{4}([0.229784890455601 0.0273919081472545 0.229784890455601 0.513038310941543]  ),
  SVector{4}([0.229784890455601 0.0273919081472545 0.513038310941543 0.229784890455601]  ),
  SVector{4}([0.229784890455601 0.513038310941543 0.229784890455601 0.0273919081472545]  ),
  SVector{4}([0.229784890455601 0.513038310941543 0.0273919081472545 0.229784890455601]  ),
  SVector{4}([0.0273919081472545 0.229784890455601 0.229784890455601 0.513038310941543]  ),
  SVector{4}([0.0273919081472545 0.229784890455601 0.513038310941543 0.229784890455601]  ),
  SVector{4}([0.0273919081472545 0.513038310941543 0.229784890455601 0.229784890455601]  ),
  SVector{4}([0.513038310941543 0.229784890455601 0.229784890455601 0.0273919081472545]  ),
  SVector{4}([0.513038310941543 0.229784890455601 0.0273919081472545 0.229784890455601]  ),
  SVector{4}([0.513038310941543 0.0273919081472545 0.229784890455601 0.229784890455601]  ),
  SVector{4}([0.222839771252627 0.222839771252627 0.124260322190910 0.430060135303837]   ),
  SVector{4}([0.222839771252627 0.222839771252627 0.430060135303837 0.124260322190910]   ),
  SVector{4}([0.222839771252627 0.124260322190910 0.222839771252627 0.430060135303837]   ),
  SVector{4}([0.222839771252627 0.124260322190910 0.430060135303837 0.222839771252627]   ),
  SVector{4}([0.222839771252627 0.430060135303837 0.222839771252627 0.124260322190910]   ),
  SVector{4}([0.222839771252627 0.430060135303837 0.124260322190910 0.222839771252627]   ),
  SVector{4}([0.124260322190910 0.222839771252627 0.222839771252627 0.430060135303837]   ),
  SVector{4}([0.124260322190910 0.222839771252627 0.430060135303837 0.222839771252627]   ),
  SVector{4}([0.124260322190910 0.430060135303837 0.222839771252627 0.222839771252627]   ),
  SVector{4}([0.430060135303837 0.222839771252627 0.222839771252627 0.124260322190910]   ),
  SVector{4}([0.430060135303837 0.222839771252627 0.124260322190910 0.222839771252627]   ),
  SVector{4}([0.430060135303837 0.124260322190910 0.222839771252627 0.222839771252627]   )
]

const tetrahedron56wn = [
  0.0012530310990666355,
  0.0012530310990666355,
  0.0012530310990666355,
  0.0012530310990666355,
  0.012033660585619839,
  0.012033660585619839,
  0.012033660585619839,
  0.012033660585619839,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.012884550846146703,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.013643172516067632,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.026693411503923556,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662,
  0.02568330123896662
]

const tetrahedron84an = [
    SVector{4}([0.919364576755549   0.026878474414817   0.026878474414817   0.026878474414817]),
    SVector{4}([0.026878474414817   0.919364576755549   0.026878474414817   0.026878474414817]),
    SVector{4}([0.026878474414817   0.026878474414817   0.919364576755549   0.026878474414817]),
    SVector{4}([0.026878474414817   0.026878474414817   0.026878474414817   0.919364576755549]),
    SVector{4}([0.438577972589591   0.187140675803470   0.187140675803470   0.187140675803470]),
    SVector{4}([0.187140675803470   0.438577972589591   0.187140675803470   0.187140675803470]),
    SVector{4}([0.187140675803470   0.187140675803470   0.438577972589591   0.187140675803470]),
    SVector{4}([0.187140675803470   0.187140675803470   0.187140675803470   0.438577972589591]),
    SVector{4}([0.473575835127937   0.473575835127937   0.026424164872063   0.026424164872063]),
    SVector{4}([0.473575835127937   0.026424164872063   0.473575835127937   0.026424164872063]),
    SVector{4}([0.473575835127937   0.026424164872063   0.026424164872063   0.473575835127937]),
    SVector{4}([0.026424164872063   0.473575835127937   0.473575835127937   0.026424164872063]),
    SVector{4}([0.026424164872063   0.473575835127937   0.026424164872063   0.473575835127937]),
    SVector{4}([0.026424164872063   0.026424164872063   0.473575835127937   0.473575835127937]),
    SVector{4}([0.352045262027356   0.352045262027356   0.147954737972644   0.147954737972644]),
    SVector{4}([0.352045262027356   0.147954737972644   0.352045262027356   0.147954737972644]),
    SVector{4}([0.352045262027356   0.147954737972644   0.147954737972644   0.352045262027356]),
    SVector{4}([0.147954737972644   0.352045262027356   0.352045262027356   0.147954737972644]),
    SVector{4}([0.147954737972644   0.352045262027356   0.147954737972644   0.352045262027356]),
    SVector{4}([0.147954737972644   0.147954737972644   0.352045262027356   0.352045262027356]),
    SVector{4}([0.225783205866940   0.732309909692947   0.020953442220056   0.020953442220056]),
    SVector{4}([0.225783205866940   0.020953442220056   0.732309909692947   0.020953442220056]),
    SVector{4}([0.225783205866940   0.020953442220056   0.020953442220056   0.732309909692947]),
    SVector{4}([0.020953442220056   0.225783205866940   0.732309909692947   0.020953442220056]),
    SVector{4}([0.020953442220056   0.225783205866940   0.020953442220056   0.732309909692947]),
    SVector{4}([0.020953442220056   0.020953442220056   0.225783205866940   0.732309909692947]),
    SVector{4}([0.732309909692947   0.225783205866940   0.020953442220056   0.020953442220056]),
    SVector{4}([0.732309909692947   0.020953442220056   0.225783205866940   0.020953442220056]),
    SVector{4}([0.732309909692947   0.020953442220056   0.020953442220056   0.225783205866940]),
    SVector{4}([0.020953442220056   0.732309909692947   0.225783205866940   0.020953442220056]),
    SVector{4}([0.020953442220056   0.732309909692947   0.020953442220056   0.225783205866940]),
    SVector{4}([0.020953442220056   0.020953442220056   0.732309909692947   0.225783205866940]),
    SVector{4}([0.158462939666092   0.647557594086975   0.096989733123466   0.096989733123466]),
    SVector{4}([0.158462939666092   0.096989733123466   0.647557594086975   0.096989733123466]),
    SVector{4}([0.158462939666092   0.096989733123466   0.096989733123466   0.647557594086975]),
    SVector{4}([0.096989733123466   0.158462939666092   0.647557594086975   0.096989733123466]),
    SVector{4}([0.096989733123466   0.158462939666092   0.096989733123466   0.647557594086975]),
    SVector{4}([0.096989733123466   0.096989733123466   0.158462939666092   0.647557594086975]),
    SVector{4}([0.647557594086975   0.158462939666092   0.096989733123466   0.096989733123466]),
    SVector{4}([0.647557594086975   0.096989733123466   0.158462939666092   0.096989733123466]),
    SVector{4}([0.647557594086975   0.096989733123466   0.096989733123466   0.158462939666092]),
    SVector{4}([0.096989733123466   0.647557594086975   0.158462939666092   0.096989733123466]),
    SVector{4}([0.096989733123466   0.647557594086975   0.096989733123466   0.158462939666092]),
    SVector{4}([0.096989733123466   0.096989733123466   0.647557594086975   0.158462939666092]),
    SVector{4}([0.322111431830857   0.322111431830857   0.322111431830857   0.033665704507429]),
    SVector{4}([0.322111431830857   0.322111431830857   0.033665704507429   0.322111431830857]),
    SVector{4}([0.322111431830857   0.033665704507429   0.322111431830857   0.322111431830857]),
    SVector{4}([0.033665704507429   0.322111431830857   0.322111431830857   0.322111431830857]),
    SVector{4}([0.792939256469618   0.097608162890442   0.097608162890442   0.011844417749498]),
    SVector{4}([0.097608162890442   0.792939256469618   0.097608162890442   0.011844417749498]),
    SVector{4}([0.097608162890442   0.097608162890442   0.792939256469618   0.011844417749498]),
    SVector{4}([0.792939256469618   0.097608162890442   0.011844417749498   0.097608162890442]),
    SVector{4}([0.097608162890442   0.792939256469618   0.011844417749498   0.097608162890442]),
    SVector{4}([0.097608162890442   0.097608162890442   0.011844417749498   0.792939256469618]),
    SVector{4}([0.792939256469618   0.011844417749498   0.097608162890442   0.097608162890442]),
    SVector{4}([0.097608162890442   0.011844417749498   0.792939256469618   0.097608162890442]),
    SVector{4}([0.097608162890442   0.011844417749498   0.097608162890442   0.792939256469618]),
    SVector{4}([0.011844417749498   0.792939256469618   0.097608162890442   0.097608162890442]),
    SVector{4}([0.011844417749498   0.097608162890442   0.792939256469618   0.097608162890442]),
    SVector{4}([0.011844417749498   0.097608162890442   0.097608162890442   0.792939256469618]),
    SVector{4}([0.541184412800237   0.133558160703568   0.296501020543124   0.028756405953071]),
    SVector{4}([0.296501020543124   0.541184412800237   0.133558160703568   0.028756405953071]),
    SVector{4}([0.133558160703568   0.296501020543124   0.541184412800237   0.028756405953071]),
    SVector{4}([0.541184412800237   0.296501020543124   0.133558160703568   0.028756405953071]),
    SVector{4}([0.133558160703568   0.541184412800237   0.296501020543124   0.028756405953071]),
    SVector{4}([0.296501020543124   0.133558160703568   0.541184412800237   0.028756405953071]),
    SVector{4}([0.541184412800237   0.133558160703568   0.028756405953071   0.296501020543124]),
    SVector{4}([0.296501020543124   0.541184412800237   0.028756405953071   0.133558160703568]),
    SVector{4}([0.133558160703568   0.296501020543124   0.028756405953071   0.541184412800237]),
    SVector{4}([0.541184412800237   0.296501020543124   0.028756405953071   0.133558160703568]),
    SVector{4}([0.133558160703568   0.541184412800237   0.028756405953071   0.296501020543124]),
    SVector{4}([0.296501020543124   0.133558160703568   0.028756405953071   0.541184412800237]),
    SVector{4}([0.541184412800237   0.028756405953071   0.133558160703568   0.296501020543124]),
    SVector{4}([0.296501020543124   0.028756405953071   0.541184412800237   0.133558160703568]),
    SVector{4}([0.133558160703568   0.028756405953071   0.296501020543124   0.541184412800237]),
    SVector{4}([0.541184412800237   0.028756405953071   0.296501020543124   0.133558160703568]),
    SVector{4}([0.133558160703568   0.028756405953071   0.541184412800237   0.296501020543124]),
    SVector{4}([0.296501020543124   0.028756405953071   0.133558160703568   0.541184412800237]),
    SVector{4}([0.028756405953071   0.541184412800237   0.133558160703568   0.296501020543124]),
    SVector{4}([0.028756405953071   0.296501020543124   0.541184412800237   0.133558160703568]),
    SVector{4}([0.028756405953071   0.133558160703568   0.296501020543124   0.541184412800237]),
    SVector{4}([0.028756405953071   0.541184412800237   0.296501020543124   0.133558160703568]),
    SVector{4}([0.028756405953071   0.133558160703568   0.541184412800237   0.296501020543124]),
    SVector{4}([0.028756405953071   0.296501020543124   0.133558160703568   0.541184412800237])
]

const tetrahedron84wn = [
    0.002144935144316,
    0.002144935144316,
    0.002144935144316,
    0.002144935144316,
    0.020826641690769,
    0.020826641690769,
    0.020826641690769,
    0.020826641690769,
    0.007210136064455,
    0.007210136064455,
    0.007210136064455,
    0.007210136064455,
    0.007210136064455,
    0.007210136064455,
    0.030798919159712,
    0.030798919159712,
    0.030798919159712,
    0.030798919159712,
    0.030798919159712,
    0.030798919159712,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.004357844813864,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.008593530677833,
    0.023000681669286,
    0.023000681669286,
    0.023000681669286,
    0.023000681669286,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.004863063904912,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259,
    0.015595140078259
]

const tetrahedronAn = [
  tetrahedron1an,
  tetrahedron4an,
  tetrahedron10an,
  tetrahedron20an,
  tetrahedron35an,
  tetrahedron56an,
  tetrahedron84an
]

const tetrahedronWn = [
  tetrahedron1wn,
  tetrahedron4wn,
  tetrahedron10wn,
  tetrahedron20wn,
  tetrahedron35wn,
  tetrahedron56wn,
  tetrahedron84wn
]


text:
module Interface

import Distributions: sample, Sampleable
import Random: GLOBAL_RNG, AbstractRNG, seed!
import MCMCChains: Chains, chainscat
import ProgressMeter

export AbstractSampler,
       AbstractTransition,
       AbstractCallback,
       init_callback,
       callback,
       transitions_init,
       transition_type,
       sample_init!,
       sample_end!,
       sample,
       psample,
       Sampleable,
       AbstractRNG,
       Chains,
       step!

"""
    AbstractSampler

The `AbstractSampler` type is intended to be inherited from when
implementing a custom sampler. Any persistent state information should be
saved in a subtype of `AbstractSampler`.

When defining a new sampler, you should also overload the function
`transition_type`, which tells the `sample` function what type of parameter
it should expect to receive.
"""
abstract type AbstractSampler end

"""
    AbstractTransition

The `AbstractTransition` type describes the results of a single step
of a given sampler. As an example, one implementation of an
`AbstractTransition` might include be a vector of parameters sampled from
a prior distribution.

Transition types should store a single draw from any sampler, since the
interface will sample `N` times, and store the results of each step in an
array of type `Array{Transition<:AbstractTransition, 1}`. If you were
using a sampler that returned a `NamedTuple` after each step, your
implementation might look like:

```
struct MyTransition <: AbstractTransition
    draw :: NamedTuple
end
```
"""
abstract type AbstractTransition end

"""
    AbstractCallback

An `AbstractCallback` types is a supertype to be inherited from if you want to use custom callback
functionality. This is used to report sampling progress such as parameters calculated, remaining
samples to run, or even plot graphs if you so choose.

In order to implement callback functionality, you need the following:

- A mutable struct that is a subtype of `AbstractCallback`
- An overload of the `init_callback` function
- An overload of the `callback` function
"""
abstract type AbstractCallback end

"""
    NoCallback()

This disables the callback functionality in the event that you wish to
implement your own callback or reporting.
"""
mutable struct NoCallback <: AbstractCallback end

"""
    DefaultCallback(N::Int)

The default callback struct which uses `ProgressMeter`.
"""
mutable struct DefaultCallback{
    ProgType<:ProgressMeter.AbstractProgress
} <: AbstractCallback
    p :: ProgType
end

DefaultCallback(N::Int) = DefaultCallback(ProgressMeter.Progress(N, 1))

function init_callback(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    return DefaultCallback(N)
end

"""
    _generate_callback(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        progress_style=:default,
        kwargs...
    )

`_generate_callback` uses a `progress_style` keyword argument to determine
which progress meter style should be used. This function is strictly internal
and is not meant to be overloaded. If you intend to add a custom `AbstractCallback`,
you should overload `init_callback` instead.

Options for `progress_style` include:

    - `:default` which returns the result of `init_callback`
    - `false` or `:disable` which returns a `NoCallback`
    - `:plain` which returns the default, simple `DefaultCallback`.
"""
function _generate_callback(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    progress_style=:default,
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    if progress_style == :default
        return init_callback(rng, ℓ, s, N; kwargs...)
    elseif progress_style == false || progress_style == :disable
        return NoCallback()
    elseif progress_style == :plain
        return DefaultCallback(N)
    else
        throw(ArgumentError("Keyword argument $progress_style is not recognized."))
    end
end

"""
    sample(
        rng::AbstractRNG,
        ℓ::Sampleable,
        s::AbstractSampler,
        N::Integer;
        kwargs...
    )

    sample(
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

`sample` returns an `MCMCChains.Chains` object containing `N` samples from a given model and
sampler. You may pass in any additional arguments through the use of keyword arguments.
"""

function sample(
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    return sample(GLOBAL_RNG, ℓ, s, N; kwargs...)
end

function sample(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    progress::Bool=true,
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    # Perform any necessary setup.
    sample_init!(rng, ℓ, s, N; kwargs...)

    # Preallocate the TransitionType vector.
    ts = transitions_init(rng, ℓ, s, N; kwargs...)

    # Add a progress meter.
    cb = progress ? _generate_callback(rng, ℓ, s, N; kwargs...) : nothing

    # Step through the sampler.
    for i=1:N
        if i == 1
            ts[i] = step!(rng, ℓ, s, N; iteration=i, kwargs...)
        else
            ts[i] = step!(rng, ℓ, s, N, ts[i-1]; iteration=i, kwargs...)
        end

        # Run a callback function.
        progress && callback(rng, ℓ, s, N, i, ts[i], cb; kwargs...)
    end

    # Wrap up the sampler, if necessary.
    sample_end!(rng, ℓ, s, N, ts; kwargs...)

    return Chains(rng, ℓ, s, N, ts; kwargs...)
end

"""
    sample_init!(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

Performs whatever initial setup is required for your sampler. This function is not intended
to return any value -- any set up should mutate the sampler or the model type in-place.

A common use for `sample_init!` might be to instantiate a particle field for later use,
or find an initial step size for a Hamiltonian sampler.
"""
function sample_init!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    debug::Bool=false,
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    # Do nothing.
    debug && @warn "No sample_init! function has been implemented for objects
           of types $(typeof(ℓ)) and $(typeof(s))"
end

"""
    sample_end!(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer,
        ts::Vector{TransitionType};
        kwargs...
    )

Performs whatever finalizing the sampler requires. This function is not intended
to return any value -- any set up should mutate the sampler or the model type in-place.

`sample_end!` is useful in cases where you might like to perform some transformation
on your vector of `AbstractTransitions`, save your sampler struct to disk, or otherwise
perform any clean-up or finalization.
"""
function sample_end!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    ts::Vector{TransitionType};
    debug::Bool=false,
    kwargs...
) where {
    ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    TransitionType<:AbstractTransition
}
    # Do nothing.
    debug && @warn "No sample_end! function has been implemented for objects
           of types $(typeof(ℓ)) and $(typeof(s))"
end

"""
    step!(
        rng::AbstractRNG,
        ℓ::Sampleable,
        s::AbstractSampler,
        N::Integer;
        kwargs...
    )

    step!(
        rng::AbstractRNG,
        ℓ::Sampleable,
        s::AbstractSampler;
        kwargs...
    )

    step!(
        rng::AbstractRNG,
        ℓ::Sampleable,
        s::AbstractSampler,
        N::Integer,
        t::AbstractTransition;
        kwargs...
    )

Returns a single `AbstractTransition` drawn using the provided random number generator,
model, and sampler. `step!` is the function that performs inference, and it is how
a model moves from one sample to another.

`step!` may modify the model or the sampler in-place. As an example, you may have a state
variable in your sampler that contains a vector of particles or some other value that
does not need to be included in the `AbstractTransition` struct returned.

Every `step!` call after the first has access to the previous `AbstractTransition`.
"""
function step!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    debug::Bool=false,
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    # Do nothing.
    debug && @warn "No step! function has been implemented for objects of types \n- $(typeof(ℓ)) \n- $(typeof(s))"
end

function step!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    return step!(rng, ℓ, s, 1; kwargs...)
end

function step!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    t::TransitionType;
    kwargs...
) where {ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    TransitionType<:AbstractTransition
}
    # Do nothing.
    # @warn "No step! function has been implemented for objects
    #        of types $(typeof(ℓ)) and $(typeof(s))"
    return step!(rng, ℓ, s, N; kwargs...)
end

function step!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    t::Nothing;
    debug::Bool=true,
    kwargs...
) where {ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    TransitionType<:AbstractTransition
}
    debug && @warn "No transition type passed in, running normal step! function."
    return step!(rng, ℓ, s, N; kwargs...)
end

"""
    transitions_init(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

Generates a vector of `AbstractTransition` types of length `N`.
"""
function transitions_init(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    return Vector{transition_type(s)}(undef, N)
end

"""
    callback(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer,
        iteration::Integer,
        cb::CallbackType;
        kwargs...
    )

`callback` is called after every sample run, and allows you to run some function on a
subtype of `AbstractCallback`. Typically this is used to increment a progress meter, show a
plot of parameter draws, or otherwise provide information about the sampling process to the user.

By default, `ProgressMeter` is used to show the number of samples remaning.
"""
function callback(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    iteration::Integer,
    t::TransitionType,
    cb::CallbackType;
    kwargs...
) where {
    ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    CallbackType<:AbstractCallback,
    TransitionType<:AbstractTransition
}
    # Default callback behavior.
    ProgressMeter.next!(cb.p)
end

function callback(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    iteration::Integer,
    t::TransitionType,
    cb::NoCallback;
    kwargs...
) where {
    ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    TransitionType<:AbstractTransition
}
    # Do nothing.
end

"""
    transition_type(s::AbstractSampler)

Return the type of `AbstractTransition` that is to be returned by an
`AbstractSampler` after each `step!` call.
"""
transition_type(s::AbstractSampler) = AbstractTransition

"""
    psample(ℓ::ModelType, s::SamplerType, N::Integer, n_chains::Integer; kwargs...)
    psample(rng::AbstractRNG, ℓ::ModelType, s::SamplerType, N::Integer, n_chains::Integer; kwargs...)

Samples `n_chains` using the threads available, and combines the returned chains into a single
`Chains` struct.

By default, the model and sampler types are deep copied for each thread to prevent contamination
between threads.
"""
function psample(
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    n_chains::Integer;
    kwargs...
) where {
    ModelType <: Sampleable,
    SamplerType <: AbstractSampler
}
    return psample(GLOBAL_RNG, ℓ, s, N, n_chains; kwargs...)
end

function psample(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    n_chains::Integer;
    kwargs...
) where {
    ModelType <: Sampleable,
    SamplerType <: AbstractSampler
}
    # Set up a chains vector.
    chains = []

    # Create a thread lock. Used to add append chains to the vector.
    c_lock = Threads.SpinLock()

    # Create a seed for each chain using rng.
    seeds = rand(rng, UInt, n_chains)

    Threads.@threads for i in 1:n_chains
        # Generate a new RNG using the pre-made seeds.
        sub_rng = seed!(seeds[i])

        # Sample a chain.
        chain = sample(sub_rng, deepcopy(ℓ), deepcopy(s), N; progress=false, kwargs...)

        # Acquire the lock, add the chain to the vector.
        lock(c_lock) do
            push!(chains, (order=i, chain=chain))
        end
    end

    # Sort the chains by seed order, to ensure that all chains are always in the same order.
    sort!(chains, by=x -> x.order)

    # Concatenate the chains together.
    return chainscat([c.chain for c in chains]...)
end

end # module Interface


text:
src = [rand(10^4) for i in 1:10^4]

@testset "sort_batch" begin
    println("--------")
    println("tmap!")
    dst = similar(src)
    @time tmap!(sort, dst, src)
    @test all(issorted.(dst))

    println("simple")
    function simple(src, dst)
        @Threads.threads for k in 1:length(src)
            dst[k] = sort(src[k])
        end
    end

    dst = similar(src)
    @time simple(src, dst)
    @test all(issorted.(dst))
end


text:
const LIMITERS = (:VanLeer,)

include(joinpath(dirname(@__FILE__),"..","..","src","VlasovSolver","Limiters.jl"))

function test_limiter(limiter_name)
    limiter = Limiters.generate_limiter(limiter_name)

    # test special points
    @test limiter(0.0) ≈ 0.0
    @test limiter(1.0) ≈ 1.0

    # test symmetricity
    for r in (2.0, 3.0, 5.0, 11.0)
        @test limiter(r) ≈ r*limiter(1/r)
    end

    # test TVD condition
    for r in range(0.1, 0.9; step=0.1)
        @test r ≤ limiter(r) ≤ 2*r
    end
    for r in range(1.1, 1.9; step=0.1)
        @test 1 ≤ limiter(r) ≤ r
    end
    for r in range(2.0, 10.0; step=1.0)
        @test 1.0 ≤ limiter(r) ≤ 2.0
    end
end

for limiter_name in LIMITERS
    test_limiter(limiter_name)
end


text:
module FaceLabelsTests

using Test
using Gridap

vertex_to_geolabel = [1,1,2,2,2,1,1,3,3]
edge_to_geolabel = [4,4,5,5,5,5,6,6,4]
physlabel_1 = [1,3,4]
physlabel_2 = [5,3,6,2]
tag_to_name = ["label1","label2"]

labels = FaceLabels(
  [vertex_to_geolabel, edge_to_geolabel],
  [physlabel_1, physlabel_2],
  tag_to_name)

@test isa(labels,FaceLabels)
@test labels_on_dim(labels,0) == vertex_to_geolabel
@test labels_on_dim(labels,1) == edge_to_geolabel
@test labels_on_tag(labels,1) == physlabel_1
@test labels_on_tag(labels,2) == physlabel_2
@test tag_from_name(labels,"label1")==1

r = [1, 1, 2, 2, 2, 2, 2, 2, 1]
@test collect(first_tag_on_face(labels,1)) == r

add_tag_from_tags!(labels,"label3",[1,2])

@test labels_on_tag(labels,3) == [1, 2, 3, 4, 5, 6]
@test tag_from_name(labels,"label3") == 3

add_tag_from_tags!(labels,"label4",["label1","label2"])

@test labels_on_tag(labels,4) == [1, 2, 3, 4, 5, 6]
@test tag_from_name(labels,"label4") == 4

@test string(labels) == "FaceLabels object"
s = "FaceLabels object:\n 0-faces: 9\n 1-faces: 9\n tags: 4"
@test sprint(show,"text/plain",labels) == s

end # module


text:
#=
You are given an array of non-negative integers that represents a two-dimensional elevation map where each element is unit-width wall and the integer is the height. Suppose it will rain and all spots between two walls get filled up.

Compute how many units of water remain trapped on the map in O(N) time and O(1) space.

For example, given the input [2, 1, 2], we can hold 1 unit of water in the middle.

Given the input [3, 0, 1, 3, 0, 5], we can hold 3 units in the first index, 2 in the second, and 3 in the fourth index (we cannot hold 5 since it would run off to the left), so we can trap 8 units of water.
=#
using Test
include("Solutions/problem30_watertrapped_in_elevation_map.jl")

@test watertrapped([2, 1, 2]) == 1
@test watertrapped([3, 0, 1, 3, 0, 5]) == 8

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
abstract type AbstractSymmetry <: AbstractShape end

abstract type AbstractOctantSymmetry <: AbstractSymmetry end

abstract type AbstractOddOctantSymmetry <: AbstractOctantSymmetry end

abstract type AbstractEvenOctantSymmetry <: AbstractOctantSymmetry end

struct OddOctantSymmetricShape{I <: Integer, S <: AbstractShape} <: AbstractOddOctantSymmetry
    origin::Point{I}
    shape::S
end

struct EvenOctantSymmetricShape{I <: Integer, S <: AbstractShape} <: AbstractOddOctantSymmetry
    origin::Point{I}
    shape::S
end

struct OddSymmetricVerticalLines4{I <: Integer} <: AbstractSymmetry
    origin::Point{I}
    point::Point{I}
end

struct EvenSymmetricVerticalLines4{I <: Integer} <: AbstractSymmetry
    origin::Point{I}
    point::Point{I}
end

#####
##### reflections
#####

get_odd_reflection_horizontal_line(shape::Point, i) = Point(i + i - shape.i, shape.j)
get_even_reflection_horizontal_line(shape::Point, i) = Point(i + i - shape.i - one(i), shape.j)

get_odd_reflection_vertical_line(shape::Point, j) = Point(shape.i, j + j - shape.j)
get_even_reflection_vertical_line(shape::Point, j) = Point(shape.i, j + j - shape.j - one(j))

get_reflection_diagonal_line(shape::Point, i, j) = Point(i + shape.j - j, j + shape.i - i)

get_odd_reflection_horizontal_line(shape::VerticalLine, i) = VerticalLine(i + i - shape.i_max, i + i - shape.i_min, shape.j)
get_even_reflection_horizontal_line(shape::VerticalLine, i) = VerticalLine(i + i - shape.i_max - one(i), i + i - shape.i_min - one(i), shape.j)

get_odd_reflection_vertical_line(shape::VerticalLine, j) = VerticalLine(shape.i_min, shape.i_max, j + j - shape.j)
get_even_reflection_vertical_line(shape::VerticalLine, j) = VerticalLine(shape.i_min, shape.i_max, j + j - shape.j - one(j))

function get_reflection_diagonal_line(shape::VerticalLine, i, j)
    new_i = i + shape.j - j
    j_min = j + shape.i_min - i
    j_max = j + shape.i_max - i
    return HorizontalLine(new_i, j_min, j_max)
end

#####
##### AbstractOctantSymmetry
#####

is_valid(shape::AbstractOctantSymmetry) = is_valid(shape.shape)

#####
##### OddOctantSymmetricShape
#####

function get_i_min(shape::OddOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_odd_reflection_horizontal_line(generator, origin.i)
    return min(get_i_min(generator), get_i_min(generator_reflected))
end

function get_i_max(shape::OddOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_odd_reflection_horizontal_line(generator, origin.i)
    return max(get_i_max(generator), get_i_max(generator_reflected))
end

function get_i_extrema(shape::OddOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_odd_reflection_horizontal_line(generator, origin.i)
    return (min(get_i_min(generator), get_i_min(generator_reflected)), max(get_i_max(generator), get_i_max(generator_reflected)))
end

function get_j_min(shape::OddOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_odd_reflection_horizontal_line(generator, origin.j)
    return min(get_j_min(generator), get_j_min(generator_reflected))
end

function get_j_max(shape::OddOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_odd_reflection_horizontal_line(generator, origin.j)
    return max(get_j_max(generator), get_j_max(generator_reflected))
end

function get_j_extrema(shape::OddOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_odd_reflection_horizontal_line(generator, origin.j)
    return (min(get_j_min(generator), get_j_min(generator_reflected)), max(get_j_max(generator), get_j_max(generator_reflected)))
end

function draw!(f::F, image::AbstractMatrix, shape::OddOctantSymmetricShape{I, Point{I}}, color) where {F <: Function, I}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    origin = shape.origin
    point = shape.shape

    i_origin = origin.i
    j_origin = origin.j
    i = point.i
    j = point.j

    i_diff = i - i_origin
    j_diff = j - j_origin

    f(image, i_origin - j_diff, j_origin - i_diff, color)
    f(image, i_origin + j_diff, j_origin - i_diff, color)
    f(image, i_origin - i_diff, j_origin - j_diff, color)
    f(image, i_origin + i_diff, j_origin - j_diff, color)
    f(image, i_origin - i_diff, j_origin + j_diff, color)
    f(image, i, j, color)
    f(image, i_origin - j_diff, j_origin + i_diff, color)
    f(image, i_origin + j_diff, j_origin + i_diff, color)

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::OddOctantSymmetricShape{I, VerticalLine{I}}, color) where {F <: Function, I}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    origin = shape.origin
    line = shape.shape

    i_origin = origin.i
    j_origin = origin.j

    i_min_line = line.i_min
    i_max_line = line.i_max
    j_line = line.j

    i_inner_relative = i_min_line - i_origin
    i_outer_relative = i_max_line - i_origin
    j_relative = j_line - j_origin

    draw!(f, image, HorizontalLine(i_origin - j_relative, j_origin - i_outer_relative, j_origin - i_inner_relative), color)
    draw!(f, image, HorizontalLine(i_origin + j_relative, j_origin - i_outer_relative, j_origin - i_inner_relative), color)
    draw!(f, image, VerticalLine(i_origin - i_outer_relative, i_origin - i_inner_relative, j_origin - j_relative), color)
    draw!(f, image, VerticalLine(i_origin + i_inner_relative, i_origin + i_outer_relative, j_origin - j_relative), color)
    draw!(f, image, VerticalLine(i_origin - i_outer_relative, i_origin - i_inner_relative, j_origin + j_relative), color)
    draw!(f, image, VerticalLine(i_origin + i_inner_relative, i_origin + i_outer_relative, j_origin + j_relative), color)
    draw!(f, image, HorizontalLine(i_origin - j_relative, j_origin + i_inner_relative, j_origin + i_outer_relative), color)
    draw!(f, image, HorizontalLine(i_origin + j_relative, j_origin + i_inner_relative, j_origin + i_outer_relative), color)

    return nothing
end

#####
##### EvenOctantSymmetricShape
#####

function get_i_min(shape::EvenOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_even_reflection_horizontal_line(generator, origin.i)
    return min(get_i_min(generator), get_i_min(generator_reflected))
end

function get_i_max(shape::EvenOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_even_reflection_horizontal_line(generator, origin.i)
    return max(get_i_max(generator), get_i_max(generator_reflected))
end

function get_i_extrema(shape::EvenOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_even_reflection_horizontal_line(generator, origin.i)
    return (min(get_i_min(generator), get_i_min(generator_reflected)), max(get_i_max(generator), get_i_max(generator_reflected)))
end

function get_j_min(shape::EvenOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_even_reflection_horizontal_line(generator, origin.j)
    return min(get_j_min(generator), get_j_min(generator_reflected))
end

function get_j_max(shape::EvenOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_even_reflection_horizontal_line(generator, origin.j)
    return max(get_j_max(generator), get_j_max(generator_reflected))
end

function get_j_extrema(shape::EvenOctantSymmetricShape)
    origin = shape.origin
    generator = shape.shape
    generator_reflected = get_even_reflection_horizontal_line(generator, origin.j)
    return (min(get_j_min(generator), get_j_min(generator_reflected)), max(get_j_max(generator), get_j_max(generator_reflected)))
end

function draw!(f::F, image::AbstractMatrix, shape::EvenOctantSymmetricShape{I, Point{I}}, color) where {F <: Function, I}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    origin = shape.origin
    point = shape.shape

    i_origin = origin.i
    j_origin = origin.j
    i = point.i
    j = point.j

    i_diff = i - i_origin
    j_diff = j - j_origin

    f(image, i_origin - j_diff, j_origin - i_diff, color)
    f(image, i_origin + j_diff - one(I), j_origin - i_diff, color)
    f(image, i_origin - i_diff, j_origin - j_diff, color)
    f(image, i_origin + i_diff - one(I), j_origin - j_diff, color)
    f(image, i_origin - i_diff, j_origin + j_diff - one(I), color)
    f(image, i - one(I), j - one(I), color)
    f(image, i_origin - j_diff, j_origin + i_diff - one(I), color)
    f(image, i_origin + j_diff - one(I), j_origin + i_diff - one(I), color)

    return nothing
end

function draw!(f::F, image::AbstractMatrix, shape::EvenOctantSymmetricShape{I, VerticalLine{I}}, color) where {F <: Function, I}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    origin = shape.origin
    line = shape.shape

    i_origin = origin.i
    j_origin = origin.j

    i_min_line = line.i_min
    i_max_line = line.i_max
    j_line = line.j

    i_inner_relative = i_min_line - i_origin
    i_outer_relative = i_max_line - i_origin
    j_relative = j_line - j_origin

    draw!(f, image, HorizontalLine(i_origin - j_relative, j_origin - i_outer_relative, j_origin - i_inner_relative), color)
    draw!(f, image, HorizontalLine(i_origin + j_relative - one(I), j_origin - i_outer_relative, j_origin - i_inner_relative), color)
    draw!(f, image, VerticalLine(i_origin - i_outer_relative, i_origin - i_inner_relative, j_origin - j_relative), color)
    draw!(f, image, VerticalLine(i_origin + i_inner_relative - one(I), i_origin + i_outer_relative - one(I), j_origin - j_relative), color)
    draw!(f, image, VerticalLine(i_origin - i_outer_relative, i_origin - i_inner_relative, j_origin + j_relative - one(I)), color)
    draw!(f, image, VerticalLine(i_origin + i_inner_relative - one(I), i_origin + i_outer_relative - one(I), j_origin + j_relative - one(I)), color)
    draw!(f, image, HorizontalLine(i_origin - j_relative, j_origin + i_inner_relative - one(I), j_origin + i_outer_relative - one(I)), color)
    draw!(f, image, HorizontalLine(i_origin + j_relative - one(I), j_origin + i_inner_relative - one(I), j_origin + i_outer_relative - one(I)), color)

    return nothing
end

#####
##### OddSymmetricVerticalLines4
#####

function is_valid(shape::OddSymmetricVerticalLines4)
    origin = shape.origin
    point = shape.point

    i_origin = origin.i
    j_origin = origin.j

    i_point = point.i
    j_point = point.j

    return (i_point >= i_origin) && (j_point >= j_origin) && (i_point - i_origin >= j_point - j_origin)
end

get_i_min(shape::OddSymmetricVerticalLines4) = get_odd_reflection_horizontal_line(shape.point, shape.origin.i).i
get_i_max(shape::OddSymmetricVerticalLines4) = shape.point.i

function get_j_min(shape::OddSymmetricVerticalLines4)
    diagonally_reflected_point = get_reflection_diagonal_line(shape.point, shape.origin.i, shape.origin.j)
    return get_odd_reflection_vertical_line(diagonally_reflected_point, shape.origin.j).j
end

get_j_max(shape::OddSymmetricVerticalLines4) = get_reflection_diagonal_line(shape.point, shape.origin.i, shape.origin.j).j

function get_j_extrema(shape::OddSymmetricVerticalLines4)
    diagonally_reflected_point = get_reflection_diagonal_line(shape.point, shape.origin.i, shape.origin.j)
    return (get_odd_reflection_vertical_line(diagonally_reflected_point, shape.origin.j).j, diagonally_reflected_point.j)
end

function draw!(f::F, image::AbstractMatrix, shape::OddSymmetricVerticalLines4, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    origin = shape.origin
    point = shape.point

    i_origin = origin.i
    j_origin = origin.j
    i = point.i
    j = point.j

    i_diff = i - i_origin
    j_diff = j - j_origin

    draw!(f, image, VerticalLine(i_origin - j_diff, i_origin + j_diff, j_origin - i_diff), color)
    draw!(f, image, VerticalLine(i_origin - i_diff, i_origin + i_diff, j_origin - j_diff), color)
    draw!(f, image, VerticalLine(i_origin - i_diff, i_origin + i_diff, j_origin + j_diff), color)
    draw!(f, image, VerticalLine(i_origin - j_diff, i_origin + j_diff, j_origin + i_diff), color)

    return nothing
end

#####
##### EvenSymmetricVerticalLines4
#####

function is_valid(shape::EvenSymmetricVerticalLines4)
    origin = shape.origin
    point = shape.point

    i_origin = origin.i
    j_origin = origin.j

    i_point = point.i
    j_point = point.j

    return (i_point > i_origin) && (j_point >= j_origin) && (i_point - i_origin >= j_point - j_origin)
end

get_i_min(shape::EvenSymmetricVerticalLines4) = get_even_reflection_horizontal_line(shape.point, shape.origin.i).i
get_i_max(shape::EvenSymmetricVerticalLines4) = shape.point.i

function get_j_min(shape::EvenSymmetricVerticalLines4)
    diagonally_reflected_point = get_reflection_diagonal_line(shape.point, shape.origin.i, shape.origin.j)
    return get_even_reflection_vertical_line(diagonally_reflected_point, shape.origin.j).j
end

get_j_max(shape::EvenSymmetricVerticalLines4) = get_reflection_diagonal_line(shape.point, shape.origin.i, shape.origin.j).j

function get_j_extrema(shape::EvenSymmetricVerticalLines4)
    diagonally_reflected_point = get_reflection_diagonal_line(shape.point, shape.origin.i, shape.origin.j)
    return (get_even_reflection_vertical_line(diagonally_reflected_point, shape.origin.j).j, diagonally_reflected_point.j)
end

function draw!(f::F, image::AbstractMatrix, shape::EvenSymmetricVerticalLines4, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    origin = shape.origin
    point = shape.point

    i_origin = origin.i
    j_origin = origin.j
    i = point.i
    j = point.j

    I = typeof(i_origin)

    i_diff = i - i_origin
    j_diff = j - j_origin

    draw!(f, image, VerticalLine(i_origin - j_diff, i_origin + j_diff - one(I), j_origin - i_diff), color)
    draw!(f, image, VerticalLine(i_origin - i_diff, i_origin + i_diff - one(I), j_origin - j_diff), color)
    draw!(f, image, VerticalLine(i_origin - i_diff, i_origin + i_diff - one(I), j_origin + j_diff - one(I)), color)
    draw!(f, image, VerticalLine(i_origin - j_diff, i_origin + j_diff - one(I), j_origin + i_diff - one(I)), color)

    return nothing
end


text:
#
# Layout for Expr nodes.
#

# Pattern matching for Expr objects.

# Supported syntax:
#   @isexpr ex _                => true
#   @isexpr ex var              => (local var = ex; true)
#   @isexpr ex :sym             => ex === :sym
#   @isexpr ex var := ...       => (local var = ex; ...)
#   @isexpr ex ::T              => ex isa T
#   @isexpr ex (... && ...)     => (... && ...)
#   @isexpr ex (... || ...)     => (... || ...)
#   @isexpr ex QuoteNode(...)   => (ex isa QuoteNode && (ex.value ...))
#   @isexpr ex Expr(...)        => (ex isa Expr && (ex.head ...) && (ex.args ...))

macro isexpr(val, pat)
    esc(_isexpr(val, pat))
end

function _isexpr(@nospecialize(val), @nospecialize(pat))
    pat !== :_ || return :(true)
    scr = gensym(:scr)
    ex = _isexpr(scr, pat)
    :(local $scr = $val; $ex)
end

function _isexpr(scr::Symbol, @nospecialize(pat))
    if pat isa Symbol
        if pat === :_
            :(true)
        else
            :(local $pat = $scr; true)
        end
    elseif pat isa QuoteNode && pat.value isa Symbol
        :($scr === $pat)
    elseif pat isa Expr
        nargs = length(pat.args)
        if pat.head === :(:=) && nargs == 2
            ex1 = _isexpr(scr, pat.args[1])
            ex2 = _isexpr(scr, pat.args[2])
            :($ex2 && $ex1)
        elseif pat.head === :(::) && nargs == 1
            T = pat.args[1]
            :($scr isa $T)
        elseif pat.head === :(::) && nargs == 2
            T = pat.args[2]
            ex = _isexpr(scr, pat.args[1])
            :($scr isa $T && $ex)
        elseif pat.head === :&& || pat.head === :||
            Expr(pat.head, Any[_isexpr(scr, argpat) for argpat in pat.args]...)
        elseif pat.head === :call && nargs >= 1
            call = pat.args[1]
            if call === :QuoteNode && nargs == 2
                ex = _isexpr(:($scr.value), pat.args[2])
                :($scr isa QuoteNode && $ex)
            elseif call === :Expr && nargs >= 2
                exs = Any[:($scr isa Expr)]
                push!(exs, _isexpr(:($scr.head), pat.args[2]))
                minlen = 0
                varlen = false
                for k = 3:lastindex(pat.args)
                    argpat = pat.args[k]
                    if argpat isa Expr && argpat.head === :... && length(argpat.args) == 1
                        !varlen || error("duplicate vararg patttern: $pat")
                        varlen = true
                    else
                        minlen += 1
                    end
                end
                if !varlen
                    push!(exs, :(length($scr.args) == $minlen))
                else
                    push!(exs, :(length($scr.args) >= $minlen))
                end
                j = 1
                seen_vararg = false
                for k = 3:lastindex(pat.args)
                    argpat = pat.args[k]
                    if argpat isa Expr && argpat.head === :... && length(argpat.args) == 1
                        argpat = argpat.args[1]
                        ex = _isexpr(:($scr.args[$j : $lastindex($scr.args) - $(minlen-j+1)]), argpat)
                        seen_vararg = true
                    elseif seen_vararg
                        ex = _isexpr(:($scr.args[$lastindex($scr.args) - $(minlen-j+1)]), argpat)
                    else
                        ex = _isexpr(:($scr.args[$j]), argpat)
                    end
                    push!(exs, ex)
                    j += 1
                end
                Expr(:&&, exs...)
            else
                error("invalid Expr pattern: $pat")
            end
        end
    else
        error("invalid Expr pattern: $pat")
    end
end

# Represent an object as an Expr node.

quoteof(obj) =
    obj

quoteof(ex::Union{Symbol,Expr,QuoteNode}) =
    QuoteNode(ex)

# Tiling Expr nodes.

tile_expr_or_repr(ex::Union{Symbol,Expr,QuoteNode}, pr = -1) =
    tile_expr(ex, pr)

tile_expr_or_repr(obj, pr = -1) =
    tile_repr(obj)

tile_expr(obj, pr = -1) =
    tile(obj)

tile_expr(qn::QuoteNode, pr = -1) =
    tile_expr_quoted(qn.value, pr)

identifier(s::Symbol) =
    Base.isidentifier(s) || Base.isoperator(s) || string(s)[1] == '@' ?
        string(s) :
        "var" * repr(string(s))

identifier(data) =
    repr(data)

tile_expr(s::Symbol, pr = -1) =
    literal(identifier(s))

function strip_line_number_nodes(ex::Expr)
    if any(arg -> arg isa LineNumberNode, ex.args)
        ex = Expr(ex.head, Any[arg for arg in ex.args if !(arg isa LineNumberNode)]...)
    end
    ex
end

function tile_expr(ex::Expr, pr = -1)
    ex = strip_line_number_nodes(ex)
    if @isexpr ex Expr(:block, args...)
        return tile_expr_block(args, pr)
    elseif @isexpr ex Expr(:module, notbare::Bool, name::Symbol, body)
        return tile_expr_module(notbare ? :module : :baremodule, name, body)
    elseif @isexpr ex Expr(:export, names...)
        return tile_expr_export(names)
    elseif @isexpr ex Expr(head := :using || :import, Expr(:(:), from, args...))
        return tile_expr_import(head, from, args)
    elseif @isexpr ex Expr(head := :using || :import, arg)
        return tile_expr_import(head, arg)
    elseif @isexpr ex Expr(:struct, mut::Bool, decl, body)
        return tile_expr_struct(mut ? Symbol("mutable struct") : :struct, decl, body)
    elseif @isexpr ex Expr(:abstract, decl)
        return tile_expr_abstract(decl)
    elseif @isexpr ex Expr(head := :if || :while || :for, cond, body)
        return tile_expr_stmt(head, cond, body, pr)
    elseif @isexpr ex Expr(head := :if, cond, body, else_body)
        return tile_expr_stmt(head, cond, body, else_body, pr)
    elseif @isexpr ex Expr(:let, decl, body)
        return tile_expr_let(decl, body)
    elseif @isexpr ex Expr(:try, body, ::Bool, ::Bool, finally_body)
        return tile_expr_try(body, nothing, nothing, finally_body)
    elseif @isexpr ex Expr(:try, body, ::Bool, catch_body, finally_body)
        return tile_expr_try(body, nothing, catch_body, finally_body)
    elseif @isexpr ex Expr(:try, body, name, catch_body, finally_body)
        return tile_expr_try(body, name, catch_body, finally_body)
    elseif @isexpr ex Expr(:try, body, ::Bool, catch_body)
        return tile_expr_try(body, nothing, catch_body, nothing)
    elseif @isexpr ex Expr(:try, body, name, catch_body)
        return tile_expr_try(body, name, catch_body, nothing)
    elseif @isexpr ex Expr(head := :function || :macro, decl, body)
        return tile_expr_function(head, decl, body)
    elseif (@isexpr ex Expr(:macrocall, ref::GlobalRef, docstr::String, def)) && ref === GlobalRef(Core, Symbol("@doc"))
        return tile_expr(docstr, pr) / tile_expr(def, pr)
    elseif @isexpr ex Expr(head := :break || :continue || :return)
        return literal(head)
    elseif @isexpr ex Expr(head := :return || :const, arg)
        return literal(head) * indent(1) * tile_expr(arg, 0)
    elseif @isexpr ex Expr(head := :local || :global, args...)
        return list_layout(Layout[tile_expr(arg, 0) for arg in args], prefix = head, par = (" ", ""))
    elseif @isexpr ex Expr(:do, call, Expr(:->, Expr(:tuple, args...), body))
        return tile_expr_do(call, args, body)
    elseif @isexpr ex Expr(:(=) || :kw, larg, rarg)
        return tile_expr_assign(:(=), larg, rarg, pr)
    elseif @isexpr ex Expr(head := :(:=) || :-> || :(<:) || :(::) || :(+=) || :(-=), larg, rarg)
        return tile_expr_assign(head, larg, rarg, pr)
    elseif @isexpr ex Expr(:call, fn := :(=>), larg, rarg)
        return tile_expr_assign(fn, larg, rarg, pr)
    elseif @isexpr ex Expr(:call, fn, args...)
        return tile_expr_call(fn, args, pr)
    elseif @isexpr ex Expr(fn := :|| || :&&, args...)
        return tile_expr_call(fn, args, pr)
    elseif @isexpr ex Expr(fn := :(::), arg)
        return tile_expr_call(fn, Any[arg], pr)
    elseif @isexpr ex Expr(fn := :(.), arg1, QuoteNode(arg2))
        return tile_expr_call(fn, [arg1, arg2], pr)
    elseif (@isexpr ex Expr(:(.), op::Symbol, Expr(:tuple, arg1, arg2))) && Meta.isoperator(op)
        return tile_expr_call(Symbol("." * string(op)), [arg1, arg2], pr)
    elseif @isexpr ex Expr(fn := :(.), arg1, arg2)
        return tile_expr_call(fn, [arg1, arg2], pr)
    elseif (@isexpr ex Expr(:comparison, args...)) && length(args) >= 3 && isodd(length(args))
        return tile_expr_comparison(args, pr)
    elseif @isexpr ex Expr(:..., arg)
        return pr > 0 ? literal("(") * tile_expr(arg, 0) * literal("...)") : tile_expr(arg, 0) * literal("...")
    elseif @isexpr ex Expr(:tuple, arg)
        return literal("(") * tile_expr(arg, 0) * literal(",)")
    elseif @isexpr ex Expr(:macrocall, name, args...)
        return tile_expr_macro(name, args, pr)
    elseif @isexpr ex Expr(:tuple, args...)
        return list_layout(Layout[tile_expr(arg, 0) for arg in args])
    elseif @isexpr ex Expr(:generator, call, loop)
        return list_layout(Layout[tile_expr(call, 0), tile_expr(loop, 0)],
                           sep = " for ", sep_brk = :start, par = ("", ""))
    elseif @isexpr ex Expr(:filter, pred, loop)
        return list_layout(Layout[tile_expr(loop, 0), tile_expr(pred, 0)],
                           sep = " if ", sep_brk = :start, par = ("", ""))
    elseif @isexpr ex Expr(:vect || :comprehension, args...)
        return list_layout(Layout[tile_expr(arg, 0) for arg in args], par = ("[", "]"))
    elseif @isexpr ex Expr(:ref || :typed_comprehension, t, args...)
        return tile_expr_ref(t, args)
    elseif @isexpr ex Expr(:curly, t, args...)
        return tile_expr_ref(t, args, par = ("{", "}"))
    elseif @isexpr ex Expr(:string, args...)
        return tile_expr_string(args)
    elseif @isexpr ex Expr(:quote, arg)
        return tile_expr_quoted(arg, pr)
    elseif @isexpr ex Expr(:$, arg)
        return literal("\$(") * tile_expr(arg) * literal(")")
    end
    tile_expr_fallback(ex)
end

tile_exprs(data) =
    tile_expr(data)

function tile_exprs(ex::Expr, parskip::Union{Bool,Nothing}=nothing)
    ex = strip_line_number_nodes(ex)
    if parskip === nothing
        parskip = @isexpr ex Expr(:toplevel, _...)
    end
    if @isexpr ex Expr(:toplevel || :block, args...)
        lts = Layout[]
        for arg in args
            if parskip && !isempty(lts)
                push!(lts, ZERO)
            end
            push!(lts, tile_expr(arg))
        end
        !isempty(lts) ? /(lts) : ZERO
    else
        tile_expr(ex)
    end
end

tile_expr_quoted(s::Symbol, pr) =
    literal(repr(s))

function tile_expr_quoted(ex, pr)
    if pr < 0
        lt = literal("quote")
        ex_lt = tile_exprs(ex)
        if ex_lt !== ZERO
            lt = lt / (indent(4) * ex_lt)
        end
        lt = lt / literal("end")
    else
        lt = literal(":(") * tile_expr(ex, 0) * literal(")")
    end
    lt
end

function tile_expr_block(args, pr)
    if pr < 0
        lt = literal("begin")
        if !isempty(args)
            lts = Layout[tile_expr(arg) for arg in args]
            lt = lt / (indent(4) * (/(lts)))
        end
        lt = lt / literal("end")
    else
        if isempty(args)
            lt = literal("(;;)")
        elseif length(args) == 1
            lt = tile_expr(args[1], pr)
        else
            lts = Layout[tile_expr(arg, 0) for arg in args]
            lt = list_layout(lts, sep = "; ")
        end
    end
    lt
end

function tile_expr_module(head, name, body)
    lt = literal(head) * indent(1) * literal(name)
    body_lt = tile_exprs(body, true)
    if body_lt !== ZERO
        lt = lt / ZERO / body_lt
    end
    lt = lt / ZERO / literal("end")
    lt
end

tile_expr_export(args) =
    list_layout(prefix = "export ", par = ("", ""),
                Layout[tile_expr(arg, 0) for arg in args])

tile_expr_import_path(ex) =
    (@isexpr ex Expr(:., args...)) ?
        list_layout([tile_expr(arg) for arg in args], sep=".", par=("", "")) :
        tile_expr(ex)

tile_expr_import(head, arg) =
    literal(head) * indent(1) * tile_expr_import_path(arg)

function tile_expr_import(head, from, args)
    prefix = literal(head) * indent(1) * tile_expr_import_path(from)
    list_layout(prefix=prefix, par=(": ", ""),
                Layout[tile_expr_import_path(arg) for arg in args])
end

function tile_expr_struct(head, decl, body)
    lt = literal(head) * indent(1) * tile_expr(decl, 0)
    body_lt = tile_exprs(body)
    if body_lt !== ZERO
        lt = lt / (indent(4) * body_lt)
    end
    lt = lt / literal("end")
    lt
end

tile_expr_abstract(decl) =
    literal("abstract type ") * tile_expr(decl, 0) / literal("end")

function tile_expr_stmt(head, cond, body, pr)
    lt = literal(head) * indent(1) * tile_expr(cond, 0)
    body_lt = tile_exprs(body)
    if body_lt !== ZERO
        lt = lt / (indent(4) * body_lt)
    end
    lt = lt / literal("end")
end

const if_precedence = Base.operator_precedence(:?)

function tile_expr_stmt(head, cond, body, else_body, pr)
    if pr >= 0 && head === :if && !(@isexpr body Expr(:block, _...)) && !(@isexpr else_body Expr(:block || :elseif, _...))
        cond_lt = tile_expr(cond, if_precedence)
        body_lt = tile_expr(body, if_precedence)
        else_body_lt = tile_expr(else_body, if_precedence)
        vlt = (cond_lt * literal(" ?")) /
              (indent(4) * body_lt * literal(" :")) /
              (indent(4) * else_body_lt)
        hlt = nobreak(cond_lt * literal(" ? ") * body_lt * literal(" : ") * else_body_lt)
        lt = hlt !== nothing ? (vlt | hlt) : vlt
        if pr >= if_precedence
            lt = literal("(") * lt * literal(")")
        end
        return lt
    end
    lt = literal(head) * indent(1) * tile_expr(cond, 0)
    body_lt = tile_exprs(body)
    if body_lt !== ZERO
        lt = lt / (indent(4) * body_lt)
    end
    if @isexpr else_body Expr(:elseif, cond′, body′)
        lt = lt / tile_expr_stmt(:elseif, cond′, body′, -1)
    elseif @isexpr else_body Expr(:elseif, cond′, body′, else_body′)
        lt = lt / tile_expr_stmt(:elseif, cond′, body′, else_body′, -1)
    else
        lt = lt / literal("else")
        else_body_lt = tile_exprs(else_body)
        if else_body_lt !== ZERO
            lt = lt / (indent(4) * else_body_lt)
        end
        lt = lt / literal("end")
    end
    lt
end

function tile_expr_let(decl, body)
    lt = literal(:let)
    if @isexpr decl Expr(:block, arg1, rest...)
        decl_lt = tile_expr(arg1, 0)
        for arg in rest
            decl_lt = decl_lt * literal(",") / tile_expr(arg, 0)
        end
    else
        decl_lt = tile_expr(decl)
    end
    if decl_lt !== ZERO
        lt = lt * indent(1) * decl_lt
    end
    body_lt = tile_exprs(body)
    if body_lt !== ZERO
        lt = lt / (indent(4) * body_lt)
    end
    lt = lt / literal("end")
    lt
end

function tile_expr_try(body, name, catch_body, finally_body)
    lt = literal("try")
    body_lt = tile_exprs(body)
    if body_lt !== ZERO
        lt = lt / (indent(4) * body_lt)
    end
    if catch_body !== nothing
        catch_lt = literal("catch")
        if name !== nothing
            catch_lt = catch_lt * indent(1) * tile_expr(name, 0)
        end
        lt = lt / catch_lt
        catch_body_lt = tile_exprs(catch_body)
        if catch_body_lt !== ZERO
            lt = lt / (indent(4) * catch_body_lt)
        end
    end
    if finally_body !== nothing
        lt = lt / literal("finally")
        finally_body_lt = tile_exprs(finally_body)
        if finally_body_lt !== ZERO
            lt = lt / (indent(4) * finally_body_lt)
        end
    end
    lt = lt / literal("end")
    lt
end

function tile_expr_function(head, decl, body)
    lt = literal(head) * indent(1) * tile_expr(decl)
    body_lt = tile_exprs(body)
    if body_lt !== ZERO
        lt = lt / (indent(4) * body_lt)
    end
    lt = lt / literal("end")
    lt
end

function tile_expr_do(call, args, body)
    call_lt = tile_expr(call, 0)
    body_lt = tile_exprs(body)
    if isempty(args)
        args_lt = literal(" do")
    else
        args_lt = literal(" do ") *
                  list_layout(Layout[tile_expr(arg, 0) for arg in args],
                              par = ("", ""))
    end
    if body_lt === ZERO
        call_lt * args_lt /
        literal("end")
    else
        call_lt * args_lt /
        (indent(4) * body_lt) /
        literal("end")
    end
end

function tile_expr_assign(op, larg, rarg, pr)
    pr′ = Base.operator_precedence(op)
    sep = op !== :(::) ? " $op " : string(op)
    lt = pair_layout(tile_expr(larg, pr′), tile_expr(rarg, pr′), sep = sep)
    if pr >= pr′
        lt = literal("(") * lt * literal(")")
    end
    lt
end

function tile_expr_call(fn, args, pr)
    nargs = length(args)
    if fn isa Symbol
        pr′ = Base.operator_precedence(fn)
        prefix = ""
        par = ("(", ")")
        if pr′ > pr && pr′ > 0
            par = ("", "")
        end
        sep = ", "
        if fn === :(:) && 2 <= nargs <= 3
            sep = ":"
        elseif (Base.isunaryoperator(fn) || fn === :(::)) && nargs == 1
            prefix = string(fn)
        elseif (Base.isbinaryoperator(fn) || fn === :in || fn === :isa) && nargs == 2
            if fn === :. || fn === :(::)
                sep = string(fn)
            else
                sep = " $fn "
            end
            assoc = Base.operator_associativity(fn)
            if assoc === :right
                while ((@isexpr args[end] Expr(:call, fn′::Symbol, arg1, arg2)) && fn′ === fn) ||
                      ((@isexpr args[end] Expr(fn′ := :|| || :&&, arg1, arg2)) && fn′ === fn)
                    pop!(args)
                    push!(args, arg1)
                    push!(args, arg2)
                end
            elseif assoc === :left
                while ((@isexpr args[1] Expr(:call, fn′::Symbol, arg1, arg2)) && fn′ === fn) ||
                      ((@isexpr args[1] Expr(fn′ := :., arg1, QuoteNode(arg2))) && fn′ === fn) ||
                      ((@isexpr args[1] Expr(fn′ := :. || :(::), arg1, arg2)) && fn′ === fn)
                    popfirst!(args)
                    pushfirst!(args, arg2)
                    pushfirst!(args, arg1)
                end
            end
        else
            prefix = identifier(fn)
            par = ("(", ")")
        end
        if length(args) >= 1 && @isexpr args[1] Expr(:parameters, params...)
            args = args[2:end]
            args_lt = list_layout(Layout[tile_expr(arg, 0) for arg in args], par = ("", ""))
            params_lt = list_layout(Layout[tile_expr(param, 0) for param in params], par = ("", ""))
            list_layout(Layout[args_lt, params_lt], sep = "; ", prefix = prefix, par = ("(", ")"))
        else
            list_layout(Layout[tile_expr(arg, pr′) for arg in args],
                        prefix = prefix,
                        par = par,
                        sep = sep)
        end
    else
        list_layout(Layout[tile_expr(arg, 0) for arg in args], prefix = tile_expr(fn))
    end
end

const comparison_precedence = Base.operator_precedence(:(==))

function tile_expr_comparison(args, pr)
    arg = pop!(args)
    vlt = hlt = tile_expr(arg, comparison_precedence)
    while !isempty(args)
        op = pop!(args)
        arg = pop!(args)
        arg_lt = tile_expr(arg, comparison_precedence)
        vlt = (arg_lt * literal(" $op")) / vlt
        hlt = arg_lt * literal(" $op ") * hlt
    end
    hlt = nobreak(hlt)
    lt = hlt !== nothing ? vlt | hlt : vlt
    if pr >= comparison_precedence
        lt = literal("(") * lt * literal(")")
    end
    lt
end

function tile_expr_macro(fn, args, pr)
    prefix = tile_expr(fn, pr)
    if pr < 0 && length(args) <= 1
        par = (" ", "")
        sep = " "
    else
        par = ("(", ")")
        sep = ", "
    end
    list_layout(Layout[tile_expr(arg, 0) for arg in args], prefix = prefix, par = par, sep = sep)
end

const dot_precedence = Base.operator_precedence(:.)

function tile_expr_ref(t, args; par = ("[", "]"))
    prefix = tile_expr(t, dot_precedence - 1)
    list_layout(Layout[tile_expr(arg, 0) for arg in args],
                prefix = prefix,
                par = par)
end

function tile_expr_string(chunks)
    lts = Layout[literal("\"")]
    for chunk in chunks
        if chunk isa AbstractString
            push!(lts, literal(escape_string(chunk, "\"\$")))
        else
            push!(lts, literal("\$("))
            push!(lts, tile_expr(chunk))
            push!(lts, literal(")"))
        end
    end
    push!(lts, literal("\""))
    *(lts...)
end

function tile_expr_fallback(ex::Expr)
    items = Layout[tile(ex.head)]
    for arg in ex.args
        push!(items, tile(arg))
    end
    list_layout(prefix = "\$", Layout[list_layout(prefix = "Expr", items)])
end


text:
using CSV
using DataFrames
using DelimitedFiles
using Glob
using Plots; pgfplotsx()

# Create plots
function plot_gmres_orth()

	# Plot labels
	labels = Dict(
		# kernels
		"ser" => "SparseArrays", "par" =>  "MtSpMV.jl", "mkl" => "",
		# formats
		"csr" => "", "csc" => "",
		# orth
		"mgs" => ", MGS", "cgs" => ", CGS", "dgks" => ", DGKS",
		# precon
		"none" => "", "jacobi" => "", "ilu0" => "",
		"parjacobi" => "", "parilu0" => "",
	)

   	# Plot linestyles
   	linestyles = [:solid, :dash, :dot]

   	# Benchmarks
   	benchs = [
   		"bench_gmres_orth_par_ser",
   		"bench_gmres_orth_ser_par",
   		"bench_gmres_orth_par_par"
   	]

   	# Folders
   	folders = [
   		joinpath("bench_par_ser", readdir("bench_par_ser")[1]),
   		joinpath("bench_ser_par", readdir("bench_ser_par")[1]),
		joinpath("bench_par_par", readdir("bench_par_par")[1])
	]

   	# Get result for serial spmv and serial dense (orth=mgs)
   	file = joinpath(folders[2], "csv", "runtimes_ser_csr_mgs_none.csv")
   	df = DataFrame!(CSV.File(file, header=true))
   	y = convert(Vector, df[!, 5])
   	y_ser_ser = similar(y)
   	y_ser_ser .= y[1]

   	# Iterate through folders
   	for (ifolder, folder) in enumerate(folders)

   		# Create plots
   		plt1 = plot(ylabel="Runtime in seconds", xlabel="Number of threads",
   					legend=:topright, framestyle=:box)
   		plt2 = plot(ylabel="Parallel speedup", xlabel="Number of threads",
   					legend=:none, framestyle=:box)

		# Get all csv files with runtimes in sorted order
		files = glob(string("runtimes*"), joinpath(folder, "csv"))
		if length(files) == 3
			files = [files[3], files[1], files[2]]
		end

		# Add results from all files to plots
		for (ifile, file) in enumerate(files)

			# Get benchmark results
		    df = DataFrame!(CSV.File(file, header=true))
		    x = convert(Vector, df[!, 1])
		    y = convert(Vector, df[!, 5])

		    # Create label from filename
		    infos = file[findlast('/', file) + 1 : findlast('.', file) - 1]
		    infos = infos[findfirst('_', infos) + 1 : length(infos)]
		    infos = split(infos, "_")
		    label = ""
		    for info in infos
				label = string(label, labels[info])
			end

			# Add runtimes to first plot
			plot!(plt1, x, y, label=label, linestyles=linestyles[ifile])

			# Add speedups to second plot
			y_speedups = y_ser_ser ./ y
			plot!(plt2, x, y_speedups, linestyles=linestyles[ifile])

		end # for

		# Plot side by side and save
		width = 530
		plot(plt1, plt2, layout=(1,2), size=(width, width/2))
		savefig(string(benchs[ifolder], ".png"))
		savefig(string(benchs[ifolder], ".tex"))

	end # for

end # function

max_stars_repo_name: jkrch/JuliaParGMRES-Plots

max_issues_repo_name: jkrch/JuliaParGMRES-Plots

max_forks_repo_name: jkrch/JuliaParGMRES-Plots


text:
using PkgServer
using Test

@testset "PkgServer.jl" begin
    # Write your own tests here.
end


text:
##############################
## Timedata display methods ##
##############################

import Base.Multimedia.display
function display(tn::AbstractTimedata)
    ## display information about an array

    ## set display parameters
    maxDispCols = 5;

    ## get type and field information
    typ = typeof(tn)
    println("\ntype: $typ")
    print("dimensions: ")
    print(size(tn))
    print("\n")

    ## get first entries
    (nrow, ncol) = size(tn)

    showCols = minimum([maxDispCols ncol]);

    Peekidx = DataFrame(idx = tn.idx);
    Peek = [Peekidx tn.vals[:, 1:showCols]];
    display(Peek)
end

function display(tn::AbstractTimedata, nCols::Int)
    ## display information about an array

    ## set display parameters
    maxDispCols = nCols;

    ## get type and field information
    typ = typeof(tn)
    println("\ntype: $typ")
    print("dimensions: ")
    print(size(tn))
    print("\n")

    ## get first entries
    (nrow, ncol) = size(tn)

    showCols = minimum([maxDispCols ncol]);

    Peekidx = DataFrame(idx = tn.idx);
    Peek = [Peekidx tn.vals[:, 1:showCols]];
    display(Peek)
end

#######################
## display in IJulia ##
#######################

### Remark: timematr has its own display function!!
import Base.writemime
function writemime(io::IO,
                        ::MIME"text/html",
                        td::AbstractTimedata)

    typ = typeof(td)

    ## set display parameters
    maxDispCols = 5;
    (nrow, ncol) = size(td)

    write(io, "
$typ
")
    ## write(io, "
Dimensions: ($nrow, $ncol)
")
    write(io, "
Dimensions: ($nrow, $ncol)
")

    fromDate = td.idx[1];
    toDate = td.idx[end];

    ## write(io, "
-------------------------------------------
")
    ## write(io, "
From: $fromDate,
    ## To: $toDate
")
    write(io, "
From: $fromDate, To: $toDate
")
    ## write(io, "
-------------------------------------------
")

    showCols = minimum([maxDispCols (ncol+1)]);

    df = convert(DataFrame, td)

    writemime(io, "text/html", df[:, 1:showCols])
end

function Base.writemime(io::IO,
                        ::MIME"text/html",
                        tm::AbstractTimematr)

    # show only rounded values
    signDigits = 3

    typ = typeof(tm)

    ## set display parameters
    maxDispCols = 8;
    (nrow, ncol) = size(tm)

    write(io, "
$typ
")
    ## write(io, "
Dimensions: ($nrow, $ncol)
")
    write(io, "
Dimensions: ($nrow, $ncol)
")

    fromDate = tm.idx[1];
    toDate = tm.idx[end];

    ## write(io, "
-------------------------------------------
")
    ## write(io, "
From: $fromDate,
    ## To: $toDate
")
    write(io, "
From: $fromDate, To: $toDate
")
    ## write(io, "
-------------------------------------------
")

    showCols = minimum([maxDispCols (ncol+1)]);

    rndVals = round(core(tm), signDigits)
    tm2 = Timematr(rndVals, names(tm), idx(tm))

    df = convert(DataFrame, tm2)

    writemime(io, "text/html", df[:, 1:showCols])
end

##################
## html display ##
##################

type HTML
   s::String
end
import Base.writemime
writemime(io::IO, ::MIME"text/html", x::HTML) = print(io, x.s)

macro table(title::String, expr::Union(Expr, Symbol))
    quote
        titleString = $title
        htmlTitleString = string("
",
                                 titleString,
                                 "
")
        display(HTML(htmlTitleString))
        $(esc(expr))
    end
end

#########
## str ##
#########

function str(tn::AbstractTimedata)
    ## display information about an array

    ## set display parameters
    maxDispRows = Inf;
    maxDispCols = 5;

    ## get type and field information
    typ = typeof(tn)
    println("\ntype: $typ")
    nms = ["vals", "idx"]
    ## types = DataType[]
    for fieldname in nms
        tp = fieldtype(tn, symbol(fieldname))
        ## push!(types, fieldtype(tn, fieldname))
        println(":$fieldname  \t\t  $tp")
    end

    print("\ndimensions of vals: ")
    print(size(tn))

    ## get min and max
    ## minVal = minimum(tn.vals);
    ## maxVal = maximum(tn.vals);

    fromDate = tn.idx[1];
    toDate = tn.idx[end];

    println("\n\n-------------------------------------------")
    println("From: $fromDate, To: $toDate")
    println("-------------------------------------------\n")

    ## get first entries
    (nrow, ncol) = size(tn)

    showCols = minimum([maxDispCols ncol]);

    Peekidx = DataFrame(idx = tn.idx);
    Peek = [Peekidx tn.vals[:, 1:showCols]];
    display(Peek)
end