text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

# -------------
# COMPOSITIONS
# -------------

"""
    clr(c)

Centered log-ratio transformation of composition `c`.
"""
function clr(c::Composition{D}) where {D}
  w = components(c) .+ eps()
  μ = geomean(w)
  SVector(ntuple(i -> log(w[i] / μ), D))
end

"""
    clrinv(x)

Inverse clr transformation of coordinates `x`.
"""
clrinv(x::SVector{D,T}) where {D,T<:Real} =
  Composition(𝓒(exp.(x)))

# -------
# TABLES
# -------

"""
    CLR()

Centered log-ratio transform following the
[TableTransforms.jl](https://github.com/JuliaML/TableTransforms.jl)
interface.
"""
struct CLR <: LogRatio end

refvar(::CLR, vars) = last(vars)

newvars(::CLR, n) = collect(n)

oldvars(::CLR, vars, rvar) = collect(vars)

function applymatrix(::CLR, X)
  μ = geomean.(eachrow(X))
  L = log.(X .+ eps())
  l = log.(μ .+ eps())
  L .- l
end

function revertmatrix(::CLR, Y)
  E = exp.(Y)
  mapslices(𝓒, E, dims=2)
end


text:
function part1(input)
    rules = parse_rules(input)
    return count(x -> "shiny gold" in x,
                 values(transitive_contained_colors(rules)))
end

function part2(input)
    rules = parse_rules(input)
    required_bags = 0
    queue = Dict("shiny gold" => 1)
    while !isempty(queue)
        color, quantity = pop!(queue)
        required_bags += quantity
        for (color2, quantity2) in rules[color]
            queue[color2] = get(queue, color2, 0) + quantity * quantity2
        end
    end
    return required_bags - 1
end

function parse_rules(input)
    rules = Dict{String, Dict{String, Int}}()
    for line in readlines(input)
        parse_rule!(rules, line)
    end
    return rules
end

function parse_rule!(rules, line)
    outer, inner = match(r"(.+) bags contain (.+)", line).captures
    rule = get!(rules, outer, Dict{String, Int}())
    inner == "no other bags." && return
    for contained in split(inner, ", ")
        quantity, color = match(r"(\d)+ (.+) bag", contained).captures
        rule[color] = parse(Int, quantity)
    end
end

function transitive_contained_colors(rules)
    contained_colors = Dict{String, Vector{String}}()
    for color in keys(rules)
        found_colors = collect(keys(rules[color]))
        i = 1
        while i <= length(found_colors)
            new_colors = setdiff(keys(rules[found_colors[i]]), found_colors)
            append!(found_colors, new_colors)
            i += 1
        end
        contained_colors[color] = found_colors
    end
    return contained_colors
end

ext: jl


text:
"""
    is_choi_matrix(JN :: AbstractMatrix, dimA :: Int, dimB :: Int) :: Bool

Returns `true` if the supplied matrix `JN` is a Choi operator.
This function returns `false` if
* `size(JN) != (dimA * dimB, dimA * dimB)`
"""
function is_choi_matrix(JN :: AbstractMatrix, dimA :: Int, dimB :: Int) :: Bool
    JN_dim = dimA*dimB
    if size(JN) != (JN_dim, JN_dim)
        return false
    end

    return true
end

"""
    Choi( JN :: AbstractMatrix, in_dim :: Int, out_dim :: Int ) :: Choi{<:Number}

    Choi( N :: Function, in_dim :: Int, out_dim :: Int ) :: Choi{ComplexF64}

    Choi( kraus_ops :: Vector{<:AbstractMatrix} ) :: Choi{ComplexF64}

Constructs the Choi matrix representation of a quantum channel.
If either a function `N` or set of kraus operators is provided as input, the
Choi matrix is constructed with the [`choi`](@ref) method.

The `Choi` type contains the fields:
* `JN :: Matrix{<:Number}` - The choi matrix.
* `in_dim :: Int` - The channel's input dimension.
* `out_dim :: Int` - The Channel's output dimension.

A `DomainError` is thrown if [`is_choi_matrix`](@ref) returns `false`.
"""
struct Choi{T<:Number}
    JN :: Matrix{T}
    in_dim :: Int
    out_dim :: Int
    Choi(
        JN :: AbstractMatrix,
        in_dim :: Int,
        out_dim :: Int
    ) = is_choi_matrix(JN, in_dim, out_dim) ? new{eltype(JN)}(JN, in_dim, out_dim) : throw(
        DomainError(JN, "The Choi matrix dimensions are not valid.")
    )
    Choi(
        N :: Function,
        in_dim :: Int,
        out_dim :: Int
    ) = Choi( choi(N, in_dim, out_dim), in_dim, out_dim)
    Choi(
        kraus_ops :: Vector{<:AbstractMatrix},
    ) = Choi(choi(kraus_ops), size(kraus_ops[1])...)
end

# print out matrix forms when Choi types are displayed
function show(io::IO, mime::MIME{Symbol("text/plain")}, choi_op :: Choi)
    summary(io, choi_op)
    print(io, "\nin_dim : ", choi_op.in_dim, ", out_dim : ", choi_op.out_dim)
    print(io, "\nJN : ")
    show(io, mime, choi_op.JN)
end

"""
    parChoi(chan1 :: Choi, chan2 :: Choi) :: Choi

Returns the tensor product of two [`Choi`](@ref) matrices

```math
    J^{AB}_{\\mathcal{N}}\\otimes J^{A'B'}_{\\mathcal{M}} \\to
    J^{AA'BB'}_{\\mathcal{N}\\otimes\\mathcal{M}}
```

where ``J^{A:B}_{\\mathcal{N}}`` and ``J^{A':B'}_{\\mathcal{M}}`` are the Choi
matrices for `chan1` and `chan2` respectively.
Note the implicit swap between systems ``B \\leftrightarrow A'``.
"""
function parChoi(chan1 :: Choi, chan2 :: Choi) :: Choi
    par_dims = [chan1.in_dim, chan1.out_dim, chan2.in_dim, chan2.out_dim]
    par_JN = permuteSubsystems(kron(chan1.JN, chan2.JN), [1,3,2,4], par_dims)
    par_in_dim = chan1.in_dim * chan2.in_dim
    par_out_dim = chan1.out_dim * chan2.out_dim

    return Choi(par_JN, par_in_dim, par_out_dim)
end

"""
    choi(𝒩 :: Function, Σ :: Int, Λ :: Int) :: Matrix{ComplexF64}
    choi(kraus_ops :: Vector{<:AbstractMatrix}) :: Matrix{ComplexF64}

This function returns the Choi state of a channel represented either as a generic
function `𝒩` or a set of Kraus operators `kraus_ops`. It does this using that

```math
        J_{\\mathcal{N}} = \\sum_{a,b \\in \\Sigma} E_{a,b} \\otimes \\mathcal{N}(E_{a,b}) ,
```

where ``\\Sigma`` is the finite alphabet indexing the input space and ``E_{a,b}``
is a square matrix of dimension ``\\Sigma`` with a ``1`` in the ``(a,b)`` entry
and a ``0`` everywhere else. The input ``\\Lambda`` is the output dimension.
Note this assumes you have a function that calculates
``\\mathcal{N}(X)`` for arbitrary input ``X``. As many of the functions for channels
in this module have multiple parameters, please note that if you have a channel function
`𝒩(ρ, p, q)` that calculates ``\\mathcal{N}_{p,q}(\\rho)``, you can declare a function
`𝒩_xy(ρ) = 𝒩(ρ,x,y)` for fixed `(x,y)` and then call, `choi(𝒩_xy, Σ)`.
"""
function choi(𝒩 :: Function, Σ :: Int, Λ :: Int) :: Matrix{ComplexF64}
    eab_matrix = zeros(ComplexF64, Σ,Σ)
    choi_matrix = zeros(ComplexF64, Σ*Λ,Σ*Λ)
    for i in 1:Σ
        row_ids = ((i-1) * Λ + 1):((i-1) * Λ + Λ)
        for j in 1:Σ
            col_ids = ((j-1) * Λ + 1):((j-1) * Λ + Λ)

            eab_matrix[i,j] = 1
            choi_matrix[row_ids,col_ids] += 𝒩(eab_matrix)
            eab_matrix[i,j] = 0
        end
    end
    return choi_matrix
end
function choi(kraus_ops :: Vector{<:AbstractMatrix}) :: Matrix{ComplexF64}
    kraus_channel(ρ) = krausAction(kraus_ops, ρ)
    dim_out, dim_in = size(kraus_ops[1])

    choi(kraus_channel, dim_in, dim_out)
end

"""
    isometricChannel(kraus_ops :: Vector) :: Matrix

This function builds the isometric representation ``V`` of a
channel ``\\mathcal{N}: A \\to B`` from the Kraus operators ``\\{K_{i}\\}``.
It does this by calculating
```math
    V = \\sum_{i} K_{i} \\otimes |i\\rangle
```
"""
function isometricChannel(kraus_ops :: Vector) :: Matrix
    dimB , dimA = size(kraus_ops[1])
    dimE = length(kraus_ops)
    V, e_vec = zeros(dimE*dimB,dimA), zeros(dimE)
    for i in 0:dimE-1;
        V[1+i*dimB:(i+1)*dimB,1:dimA] = kraus_ops[i+1]
    end
    label_vec = [1:dimB*dimE;]
    perms_ids_vec = permuteSubsystems(label_vec,[2,1],[dimE,dimB])
    V = V[perms_ids_vec,:]

    return V
end

"""
    complementaryChannel(kraus_ops :: Vector) :: Vector

This function takes a set of Kraus operators for a channel ``\\mathcal{N}_{A \\to B}``
and returns a set of Kraus operators for the complementary channel,
``\\mathcal{N}^{c}_{A \\to E}``. It does this by generating the Kraus operators of
the isometric representation of the channel followed by partial trace on the
``B`` space.
!!! info
    If ``\\mathcal{N}_{A \\to B}`` is already isometric, the code lets dimE=2
    so that functions are well behaved.
"""
function complementaryChannel(kraus_ops :: Vector) :: Vector
    V = isometricChannel(kraus_ops)
    dimB , dimA = size(kraus_ops[1])
    dimE = length(kraus_ops)

    comp_kraus = Vector{Matrix{Complex{Float64}}}(undef, dimB)
    if dimE == 1
        #If dimE = 1, the complementary channel is the replacer channel
        #To get it to work with other functions, we let dimE=2 so it isn't scalar
        b_vec = zeros(dimB)
        for i in [1:dimB;]
            b_vec[i] = 1
            comp_kraus[i] = kron([1;0],b_vec')
            b_vec[i] = 0
        end
    else
        b_vec, e_id = zeros(dimB), Matrix(1I,dimE,dimE)
        for i in [1:dimB;]
            b_vec[i] = 1
            comp_kraus[i] = kron(b_vec,e_id)'*V
            b_vec[i] = 0
        end
    end

    return comp_kraus
end

"""
    krausAction(kraus_ops :: Vector, X)

This function takes a set of Kraus operators for a channel
``\\mathcal{N}_{A \\to B}`` and returns the output of the channel
for input ``X``. That is, given input ``X`` it returns
```math
    \\mathcal{N}_{A \\to B}(X) = \\sum_{i} K_{i} X K_{i}^{\\ast}
```
"""
function krausAction(kraus_ops :: Vector, X)
    dimB, dimA = size(kraus_ops[1])
    out = zeros(dimB,dimB)
    for i = [1:length(kraus_ops);]
        out += kraus_ops[i]*X*kraus_ops[i]'
    end
    return out
end

"""
    depolarizingChannel(ρ :: Matrix{Float64}, q :: Union{Int,Float64}) :: Matrix{ComplexF64}

This calculates the action of the depolarizing channel,

```math
\\Delta_{q}(\\rho) = (1-q)\\rho + q \\text{Tr}(\\rho) \\frac{1}{d} I_{AB} ,
```

where ``q \\in [0,1].``
Note these channels are the channels covariant with respect to the unitary group.

A `DomainError` is thrown if:
* Matrix `ρ` is not square
* Input `q` does not satisfy `0 ≤ q ≤ 1`
"""
function depolarizingChannel(ρ :: Matrix{<:Number}, q :: Union{Int,Float64}) :: Matrix{ComplexF64}
    dim = size(ρ,1)
    if !isequal(dim,size(ρ,2))
        throw(DomainError(ρ, "the input ρ is not a square matrix"))
    elseif !(0 ≤ q ≤ 1)
        throw(DomainError(q, "depolarizingChannel requires q ∈ [0,1]."))
    end

    return (1-q)*ρ  + q*(tr(ρ))*(1/(dim))*I
end

"""
    dephrasureChannel(
        ρ :: Matrix{<:Number},
        p :: Union{Int,Float64},
        q :: Union{Int,Float64}
    ) :: Matrix{ComplexF64}

This function calculates the action of the [dephrasureChannel](https://arxiv.org/abs/1806.08327),

```math
\\mathcal{N}_{p,q}( \\rho) := (1-q)((1-p) \\rho + pZ \\rho Z) + q \\text{Tr}( \\rho) |e\\rangle \\langle e|,
```

where ``p,q \\in [0,1]``, ``Z`` is the Pauli-Z matrix, and ``|e\\rangle\\langle e|``
an  error  flag orthogonal to the  Hilbert space of input state ``\\rho``.

A `DomainError` is thrown if:
* Matrix `ρ` is not `2x2`
* Inputs `p` or `q` do not satisdy `0 ≤ p,q ≤ 1`
"""
function dephrasureChannel(ρ :: Matrix{<:Number},p :: Union{Int,Float64}, q :: Union{Int,Float64}) :: Matrix{ComplexF64}
    if ((size(ρ,1)!=2)||(size(ρ,2)!=2))
        throw(DomainError(ρ, "the input ρ is not a qubit"))
    elseif !(0 ≤ q ≤ 1)
        throw(DomainError(q, "dephrasureChannel requires q ∈ [0,1]."))
    elseif !(0 ≤ p ≤ 1)
        throw(DomainError(p, "dephrasureChannel requires p ∈ [0,1]."))
    end
    pauli_Z = [1 0 ; 0 -1]
    output_ρ = zeros(ComplexF64,3,3)
    output_ρ[1:2,1:2]=(1-q)*((1-p)*ρ + p*pauli_Z*ρ*pauli_Z)
    output_ρ[3,3] = q*tr(ρ)
    return output_ρ
end

"""
    wernerHolevoChannel(ρ :: Matrix{<:Number}, λ :: Union{Int,Float64}) :: Matrix{ComplexF64}

This function calculates the action of the [generalized Werner-Holevo channels](https://arxiv.org/abs/1406.7142)

```math
    \\mathcal{W}_{d,\\lambda}(\\rho) = \\lambda \\Phi_{0}(ρ) + (1-\\lambda) \\Phi_{1}(ρ)
```

where ``\\lambda \\in [0,1]``. This means these are convex combinations of the original [Werner-Holevo channels](https://arxiv.org/abs/quant-ph/0203003)
which are defined as

```math
    \\Phi_{0}(ρ) = \\frac{1}{d+1}(\\text{Tr}[\\rho]I_{d} +ρ^{T}) \\hspace{1cm}
    \\Phi_{1}(ρ) = \\frac{1}{d-1}(\\text{Tr}[\\rho]I_{d} -ρ^{T}) .
```

Note the Choi matrices of these generalized channels are the (unnormalized) Werner states.

A `DomainError` is thrown if:
* Matrix `ρ` is not square
* `p` is not in  range `0 ≤ p ≤ 1`
"""
function wernerHolevoChannel(ρ :: Matrix{<:Number}, λ :: Union{Int,Float64}) :: Matrix{ComplexF64}
    if !isequal(size(ρ)...)
        throw(DomainError(ρ, "the input ρ is not a square matrix"))
    elseif !(0 ≤ λ ≤ 1)
        throw(DomainError(λ, "wernerHolevoChannel requires λ ∈ [0,1]."))
    end
    dim = size(ρ,1)
    term_1 = 1 / (dim+1) * (tr(ρ) * I + transpose(ρ))
    term_2 = 1 / (dim-1) * (tr(ρ) * I - transpose(ρ))
    return λ * term_1 + (1-λ) * term_2
end
#The aligned below makes the documenter space the bmatrix properly
"""
    siddhuChannel(ρ :: Matrix{<:Number}, s :: Union{Int,Float64}) :: Matrix{<:Number}

This function calculates the action of the Siddhu channel ``N_{s}`` which is defined by Kraus operators:
```math
    \\begin{aligned}
    K_{0} = \\begin{bmatrix} \\sqrt{s} & 0 & 0 \\\\ 0 & 0 & 0 \\\\ 0 & 1 & 0 \\end{bmatrix}
    \\hspace{5mm}
    K_{1} = \\begin{bmatrix} 0 & 0 & 0 \\\\ \\sqrt{1-s} & 0 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix} ,
    \\end{aligned}
```
where ``s \\in [0,1/2]``.
This channel was introduced in Equation 9 of [this paper](https://arxiv.org/abs/2003.10367).
"""
function siddhuChannel(ρ :: Matrix{<:Number}, s :: Union{Int,Float64}) :: Matrix{<:Number}
    if !isequal(size(ρ)...)
        throw(DomainError(ρ, "the input ρ is not a square matrix"))
    elseif size(ρ)[1] != 3
        throw(DomainError(ρ, "The input must be a qutrit operator."))
    elseif !(0 ≤ s ≤ 1/2)
        throw(DomainError(s, "siddhuChannel requires s ∈ [0,1/2]."))
    end

    K0 = [sqrt(s) 0 0 ; 0 0 0 ; 0 1 0]
    K1 = [0 0 0 ; sqrt(1-s) 0 0 ; 0 0 1]
    return K0*ρ*K0' + K1*ρ*K1'
end

"""
    generalizedSiddhu(
          ρ :: Matrix{<:Number},
          s :: Union{Int,Float64},
          μ :: Union{Int,Float64}
    ) :: Matrix{<:Number}

This function calculates the action of the generalized Siddhu channel
``\\mathcal{N}_{s,\\mu}`` on the qutrit state ``\\rho``.
The action of the channel is defined by Kraus operators:
```math
    \\begin{aligned}
    K_{0} = \\begin{bmatrix}
                \\sqrt{s} & 0 & 0 \\\\
                0 & \\sqrt{1-\\mu} & 0 \\\\
                0 & 0 & \\sqrt{\\mu}
             \\end{bmatrix}
    \\hspace{5mm}
    K_{1} = \\begin{bmatrix}
                0 & 0 & \\sqrt{1-\\mu} \\\\
                \\sqrt{1-s} & 0 & 0 \\\\
                0 & \\sqrt{\\mu} & 0
            \\end{bmatrix} ,
    \\end{aligned}
```
where ``s \\in [0,1/2]``, ``\\mu \\in [0,1]``.
This channel was introduced by Leditzky *et al.* (cite when on arxiv).

A `DomainError` is thrown if:
* Matrix `ρ` is not square
* `ρ` is not a 3-dimensional matrix
* `s` is not in range `0 ≤ s ≤ 1/2`
* `μ` is not in range `0 ≤ μ ≤ 1`
"""
function generalizedSiddhu(ρ :: Matrix{<:Number}, s :: Union{Int,Float64}, μ :: Union{Int,Float64}) :: Matrix{<:Number}
    if !isequal(size(ρ)...)
        throw(DomainError(ρ, "the input ρ is not a square matrix"))
    elseif size(ρ)[1] != 3
        throw(DomainError(ρ, "The input must be a qutrit operator."))
    elseif !(0 ≤ s ≤ 1/2)
        throw(DomainError(s, "generalizedSiddhu requires s ∈ [0,1/2]."))
    elseif !(0 ≤ μ ≤ 1)
        throw(DomainError(μ, "generalizedSiddhu requires μ ∈ [0,1]."))
    end

    K0 = [sqrt(s) 0 0 ; 0 sqrt(1-μ) 0 ; 0 0 sqrt(μ)]
    K1 = [0 0 sqrt(1-μ) ; sqrt(1-s) 0 0 ; 0 sqrt(μ) 0]
    return K0*ρ*K0' + K1*ρ*K1'
end

"""
    GADChannel(
        ρ :: Matrix{<:Number},
        p :: Union{Int,Float64},
        n :: Union{Int,Float64}
    ) :: Matrix{<:Number}

This function calculates the action of the generalized (qubit) amplitude damping channel ``\\mathcal{A}_{p,n}``
which is defined by Kraus operators:
```math
    \\begin{aligned}
    K_{0} =& \\sqrt{1-n} \\begin{bmatrix} 1 & 0 \\\\ 0 & \\sqrt{1-p} \\end{bmatrix}
    \\hspace{5mm}
    K_{1} =& \\sqrt{p(1-n)} \\begin{bmatrix} 0 & 1 \\\\ 0 & 0 \\end{bmatrix}  \\\\
    K_{2} =& \\sqrt{n} \\begin{bmatrix} \\sqrt{1-p} & 0 \\\\ 0 & 1 \\end{bmatrix}
    \\hspace{5mm}
    K_{3} =& \\sqrt{pn} \\begin{bmatrix} 0 & 0 \\\\ 1 & 0 \\end{bmatrix}
    \\end{aligned}
```
where ``p,n \\in [0,1]``.
This channel may be found in Section 3 of [this paper](https://arxiv.org/abs/2107.13486).
"""
function GADChannel(ρ :: Matrix{<:Number}, p :: Union{Int,Float64}, n :: Union{Int,Float64}) :: Matrix{<:Number}
    if !isequal(size(ρ)...)
        throw(DomainError(ρ, "the input ρ is not a square matrix"))
    elseif size(ρ)[1] != 2
        throw(DomainError(ρ, "The input must be a qubit operator."))
    elseif !(0 ≤ p ≤ 1)
        throw(DomainError(p, "siddhuChannel requires p ∈ [0,1]."))
    elseif !(0 ≤ n ≤ 1)
        throw(DomainError(n, "siddhuChannel requires n ∈ [0,1]."))
    end

    K0 = sqrt(1-n)*[1 0 ; 0 sqrt(1-p)]
    K1 = sqrt(p*(1-n))*[0 1 ; 0 0]
    K2 = sqrt(n)*[sqrt(1-p) 0 ; 0 1]
    K3 = sqrt(p*n)*[0 0 ; 1 0]
    return K0*ρ*K0' + K1*ρ*K1' + K2*ρ*K2' + K3*ρ*K3'
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
export CalculateSampling,
Pad5D,
IRLS,
FFTOp,
WeightingOp,
ConjugateGradients,
LinearOperator,
MatrixMultiplyOp,
InnerProduct,
DotTest,
FISTA,
power_method
include("CalculateSampling.jl")
include("Pad5D.jl")
include("IRLS.jl")
include("FFTOp.jl")
include("WeightingOp.jl")
include("ConjugateGradients.jl")
include("LinearOperator.jl")
include("MatrixMultiplyOp.jl")
include("InnerProduct.jl")
include("DotTest.jl")
include("FISTA.jl")
include("power_method.jl")


text:
# Helpers
id_or_chain(obj::AbstractString, kind::AbstractString) = kind * "id=" * obj * "&"
id_or_chain(obj::AbstractVector{<:AbstractString}, kind::AbstractString) = reduce(*, kind .* "id=" .* obj .* "&")
function period(startdate::TimeType, enddate::TimeType, dataset::AbstractString)
    starts = string.(startdate:ifelse(dataset ∈ ["GSOM", "GSOY"], Year(10), Year(1)):enddate)
    "&startdate=" .* starts .* "&enddate=" .* push!(starts[2:end], string(enddate))
end
for_assignment(::Type{Union{Missing, T}}, value::Missing) where T <: Union{AbstractString, Number, TimeType} = missing
for_assignment(::Type{Union{Missing, T}}, value::AbstractString) where T <: AbstractString = convert(T, value)
for_assignment(::Type{Union{Missing, T}}, value::Number) where T <: Number = convert(T, value)
for_assignment(::Type{Union{Missing, T}}, value::AbstractString) where T <: TimeType = T(convert(String, value))
for_assignment(::Type{T}, value::Missing) where T <: Union{AbstractString, Number, TimeType} = missing
for_assignment(::Type{T}, value::AbstractString) where T <: AbstractString = convert(T, value)
for_assignment(::Type{T}, value::Number) where T <: Number = convert(T, value)
for_assignment(::Type{T}, value::AbstractString) where T <: TimeType = T(convert(String, value))

isvalid_cdotoken(obj::AbstractString) = occursin(r"[A-Za-z]{32}", obj)
isvalid_datsets(obj::AbstractVector{<:AbstractString}) = all(isvalid_datsets, obj)
isvalid_datsets(obj::AbstractString) = obj ∈ DATASETS
isvalid_datacategories(obj::AbstractVector{<:AbstractString}) = all(isvalid_datacategories, obj)
isvalid_datacategories(obj::AbstractString) = obj ∈ DATACATEGORIES
isvalid_datatypes(obj::AbstractVector{<:AbstractString}) = all(elem -> occursin(r"^[A-Z0-9-]{3,}$", elem), obj)
isvalid_datatypes(obj::AbstractString) = occursin(r"^[A-Z0-9-]{3,}$", obj)
isvalid_locationcategories(obj::AbstractVector{<:AbstractString}) = all(isvalid_locationcategories, obj)
isvalid_locationcategories(obj::AbstractString) = obj ∈ LOCATIONCATEGORIES
isvalid_locations(obj::AbstractVector{<:AbstractString}) = all(isvalid_locations, obj)
isvalid_locations(obj::AbstractString) = occursin(r"^(CITY|CLIM_DIV|CLIM_REG|CNTRY|CNTY|FIPS|HYD_ACC|HYD_CAT|HYD_REG|HYD_SUB|ST|US_TERR|ZIP):[A-Z0-9]{1,}$", obj)
isvalid_stations(obj::AbstractVector{<:AbstractString}) = all(isvalid_stations, obj)
isvalid_stations(obj::AbstractString) = occursin(r"^(COOP|GHCND|NEXRAD|WBAN):[a-zA-Z0-9-_]{1,}$", obj)
isvalid_extent(obj::AbstractVector{<:AbstractFloat}) = length(obj) ∈ [0, 4]

# Struct
abstract type Endpoint end
struct CDO_Data <: Endpoint
    token::String
    url::Vector{String}
    CDO_Data(CDO_token::AbstractString,
             dataset::AbstractString,
             startdate::TimeType,
             enddate::TimeType,
             datatypes::Union{AbstractString, AbstractVector{<:AbstractString}},
             locations::Union{AbstractString, AbstractVector{<:AbstractString}},
             stations::Union{AbstractString, AbstractVector{<:AbstractString}},
             metric::Bool) =
        new(CDO_token,
            string.("https://www.ncdc.noaa.gov/cdo-web/api/v2/data?datasetid=",
                    dataset,
                    "&",
                    id_or_chain(datatypes, "datatype"),
                    id_or_chain(locations, "location"),
                    id_or_chain(stations, "station"),
                    "units=",
                    ifelse(metric, "metric", "standard"),
                    period(startdate, enddate, dataset),
                    "&limit=1000&offset=1"))
end
struct CDO_DataCategory <: Endpoint
    token::String
    url::String
    CDO_DataCategory(CDO_token::AbstractString, datacategory::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/datacategories/" * datacategory)
end
struct CDO_DataCategories <: Endpoint
    token::String
    url::String
    CDO_DataCategories(CDO_token::AbstractString,
                       datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                       locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                       stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                       startdate::Date,
                       enddate::Date) =
        new(CDO_token,
            string("https://www.ncdc.noaa.gov/cdo-web/api/v2/datacategories?",
                   id_or_chain(datasets, "dataset"),
                   id_or_chain(locations, "location"),
                   id_or_chain(stations, "station"),
                   "startdate=",
                   startdate,
                   "&enddate=",
                   enddate,
                   "&limit=1000&offset=1"))
end
struct CDO_Dataset <: Endpoint
    token::String
    url::String
    CDO_Dataset(CDO_token::AbstractString, dataset::AbstractString) = new(CDO_token, "https://www.ncdc.noaa.gov/cdo-web/api/v2/datasets/" * dataset)
end
struct CDO_Datasets <: Endpoint
    token::String
    url::String
    CDO_Datasets(CDO_token::AbstractString,
                 datatypes::Union{AbstractString, AbstractVector{<:AbstractString}},
                 locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                 stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                 startdate::Date,
                 enddate::Date) =
        new(CDO_token,
            string("https://www.ncdc.noaa.gov/cdo-web/api/v2/datasets?",
                   id_or_chain(datatypes, "datatype"),
                   id_or_chain(locations, "location"),
                   id_or_chain(stations, "station"),
                   "startdate=",
                   startdate,
                   "&enddate=",
                   enddate,
                   "&limit=1000&offset=1"))
end
struct CDO_DataType <: Endpoint
    token::String
    url::String
    CDO_DataType(CDO_token::AbstractString, datatype::AbstractString) =
        new(CDO_token, string("https://www.ncdc.noaa.gov/cdo-web/api/v2/datatypes/", datatype))
end
struct CDO_DataTypes <: Endpoint
    token::String
    url::String
    CDO_DataTypes(CDO_token::AbstractString,
                  datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                  datacategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                  locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                  stations::Union{AbstractString, AbstractVector{<:AbstractString}},
                  startdate::Date,
                  enddate::Date) =
        new(CDO_token,
            string("https://www.ncdc.noaa.gov/cdo-web/api/v2/datatypes?",
                   id_or_chain(datasets, "dataset"),
                   id_or_chain(datacategories, "datacategory"),
                   id_or_chain(locations, "location"),
                   id_or_chain(stations, "station"),
                   "startdate=",
                   startdate,
                   "&enddate=",
                   enddate,
                   "&limit=1000&offset=1"))
end
struct CDO_LocationCategory <: Endpoint
    token::String
    url::String
    CDO_LocationCategory(token::AbstractString, locationcategory::AbstractString) =
        new(token, string("https://www.ncdc.noaa.gov/cdo-web/api/v2/locationcategories/", locationcategory))
end
struct CDO_LocationCategories <: Endpoint
    token::String
    url::String
    CDO_LocationCategories(CDO_token::AbstractString,
                           datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                           startdate::Date,
                           enddate::Date) =
        new(CDO_token,
            string("https://www.ncdc.noaa.gov/cdo-web/api/v2/locationcategories?",
                   id_or_chain(datasets, "dataset"),
                   "startdate=",
                   startdate,
                   "&enddate=",
                   enddate,
                   "&limit=1000&offset=1"))
end
struct CDO_Location <: Endpoint
    token::String
    url::String
    CDO_Location(CDO_token::AbstractString, location::AbstractString) =
        new(CDO_token, string("https://www.ncdc.noaa.gov/cdo-web/api/v2/locations/", location))
end
struct CDO_Locations <: Endpoint
    token::String
    url::String
    CDO_Locations(CDO_token::AbstractString,
                  datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                  locationcategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                  datacategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                  startdate::Date,
                  enddate::Date) =
        new(CDO_token,
            string("https://www.ncdc.noaa.gov/cdo-web/api/v2/locations?",
                   id_or_chain(datasets, "dataset"),
                   id_or_chain(locationcategories, "locationcategory"),
                   id_or_chain(datacategories, "datacategory"),
                   "startdate=",
                   startdate,
                   "&enddate=",
                   enddate,
                   "&limit=1000&offset=1"))
end
struct CDO_Station <: Endpoint
    token::String
    url::String
    CDO_Station(CDO_token::AbstractString, station::AbstractString) =
        new(CDO_token, string("https://www.ncdc.noaa.gov/cdo-web/api/v2/stations/", station))
end
struct CDO_Stations <: Endpoint
    token::String
    url::String
    CDO_Stations(CDO_token::AbstractString,
                 datasets::Union{AbstractString, AbstractVector{<:AbstractString}},
                 locations::Union{AbstractString, AbstractVector{<:AbstractString}},
                 datacategories::Union{AbstractString, AbstractVector{<:AbstractString}},
                 datatypes::Union{AbstractString, AbstractVector{<:AbstractString}},
                 extent::Union{AbstractVector{<:AbstractFloat}},
                 startdate::Date,
                 enddate::Date) =
        new(CDO_token,
            string("https://www.ncdc.noaa.gov/cdo-web/api/v2/stations?",
                   id_or_chain(datasets, "dataset"),
                   id_or_chain(locations, "location"),
                   id_or_chain(datacategories, "datacategory"),
                   id_or_chain(datatypes, "datatype"),
                   "extend=",
                   join(extent, ","),
                   "&startdate=",
                   startdate,
                   "&enddate=",
                   enddate,
                   "&limit=1000&offset=1"))
end

# Constants
const CDO_Single = Union{CDO_DataCategory, CDO_Dataset, CDO_DataType, CDO_LocationCategory, CDO_Location, CDO_Station}
const CDO_Meta = Union{CDO_Data, CDO_DataCategories, CDO_Datasets, CDO_DataTypes, CDO_LocationCategories, CDO_Locations, CDO_Stations}
const DATASETS = ["GHCND", "GSOM", "GSOY", "NEXRAD2", "NEXRAD3", "NORMAL_ANN",
                  "NORMAL_DLY", "NORMAL_HLY", "NORMAL_MLY", "PRECIP_15", "PRECIP_HLY"]
const LOCATIONCATEGORIES = ["CITY", "CLIM_DIV", "CLIM_REG", "CNTRY", "CNTY", "HYD_ACC",
                            "HYD_CAT", "HYD_REG", "HYD_SUB", "ST", "US_TERR", "ZIP"]
const DATACATEGORIES = ["ANNAGR", "ANNDD", "ANNPRCP", "ANNTEMP", "AUAGR", "AUDD",
                        "AUPRCP", "AUTEMP", "COMP", "COMPAGR", "DD", "DUALPOLMOMENT",
                        "ECHOTOP", "EVAP", "HYDROMETEOR", "LAND", "MISC", "OTHER", "OVERLAY",
                        "PRCP", "PRES", "REFLECTIVITY", "SKY", "SPAGR", "SPDD", "SPPRCP", "SPTEMP",
                        "SUAGR", "SUDD", "SUN", "SUPRCP", "SUTEMP", "TEMP", "VELOCITY", "VERTINTLIQUID",
                        "WATER", "WIAGR", "WIDD", "WIND", "WIPRCP", "WITEMP", "WXTYPE"]
const CDO_NonValidToken = ArgumentError("The CDO_token is not valid.")
const CDO_NonValidDatasets = ArgumentError("The provided datasets argument is not valid. Valid datasets are: " * reduce(*, DATASETS .* ", ")[1:end - 2] * ".")
const CDO_NonValidDataCategories = ArgumentError("The data category argument is not valid. Valid data categories are: " * reduce(*, DATACATEGORIES .* ", ")[1:end - 2] * ".")
const CDO_NonValidDataTypes = ArgumentError("The provided datatypes argument is not valid. A valid data type has form r\"^[A-Z0-9-]{3,}\$\"")
const CDO_NonValidLocationCategories = ArgumentError("The provided location categories argument is not valid. Valid location categories are: CITY, CLIM_DIV, CLIM_REG, CNTRY, CNTY, HYD_ACC, HYD_CAT, HYD_REG, HYD_SUB, ST, US_TERR, ZIP.")
const CDO_NonValidLocations = ArgumentError("The provided locations argument is not valid. A valid location has form r\"^(CITY|CLIM_DIV|CLIM_REG|CNTRY|CNTY|FIPS|HYD_ACC|HYD_CAT|HYD_REG|HYD_SUB|ST|US_TERR|ZIP):[A-Z0-9]{1,}\$\"")
const CDO_NonValidStations = ArgumentError("The provided locations argument is not valid. A valid location has form r\"(COOP|GHCND|NEXRAD|WBAN):[a-zA-Z0-9-]{1,}\$\"")
const CDO_NonValidExtent = ArgumentError("Extent must be of length 0 or 4.")

const Names = [(:date, :datatype, :station, :attributes, :value),
               (:name, :id),
               (:mindate, :maxdate, :name, :datacoverage, :id),
               (:uid, :mindate, :maxdate, :name, :datacoverage, :id),
               (:mindate, :maxdate, :datacoverage, :id),
               (:mindate, :maxdate, :name, :datacoverage, :id),
               (:elevation, :mindate, :maxdate, :latitude, :name, :datacoverage, :id, :elevationUnit, :longitude)
               ]
const Types = [Tuple{DateTime, String, String, Union{String, Missing}, Float64},
               Tuple{String, String},
               Tuple{Date, Date, String, Union{Float64, Missing}, String},
               Tuple{String, Date, Date, String, Float64, String},
               Tuple{Date, Date, Union{Float64, Missing}, String},
               Tuple{Date, Date, String, Float64, String},
               Tuple{Union{Float64, Missing}, Date, Date, Float64, String, Float64, String, Union{String, Missing}, Float64}
               ]

"""
    schema(::Endpoint)::Vector{<:NamedTuple}
    schema(::Endpoint, obj)::NamedTuple

Applies schema
"""
schema(::Endpoint) = error("names is not defined for $(typeof(obj)).")
schema(::CDO_Data) = Vector{NamedTuple{Names[1],Types[1]}}()
schema(::CDO_Data, obj) = (date = DateTime(obj.date),
                           datatype = convert(String, obj.datatype),
                           station = convert(String, obj.station),
                           attributes = get(obj, "attributes", missing),
                           value = obj.value)
schema(::CDO_DataCategory) = Vector{NamedTuple{Names[2],Types[2]}}()
schema(::CDO_DataCategory, obj) = (name = convert(String, obj.name),
                                   id = convert(String, obj.id))
schema(::CDO_DataCategories) = Vector{NamedTuple{Names[2],Types[2]}}()
schema(::CDO_DataCategories, obj) = (name = convert(String, obj.name),
                                     id = convert(String, obj.id))
schema(::CDO_Dataset) = Vector{NamedTuple{Names[3],Types[3]}}()
schema(::CDO_Dataset, obj) = (mindate = Date(obj.mindate),
                              maxdate = Date(obj.maxdate),
                              name = convert(String, obj.name),
                              datacoverage = get(obj, "datacoverage", missing) |>
                                  (x -> isa(x, Missing) ? x : convert(Float64, x)),
                              id = convert(String, obj.id))
schema(::CDO_Datasets) = Vector{NamedTuple{Names[4],Types[4]}}()
schema(::CDO_Datasets, obj) = (uid = convert(String, obj.uid),
                               mindate = Date(obj.mindate),
                               maxdate = Date(obj.maxdate),
                               name = convert(String, obj.name),
                               datacoverage = convert(Float64, obj.datacoverage),
                               id = convert(String, obj.id))
schema(::CDO_DataType) = Vector{NamedTuple{Names[5],Types[5]}}()
schema(::CDO_DataType, obj) = (mindate = Date(obj.mindate),
                               maxdate = Date(obj.maxdate),
                               datacoverage = get(obj, "datacoverage", missing) |>
                                   (x -> isa(x, Missing) ? x : convert(Float64, x)),
                               id = convert(String, obj.id))
schema(::CDO_DataTypes) = Vector{NamedTuple{Names[3],Types[3]}}()
schema(::CDO_DataTypes, obj) = (mindate = Date(obj.mindate),
                                maxdate = Date(obj.maxdate),
                                name = convert(String, obj.name),
                                datacoverage = get(obj, "datacoverage", missing) |>
                                    (x -> isa(x, Missing) ? x : convert(Float64, x)),
                                id = convert(String, obj.id))
schema(::CDO_LocationCategory) = Vector{NamedTuple{Names[2],Types[2]}}()
schema(::CDO_LocationCategory, obj) = (name = convert(String, obj.name),
                                       id = convert(String, obj.id))
schema(::CDO_LocationCategories) = Vector{NamedTuple{Names[2],Types[2]}}()
schema(::CDO_LocationCategories, obj) = (name = convert(String, obj.name),
                                         id = convert(String, obj.id))
schema(::CDO_Location) = Vector{NamedTuple{Names[6],Types[6]}}()
schema(::CDO_Location, obj) = (mindate = Date(obj.mindate),
                               maxdate = Date(obj.maxdate),
                               name = convert(String, obj.name),
                               datacoverage = convert(Float64, obj.datacoverage),
                               id = convert(String, obj.id))
schema(::CDO_Locations) = Vector{NamedTuple{Names[6],Types[6]}}()
schema(::CDO_Locations, obj) = (mindate = Date(obj.mindate),
                                maxdate = Date(obj.maxdate),
                                name = convert(String, obj.name),
                                datacoverage = convert(Float64, obj.datacoverage),
                                id = convert(String, obj.id))
schema(::CDO_Station) = Vector{NamedTuple{Names[7],Types[7]}}()
schema(::CDO_Station, obj) = (elevation = get(obj, "elevation", missing) |>
                                  (x -> isa(x, Missing) ? x : convert(Float64, x)),
                              mindate = Date(obj.mindate),
                              maxdate = Date(obj.maxdate),
                              latitude = convert(Float64, obj.latitude),
                              name = convert(String, obj.name),
                              datacoverage = convert(Float64, obj.datacoverage),
                              id = convert(String, obj.id),
                              elevationUnit = get(obj, "elevationUnit", missing) |>
                                  (x -> isa(x, Missing) ? x : convert(String, x)),
                              longitude = convert(Float64, obj.longitude))
schema(::CDO_Stations) = Vector{NamedTuple{Names[7],Types[7]}}()
schema(::CDO_Stations, obj) = (elevation = get(obj, "elevation", missing) |>
                                  (x -> isa(x, Missing) ? x : convert(Float64, x)),
                               mindate = Date(obj.mindate),
                               maxdate = Date(obj.maxdate),
                               latitude = convert(Float64, obj.latitude),
                               name = convert(String, obj.name),
                               datacoverage = convert(Float64, obj.datacoverage),
                               id = convert(String, obj.id),
                               elevationUnit = get(obj, "elevationUnit", missing) |>
                                   (x -> isa(x, Missing) ? x : convert(String, x)),
                               longitude = convert(Float64, obj.longitude))

"""
    exception(::Endpoint)::DataFrame
    exception(::Endpoint)::ArgumentError(msg::AbstractString)

In case of no results for a data request, return either an empty DataFrame of
the expected type indicating no data for the query parameters or an ArgumentError
if it is the result can be deduced to be an error in the argument.
"""
exception(obj::Endpoint) = DataFrame(types(obj), names(obj), 0)
exception(obj::CDO_DataCategory) = ArgumentError(obj.url[57:end] * " is not a valid data category. For a complete list of valid data categories run `cdo_datacategories(CDO_token::AbstractString)`.")
exception(obj::CDO_Dataset) = ArgumentError(obj.url[51:end] * " is not a valid dataset. For a complete list of valid datasets run `cdo_datasets(CDO_token::AbstractString)`.")
exception(obj::CDO_DataType) = ArgumentError(obj.url[52:end] * " is not a valid data type. For a complete list of valid data types run `cdo_datatypes(CDO_token::AbstractString)`.")
exception(obj::CDO_LocationCategory) = ArgumentError(obj.url[61:end] * " is not a valid location category. For a complete list of valid location categories run `cdo_locationcategories(CDO_token::AbstractString)`.")
exception(obj::CDO_Location) = ArgumentError(obj.url[52:end] * " is not a valid location. For a complete list of valid locations run `cdo_locations(CDO_token::AbstractString)`.")
exception(obj::CDO_Station) = ArgumentError(obj.url[51:end] * " is not a valid weather station. For a complete list of valid stations run `cdo_stations(CDO_token::AbstractString)`.")

"""
    skeleton(endpoint::CDO_Single, jsontext::AbstractString)

Return a dataframe with the parse values of a JSON according to the endpoint.
In case of no observations an empty dataframe is returned.
If an empty JSON indicates an ArgumentError this one is throw.
"""
function skeleton(endpoint::CDO_Single, jsontext::AbstractString)
    jsontext == "{}" && return exception(endpoint)
    # output = schema(endpoint)
    DataFrame([schema(endpoint, value(jsontext))])
    # foreach(x -> push!(output, schema(endpoint, x)), value(jsontext))
    # DataFrame(output)
end

"""
    parse(obj::Endpoint)

Return the results for the endpoint.
"""
parse(obj::Endpoint) = error("parse is not defined for $(typeof(obj)).")
function parse(obj::CDO_Single)
    response = request("GET", obj.url, [("token", obj.token)])
    jsontext = String(response.body)
    skeleton(obj, jsontext)
end
function parse(obj::CDO_Meta)
    header = [("token", obj.token)]
    url = obj.url
    response = request("GET", url, header)
    jsontext = String(response.body)
    jsontext == "{}" && return DataFrame(Ts, Ns, 0)
    json = value(jsontext)
    Count = Int64(json["metadata"]["resultset"]["count"])
    urls = string.(url[1:end - 1], (1:1000:Count)[2:end])
    output = schema(obj)
    foreach(x -> push!(output, schema(obj, x)), json["results"])
    for suburl ∈ urls
        response = request("GET", suburl, header)
        jsontext = String(response.body)
        foreach(x -> push!(output, schema(obj, x)), value(jsontext)["results"])
    end
    DataFrame(output)
end
function parse(obj::CDO_Data)
    header = [("token", obj.token)]
    urls = obj.url
    output = schema(obj)
    for url ∈ urls
        response = request("GET", url, header)
        jsontext = String(response.body)
        jsontext == "{}" && return DataFrame(output)
        json = value(jsontext)
        Count = Int64(json["metadata"]["resultset"]["count"])
        urls = url[1:end - 1] .* string.(1:1000:Count)[2:end] .* "&includemetadata=false"
        foreach(x -> push!(output, schema(obj, x)), json["results"])
        for suburl ∈ urls
            response = request("GET", suburl, header)
            jsontext = String(response.body)
            foreach(x -> push!(output, schema(obj, x)), value(jsontext)["results"])
        end
    end
    DataFrame(output)
end


text:
println("Testing particle filter...")
# Bearing-only tracking example from Gordon et al. 1993
Phi = eye(4)
Phi[1, 2] = Phi[3, 4] = 1
Gamma = [0.5 0; 1 0; 0 0.5; 0 1]
q = 0.001
r = 0.005
x0 = [-0.05, 0.01, 0.7, -0.055]

process(x, t) = Phi * x + Gamma * q * randn(2)
bearing(x) = atan(x[2] / x[1])
obs_loglik(x, y, t) = logpdf(Normal(bearing(x), r), y[1])

mod = NonlinearSSM(process, obs_loglik, 4, 1, 1)
filt = ParticleFilter(4, 1000, () -> zeros(4))

nt = 24
srand(1234)
xx = simulate(mod, nt, x0)
yy = (bearing(xx) + randn(nt) * r)'

# using PyPlot
# subplot(111, aspect="equal")
# plot(xx[1, :]', xx[3, :]')
# scatter(0, 0)

x0_prior = MvNormal([0, 0, 0.4, -0.05], diagm([0.5, 0.005, 0.3, 0.01].^2))
ensemble_0 = rand(x0_prior, filt.nparticles)

ensemble = predict(mod, ensemble_0)

observe(mod, ensemble[:, 1], yy[:, 1])
observe(mod, ensemble, yy[:, 1])

ensemble_new = update(mod, ensemble, yy[:, 1], filt)

# scatter(xx[1,1], xx[3,1])
# plot(ensemble[1,:]', ensemble[3, :]', "k+")
# plot(ensemble_new[1,:]', ensemble_new[3, :]', "r.")

fs = filter(mod, yy, ensemble_0, filt)

fs_mean = mean(fs, 2)
# plot(fs_mean[1,:]', fs_mean[3,:]')

println("Particle filter passed.")


text:
# Abstract interface for diffusion penalties

"""
  Abstract type for the different diffusion penalties that can be expressed
  in the framework of Yan et.al. "Interior Penalties for Summation-by-Parts
  Discretizations of Linear Second-Order Differential Equations".

  Each type must implement a method of:

  ```
   applyPenalty(penalty::AbstractDiffusionPenalty, sbp, sbpface,
                diffusion::AbstractDiffusion, iface::Interface,
                delta_w::AbstractMatrix, theta::AbstractMatrix,
                wL::AbstractMatrix, wR::AbstractMatrix,
                nrm_face::AbstractMatrix,
                alphas::AbstractVector,
                jacL::AbstractVector, jacR::AbstractVector,
                res1L::AbstractMatrix, res1R::AbstractMatrix,
                res2L::AbstractMatrix, res2R::AbstractMatrix
               ) where {Tsol, Tres}

  function applyDirichletPenalty(penalty::AbstractDiffusionPenalty,
                      sbp, sbpface,diffusion::AbstractDiffusion,
                      bndry::Boundary, delta_w::AbstractMatrix
                      wL::AbstractMatrix, nrm_face::AbstractMatrix,
                      alpha::Number, jacL::AbstractVector,
                      res1L::AbstractMatrix)

  ```
  specializing the `penalty` argument.  See the docstrings for these
  functions for details.

"""
abstract type AbstractDiffusionPenalty end

"""
  Applies the interface penalty for 2nd derivative terms.
  More specifically, it applies the penalty matrix
  for both sides of the face at the same time:

  [res1L   = [T1 T3  [delta_w   and   [res1R   = [T1 T3  [-delta_w
   res2L]     T2 T4]  theta]           res2R]     T2 T4]  theta]

  and similarly for the right element.  Note that delta_w is negated for
  the right element, and this function performs the required negation.

  **Inputs**

   * penalty: the [`AbstractDiffusionPenalty`](@ref) object
   * sbp
   * params
   * sbpface
   * diffusion: the [`AbstractDiffusion`](@ref) object
   * iface: `Interface` object
   * delta_w: difference of entropy variables at the face, as compuated by
              [`getFaceVariables`](@ref)
   * theta: Dgk * wL + Dgn * wR, as compuated by `getFaceVariables`
   * qL: conservative variables for left element, `numDofPerNode` x
         `numNodesPerElement`
   * qR: conservative variables for right element`, same size as `qL`
   * wL: entropy variables for the left element, `numDofPerNode` x
         `numNodesPerElement`
   * wR: entropy variables for the right element, same size as `wR`
   * coordsL: xyz coordinates of volume nodes of left element, `dim` x
              `numNodesPerElement`
   * coordsR: xyz coordinates of volumke nodes of right element, same size
              as `coordsL`
   * nrm_face: (scaled) normal vectors at the face, `dim` x `numNodesPerFace`
   * alphas: vector of length 2 containing alpha_gk and alpha_gn
   * dxidxL: (scaled) mapping jacobian for left element, `dim` x `dim` x
             `numNodesPerElement`
   * dxidxR: (scaled) mapping jacobian for right element, same size as `dxidxL`
   * jacL: mapping jacobian determinant for the left element,
           `numNodesPerElement`
   * jacR: mapping jacobian determinant for the right element, same size as
           `jacR`

  **Inputs/Outputs**

   * res1L: `numDofPerNode` x `numNodesPerFace`
   * res1R: same size as above
   * res2L: same size as above
   * res2R: same size as above

  All output arrays are overwritten
"""
function applyPenalty(penalty::AbstractDiffusionPenalty, sbp,
                      params::AbstractParamType, sbpface,
                      diffusion::AbstractDiffusion, iface::Interface,
                      delta_w::AbstractMatrix, theta::AbstractMatrix,
                      qL_el::AbstractMatrix, qR_el::AbstractMatrix,
                      wL::AbstractMatrix, wR::AbstractMatrix,
                      coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      alphas::AbstractVector,
                      dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                      jacL::AbstractVector, jacR::AbstractVector,
                      res1L::AbstractMatrix, res1R::AbstractMatrix,
                      res2L::AbstractMatrix, res2R::AbstractMatrix
                     )

  error("generic fallback for applyPenalty() reached: did you forget to extend it with a new method for your AbstractDiffusionPenalty?")
end

"""
  **Inputs**

   * penalty: the [`AbstractDiffusionPenalty`](@ref) object
   * sbp
   * sbpface
   * diffusion: the [`AbstractDiffusion`](@ref) object
   * iface: `Interface` object
   * delta_w: difference of entropy variables at the face, as compuated by
              [`getFaceVariables`](@ref)
   * delta_w_dotL: array containing derivative of `delta_w` wrt `qL`,
                   `numDofPerNode` x `numDofPerNode` x `numNodesPerFace`
                   x `numNodesPerElement`
   * delta_w_dotR: derivative of `delta_w` wrt `qR, same shape as `delta_w_dotR`
   * theta: Dgk * wL + Dgn * wR, as compuated by `getFaceVariables`
   * theta_dotL derivative of `theta` wrt `qL`, same shape as `delta_w_dotL`
   * theta_dotR derivative of `theta` wrt `qR`, same shape as `delta_w_dotL`
   * qL: conservative variables for the left element, `numDofPerNode` x
         `numNodesPerElement`
   * qR: conservative variables for the right element, same size as `qL`
   * wL: entropy variables for the left element, `numDofPerNode` x
         `numNodesPerElement`
   * wR: entropy variables for the right element, same size as `wL`
   * coordsL: xyz coordinates of volume nodes of left element, `dim` x
              `numNodesPerElement`
   * coordsR: xyz coordinates of the volume nodes of the right element,
              same size as `coordsL`
   * nrm_face: (scaled) normal vectors at the face, `dim` x `numNodesPerFace`
   * alphas: vector of length 2 containing alpha_gk and alpha_gn
   * dxidxL: (scaled) mapping jacobian for left element, `dim` x `dim`,
             x `numNodesPerElement
   * dxidxR: (scaled) mapping jacobian for right element, same size as
             `dxidxL`
   * jacL: mapping jacobian determinant for the left element,
           `numNodesPerElement`
   * jacR: mapping jacobian determinant for the right element, same size as
           `jacR`

  **Inputs/Outputs**

   * res1L_dotL: derivative of `res1L` wrt `qL`, same size as `delta_w_dotL`
   * res1L_dotR: derivative of `res1L` wrt `qR`, same size as above`
   * res1R_dotL: derivative of `res1R` wrt `qL`, same size as `delta_w_dotL`
   * res1R_dotR: derivative of `res1R` wrt `qR`, same size as above`
   * res2L: see `applyPenalty`
   * res2R: see `applyPenalty
   * res2L_dotL: derivative of `res2L` wrt `qL`, same size as `delta_w_dotL`
   * res2L_dotR: derivative of `res2L` wrt `qR`, same size as above`
   * res2R_dotL: derivative of `res2R` wrt `qL`, same size as `delta_w_dotL`
   * res2R_dotR: derivative of `res2R` wrt `qR`, same size as above`

   All of these array are overwritten

  **Outputs**

   * T4_nonzero: true if the T4 penalty is zero, false otherwise. This is
                 used to determine the sparsity pattern of the Jacobian.
"""
function applyPenalty_diff(penalty::AbstractDiffusion, sbp,
                      params::AbstractParamType, sbpface,
                      diffusion::AbstractDiffusion, iface::Interface,
                      delta_w::AbstractMatrix,
                      delta_w_dotL::Abstract4DArray, delta_w_dotR::Abstract4DArray,
                      theta::AbstractMatrix,
                      theta_dotL::Abstract4DArray, theta_dotR::Abstract4DArray,
                      qL::AbstractMatrix, qR::AbstractMatrix,
                      wL::AbstractMatrix, wR::AbstractMatrix,
                      coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      alphas::AbstractVector,
                      dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                      jacL::AbstractVector, jacR::AbstractVector,
                      res1L_dotL::Abstract4DArray, res1L_dotR::Abstract4DArray,
                      res1R_dotL::Abstract4DArray, res1R_dotR::Abstract4DArray,
                      res2L::AbstractMatrix, res2R::AbstractMatrix,
                      res2L_dotL::Abstract4DArray, res2L_dotR::Abstract4DArray,
                      res2R_dotL::Abstract4DArray, res2R_dotR::Abstract4DArray)

  error("generic fallback for applyPenalty_diff() reached: did you forget to extend it with a new method for your AbstractDiffusionPenalty?")
end

"""
  Applies the Dirichlet penalty matrix T_gammaD.

  **Inputs**

   * penalty: [`AbstractDiffusionPenalty`](@ref)
   * sbp
   * params
   * sbpface
   * diffusion: [`AbstractDiffusion`](@ref)
   * bndry: `Boundary` object
   * delta_w: values to multiply the penalty against, `numDofPerNode` x
              `numNodesPerFace`
   * qL: entropy variables for the element, `numDofPerNode` x `numNodesPerElement`
   * wL: entropy variables for the element, `numDofPerNode` x `numNodesPerElement`
   * coordsL: xyz coordinates of volume nodes of the element, `dim` x
              `numNodesPerElement`
   * nrm_face: (scaled) normal vector at each face node, `dim` x
               `numNodesPerFace`
   * alpha: the alpha_gk parameter (Number)
   * dxidxL: (scaled) mapping jacobian for the element, `dim` x `dim` x
             `numNodesPerElement`
   * jacL: mapping jacobian determinant, `numNodesPerElement`

  **Inputs/Outputs**

   * res1L: `numDofPerNode` x `numNodesPerFace` array to overwrite with
            the result

"""
function applyDirichletPenalty(penalty::AbstractDiffusionPenalty, sbp,
                      params::AbstractParamType, sbpface,
                      diffusion::AbstractDiffusion, bndry::Boundary,
                      delta_w::AbstractMatrix,
                      qL::AbstractMatrix,
                      wL::AbstractMatrix,
                      coordsL::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      alpha::Number,
                      dxidxL::Abstract3DArray,
                      jacL::AbstractVector,
                      res1L::AbstractMatrix)

  error("generic fallback for applDirichletPenalty() reached: did you forget to extend it with a new method for your AbstractDiffusionPenalty?")

end

#------------------------------------------------------------------------------
# Diffusion penalties

include("br2_penalty.jl")
include("inviscid_penalty.jl")

global const DiffusionPenaltyDict = Dict{String, Type{T} where T <: AbstractDiffusionPenalty}(
"BR2" => BR2Penalty,
"Inviscid" => InviscidPenalty,
)

"""
  Returns a fully constructed [`AbstractDiffusionPenalty`](@ref).

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * name: name of penalty to get, defaults to opts["DiffusionPenalty"]

  **Outputs**

   * obj: `AbstractDiffusionPenalty` object
"""
function getDiffusionPenalty(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts,
                             name::String=opts["DiffusionPenalty"]
                            ) where {Tsol, Tres}

  obj = DiffusionPenaltyDict[name]{Tsol, Tres}(mesh, sbp, opts)
  assertArraysUnique(obj)

  return obj
end


text:
##
##  i n t e r p o l a t e . j l  Interpolation
##

function interp1d{T<:Real}(xi::Array{T,1}, yi::Array{T,1}, x::T;
                           method::Symbol = :cubic)
    n = length(xi)
    if length(yi) != n
        error("Vectors 'xi', 'yi' must be of same length.")
    end
    if n == 1
       error("At least two data points must be provided.")
    end
    if any(xi[1:end-1] .> xi[2:end])
        error("Elements in 'xi' must be strictly increasing.")
    end

    local k::Int = 0
    if x < xi[1]; k = 0
    elseif x >= xi[end]; k = n
    else
        k = maximum(find(xi[1:(end-1)] .<= x))
    end

    if method == :constant
        if k == 0; k = 1; end
        return yi[k]

    elseif method == :nearest
        if k == 0; k = 1
        elseif k < n
            if x >= (xi[k] + xi[k+1]) / 2.0
                k += 1
            end
        end
        return yi[k]

    elseif method == :linear
        if k == 0; k = 1
        elseif k == n; k = n-1; end
        s = x - xi[k]
        delta = (yi[k+1] - yi[k]) / (xi[k+1] - xi[k])
        return yi[k] + (x - xi[k])*delta

    elseif method == :spline
        error("Method $(method) is not yet implemented.")

    elseif method == :cubic
        return pchip(xi, yi, x)

    else
        error("Unknown interpolation method: $(method)")
    end
end

function interp1d{T<:Real}(xi::Array{T,1}, yi::Array{T,1}, x::Array{T,1};
                           method::Symbol = :cubic)
    local m::Int = length(x)
    y = zeros(m)
    for i = 1:m
        y[i] = interp1d(xi, yi, x[i], method = method)
    end
    return y
end

#-- Moler: Piecewise Cubic Hermitean Interpolation Polynomials
#
function pchip{T<:Real}(xi::Array{T,1}, yi::Array{T,1}, x::T)
    n = length(xi)
    if length(yi) != n
        error("Vectors 'xi', 'yi' must be of same length.")
    end
    if any(xi[1:end-1] .> xi[2:end])
        error("Elements in 'xi' must be strictly increasing.")
    end

    # First derivatives
    h = diff(xi)
    delta = diff(yi) ./ h
    d = pchip_slopes(h, delta)

    # Piecewise polynomial coefficients
    a = (3.0*delta .- 2.0*d[1:(n-1)] .- d[2:n]) ./ h
    b = (d[1:(n-1)] .- 2.0*delta .+ d[2:n]) ./ h.^2

    # Find subinterval index k so that xi[k] <= x < xi[k+1]
    if x < xi[1]
        k = 1
    else
        k = maximum(find(xi[1:(n-1)] .<= x))
    end

    # Evaluate interpolant
    s = x - xi[k]
    v = yi[k] + s*(d[k] + s*(a[k] + s*b[k]))

    return v
end

function pchip_slopes(h, delta)

    # Slopes at interior points
    local n = length(h) + 1
    d = zeros(n)
    k = find(sign(delta[1:(n-2)]) .* sign(delta[2:(n-1)]) .> 0.0) .+ 1
    w1 = 2.0*h[k] .+ h[k.-1]
    w2 = h[k] .+ 2.0*h[k.-1]
    d[k] = (w1+w2) ./ (w1./delta[k.-1] + w2./delta[k])

    # Slopes at endpoints
    d[1] = pchip_end(h[1], h[2], delta[1], delta[2])
    d[n] = pchip_end(h[n-1], h[n-2], delta[n-1], delta[n-2])

    return d
end

function pchip_end(h1, h2, del1, del2)

    # Noncentered, shape-preserving, three-point formula.
    d = ((2.0*h1 + h2)*del1 - h1*del2) / (h1 + h2)
    if sign(d) != sign(del1)
        d = 0
    elseif sign(del1) != sign(del2)  &&  abs(d) > abs(3*del1)
        d = 3.0*del1
    end

    return d
end


text:
using Test

include("leap.jl")

@testset "Year not divisible by 4: common year" begin
    @test !is_leap_year(2015)
end

@testset "Year divisible by 2, not divisible by 4: common year" begin
    @test !is_leap_year(1970)
end

@testset "Year divisible by 4, not divisible by 100: leap year" begin
    @test is_leap_year(1996)
end

@testset "Year divisible by 4 and 5: leap year" begin
    @test is_leap_year(1960)
end

@testset "Year divisible by 4, not divisible by 100: leap year" begin
    @test is_leap_year(2016)
end

@testset "Year divisible by 100, not divisible by 400: common year" begin
    @test !is_leap_year(2100)
end

@testset "Year divisible by 100 but not by 3: common year" begin
    @test !is_leap_year(1900)
end

@testset "Year divisible by 400: leap year" begin
    @test is_leap_year(2000)
end

@testset "Year divisible by 400 but not by 125: leap year" begin
    @test is_leap_year(2400)
end

@testset "Year divisible by 200, not divisible by 400: common year" begin
    @test !is_leap_year(1800)
end


text:
export NFFTOp
import Base.adjoint

mutable struct NFFTOp{T} <: AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod! :: Function
  tprod! :: Nothing
  ctprod! :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  args5 :: Bool
  use_prod5! :: Bool
  allocated5 :: Bool
  Mv5 :: Vector{T}
  Mtu5 :: Vector{T}
  plan
  toeplitz :: Bool
end

"""
    NFFTOp(shape::Tuple, tr::Trajectory; kargs...)
    NFFTOp(shape::Tuple, tr::AbstractMatrix; kargs...)

generates a `NFFTOp` which evaluates the MRI Fourier signal encoding operator using the NFFT.

# Arguments:
* `shape::NTuple{D,Int64}`  - size of image to encode/reconstruct
* `tr`                      - Either a `Trajectory` object, or a `ND x Nsamples` matrix for an ND-dimenensional (e.g. 2D or 3D) NFFT with `Nsamples` k-space samples
* (`nodes=nothing`)         - Array containg the trajectory nodes (redundant)
* (`kargs`)                 - additional keyword arguments
"""
function NFFTOp(shape::Tuple, tr::AbstractMatrix{T}; toeplitz=false, oversamplingFactor=1.25, kernelSize=3, kargs...) where {T}
  # AMM : adding tr = CuArray(tr);
  # tr = CuArray(tr);
  # AMM: Temp: I am changing hte oversamplingFactor from 1.25(default) to 1.3
  oversamplingFactor=1.3
  # @infiltrate
  plan = plan_nfft(tr, shape, m=kernelSize, σ=oversamplingFactor, precompute=NFFT.FULL)
  # AMM: This is if I want to run CuNFFT
  # plan = plan_nfft(typeof(tr), tr, shape, m=kernelSize, σ=oversamplingFactor, precompute=NFFT.FULL)
  # AMM: Temp: Infiltrating
  # @infiltrate
  return NFFTOp{Complex{T}}(size(tr,2), prod(shape), false, false
            , (res,x) -> (res .= produ(plan,x))
            , nothing
            , (res,y) -> (res .= ctprodu(plan,y))
            , 0, 0, 0, false, false, false, Complex{T}[], Complex{T}[]
            , plan, toeplitz)
end

function NFFTOp(shape::Tuple, tr::Trajectory; toeplitz=false, oversamplingFactor=1.25, kernelSize=3, kargs...)
  return NFFTOp(shape, kspaceNodes(tr); toeplitz=toeplitz, oversamplingFactor=oversamplingFactor, kernelSize=kernelSize, kargs...)
end

function produ(plan::NFFT.NFFTPlan, x::Vector{T}) where T<:Union{Real,Complex}
  y = plan * reshape(x[:],plan.N)
  return vec(y)
end

function ctprodu(plan::NFFT.NFFTPlan, y::Vector{T}) where T<:Union{Real,Complex}
  x = adjoint(plan) * y[:]
  return vec(x)
end

function Base.copy(S::NFFTOp{T}) where {T}
  plan = copy(S.plan)
  return NFFTOp{T}(size(plan.x,2), prod(plan.N), false, false
              , (res,x) -> (res .= produ(plan,x))
              , nothing
              , (res,y) -> (res .= ctprodu(plan,y))
              , 0, 0, 0, false, false, false, T[], T[]
              , plan, S.toeplitz)
end

#########################################################################
### Toeplitz Operator ###
#########################################################################
struct NFFTToeplitzNormalOp{T,D,W}
  shape::NTuple{D,Int}
  weights::W
  fftplan
  ifftplan
  λ::Array{T}
  xL1::Array{T,D}
  xL2::Array{T,D}
end

function NFFTToeplitzNormalOp(S::NFFTOp{T}, W::UniformScaling=I) where {T}
  shape = S.plan.N

  # plan the FFTs
  fftplan  = plan_fft( zeros(T, 2 .* shape);flags=FFTW.MEASURE)
  ifftplan = plan_ifft(zeros(T, 2 .* shape);flags=FFTW.MEASURE)

  λ = calculateToeplitzKernel(shape, S.plan.x; m = S.plan.params.m, σ = S.plan.params.σ, window = S.plan.params.window, LUTSize = S.plan.params.LUTSize, fftplan = fftplan)

  xL1 = Array{T}(undef, 2 .* shape)
  xL2 = similar(xL1)

  return NFFTToeplitzNormalOp(shape, W, fftplan, ifftplan, λ, xL1, xL2)
end

function SparsityOperators.normalOperator(S::NFFTOp, W::UniformScaling=I)
  if S.toeplitz
    return NFFTToeplitzNormalOp(S,W)
  else
    return NormalOp(S,W)
  end
end

function SparsityOperators.normalOperator(S::NFFTOp, W)
  if S.toeplitz
    @warn "Topelitz with non-Uniform scaling is currently not implemented. Using a non-Toeplitz implemenation instead."
  end
  return NormalOp(S,W)
end

function LinearAlgebra.mul!(y, S::NFFTToeplitzNormalOp, b)
  S.xL1 .= 0
  b = reshape(b, S.shape)

  S.xL1[CartesianIndices(b)] .= b
  mul!(S.xL2, S.fftplan, S.xL1)
  S.xL2 .*= S.λ
  mul!(S.xL1, S.ifftplan, S.xL2)

  y .= vec(S.xL1[CartesianIndices(b)])
  return y
end

Base.:*(S::NFFTToeplitzNormalOp, b::AbstractVector) = mul!(similar(b), S, b)
Base.size(S::NFFTToeplitzNormalOp) = S.shape
Base.size(S::NFFTToeplitzNormalOp, dim) = S.shape[dim]
Base.eltype(::Type{NFFTToeplitzNormalOp{T,D,W}}) where {T,D,W} = T

function Base.copy(A::NFFTToeplitzNormalOp{T,D,W}) where {T,D,W}
  fftplan  = plan_fft( zeros(T, 2 .* A.shape); flags=FFTW.MEASURE)
  ifftplan = plan_ifft(zeros(T, 2 .* A.shape); flags=FFTW.MEASURE)
  return NFFTToeplitzNormalOp(A.shape, A.weights, fftplan, ifftplan, A.λ, A.xL1, A.xL2)
end


text:
#### Kahan matrix

export Kahan
"""
[`Kahan` matrix](http://math.nist.gov/MatrixMarket/data/MMDELI/kahan/kahan.html)

```julia
julia> A=Kahan(5,5,1,35)
5x5 Kahan{Int64,Int64}:
 1.0  -0.540302  -0.540302  -0.540302  -0.540302
 0.0   0.841471  -0.454649  -0.454649  -0.454649
 0.0   0.0        0.708073  -0.382574  -0.382574
 0.0   0.0        0.0        0.595823  -0.321925
 0.0   0.0        0.0        0.0        0.501368

julia> A=Kahan(5,3,0.5,0)
5x3 Kahan{Float64,Int64}:
 1.0  -0.877583  -0.877583
 0.0   0.479426  -0.420735
 0.0   0.0        0.229849
 0.0   0.0        0.0
 0.0   0.0        0.0

julia> A=Kahan(3,5,0.5,1e-3)
3x5 Kahan{Float64,Float64}:
 1.0  -0.877583  -0.877583  -0.877583  -0.877583
 0.0   0.479426  -0.420735  -0.420735  -0.420735
 0.0   0.0        0.229849  -0.201711  -0.201711
```

For more details see: N. Higham, A Survey of Condition Number Estimation
for Triangular Matrices, SIMAX, Vol. 29, No. 4, pp. 588, 1987,
http://eprints.ma.man.ac.uk/695/01/covered/MIMS_ep2007_10.pdf
"""
struct Kahan{T<:Number,T1<:Number} <: AbstractMatrix{T}
    m::Int # dimension
    n::Int # dimension
    theta::T # angle
    pert::T1 # perturbation is pert*eps()
end # immutable

# Define its size
size(A::Kahan, r::Int) = r==1 ? A.m : A.n
size(A::Kahan) = A.m, A.n

# Index into a Kahan
function getindex(A::Kahan,i::Integer,j::Integer)
    m=minimum(size(A))
    t=tan(A.theta)
    c=1.0/sqrt(1.0+t^2)
    s=t*c
    if i>m; return 0.0
    elseif i>j; return 0.0
    elseif i==j; return s^(i-1)+A.pert*eps()*(m-i+1)
    else return -c*s^(i-1)
    end
end # getindex

# Dense version of Kahan
Matrix(A::Kahan) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)]


text:
type GraphNode
    index::Int
end

type Graph
    nodes::Dict{Int, GraphNode}
    edges::Dict{Int, Array{Int, 1}}
    lengths::Dict{Tuple{Int, Int}, Real}
    max_index::Int
end

function init_graph()
    return Graph(Dict{Int, GraphNode}(), Dict{Int, Array{Int, 1}}(),
        Dict{Tuple{Int, Int}, Real}(), 1)
end

function num_nodes(g::Graph)
    return length(g.nodes)
end

function add_node!(g::Graph, index::Int)
    # Adds a node given its neighbors and edge lengths.
    new_node = GraphNode(index)
    g.nodes[index] = new_node
    g.edges[index] = Array(Int, 0)
    return index
end

function add_node!(g::Graph, index::Int, n::Array{Int, 1}, l::Array)
    # Adds a node given its neighbors and edge lengths.
    new_node = GraphNode(index)
    if index in keys(g.nodes)
        println("Warning: index already exists")
    end
    g.nodes[index] = new_node
    g.max_index += 1
    g.edges[index] = n
    for (node, length) in zip(n, l)
        push!(g.edges[node], index)
        g.lengths[(node, index)] = length
        g.lengths[(index, node)] = length
    end
    return index
end

function add_edge!(g::Graph, i1::Int, i2::Int, d::Real)
    push!(g.edges[i1], i2)
    push!(g.edges[i2], i1)
    g.lengths[(i1, i2)] = d
    g.lengths[(i2, i1)] = d
end

function all_shortest_paths(g::Graph)
    # returns the all-source shortest paths as a NxN matrix
    # Using Floyd-Warshall
    n = num_nodes(g)
    dist = fill(Inf, n, n)
    for (i1, i2) in keys(g.lengths)
        dist[i1, i2] = g.lengths[i1, i2]
    end
    for i in keys(g.nodes)
        dist[i, i] = 0
    end
    for k in range(1, n)
        for i in range(1, n)
            for j in range(1, n)
                if dist[i,j] > dist[i,k] + dist[k,j]
                    dist[i,j] = dist[i,k] + dist[k,j]
                end
            end
        end
    end
    return dist
end


text:
# for loops in Julia
# Navendu Pottekkat | 2020 | Apache-2.0 License
# https://www.apache.org/licenses/LICENSE-2.0

for i in 1:5
    print(i, ", ")
end


text:
using JuPOT

             m = 6 # No. of days data
             returns_data = ARGS[1]
			 returns = parse.(split(returns_data, ","))

             meanReturns_data =  ARGS[2] # Returns a matrix of size(n) with entries between 0-1s
			 meanReturns = parse.(split(meanReturns_data, ","))

             tickers_data = ARGS[3]
			 tickers = split(tickers_data, ",")

             n = length(tickers)

             matrix = reshape(returns,m,n)
             covariance = cov(matrix) #covariance matrix of returns

             # Assets data structure containing, names, expected returns, covarariance
             assets = AssetsCollection(tickers, meanReturns, covariance)

             target_ret_test = [0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80]

             function efficientFrontierMVO(assets, target_ret_test)

                 # Create the Arrays to be returned, and set them to be the length of the
                 # target_ret_test from the input.
                 n = length(target_ret_test)
                 risk = Array(Float64,n)
                 index = Array(Float64,n)

                 for i in 1:n
                     # Get the target return we want to test.
                     target_ret = target_ret_test[i]

                     # Setup and Optimize the SimpleMVO with the new Target Return.
                     mvo = SimpleMVO(assets, target_ret; short_sale=false)
                     result = optimize(mvo)

                     # Risk is the Objective Value. Index the Target Return.
                     risk[i] = result[1]
                     index[i] = target_ret
                 end

                 # Return as a Tuple of the Arrays
                 return risk, index
             end

             efd = efficientFrontierMVO(assets, target_ret_test)


text:
using AWS.EC2
using AWS

include("config.jl")

instances = ec2_launch(ec2_ami, ec2_sshkey, insttype=ec2_insttype, n=ec2_instnum, clustername=ec2_clustername)
instances

hosts = ec2_hostnames(instances)
sshflags = `-o UserKnownHostsFile=/dev/null -o StrictHostKeyChecking=no -o LogLevel=ERROR -i $(ec2_sshkey_file)`
if ec2_install_julia
    println("Installing Julia on launched nodes")
    @sync begin
        for host in hosts
            @async begin
                run(`scp $sshflags jinstall.sh ubuntu@$(host[2]):/home/ubuntu/`)
                run(`ssh $sshflags ubuntu@$(host[2]) "./jinstall.sh"`)
            end
        end
    end

    println("Configuring ssh for passwordless login in-between hosts")
    run(`ssh $sshflags ubuntu@$(hosts[1][2]) "ssh-keygen -b 2048 -t rsa -f /home/ubuntu/.ssh/id_rsa -q -N \"\""`)
    run(`scp $sshflags ubuntu@$(hosts[1][2]):/home/ubuntu/.ssh/id_rsa ubuntu@$(hosts[1][2]):/home/ubuntu/.ssh/id_rsa.pub /tmp/`)

    @sync begin
        for host in hosts
            @async begin
                run(`scp $sshflags /tmp/id_rsa  /tmp/id_rsa.pub ubuntu@$(host[2]):/home/ubuntu/.ssh/`)
                run(`ssh $sshflags ubuntu@$(host[2]) "cat /home/ubuntu/.ssh/id_rsa.pub >> /home/ubuntu/.ssh/authorized_keys"`)
            end
        end
    end

    println("DONE installing Julia on all launched nodes")
end


text:
export ishom, ishet, heterozygosity, het

"""
```
ishom(locus::T) where T <: GenoArray
ishom(locus::Genotype)
ishom(locus::Missing)
```
A series of methods to test if a locus or loci are homozygous and return `true` if
it is, `false` if it isn't, and `missing` if it's `missing`. The vector version
simply broadcasts the function over the elements.
"""
@inline function ishom(locus::Genotype)
    # if allele 1 equels all others, return true
    return all(@inbounds locus[1] .== locus)
end

@inline function ishom(locus::Missing)
    # if the locus is Missing, return missing. no muss no fuss
    return missing
end

@inline function ishom(locus::T) where T<:GenoArray
    return @inbounds ishom.(locus)
end

"""
```
ishet(locus::T) where T <: GenoArray
ishet(locus::Genotype)
ishet(locus::Missing)
```
A series of methods to test if a locus or loci are heterozygous and return `true` if
it is, `false` if it isn't. The vector version simply broadcasts the function over the
elements. Under the hood, this function is simply `!ishom`.
"""
@inline function ishet(locus::Genotype)
    return !ishom(locus)
end

@inline function ishet(locus::Missing)
    # if the locus is Missing, return missing. no muss no fuss
    return missing
end

@inline function ishet(locus::T) where T<:GenoArray
    return @inbounds ishet.(locus)
end

"""
    gene_diversity_nei87(het_exp::Union{Missing,AbstractFloat}, het_obs::Union{Missing,AbstractFloat}, n::Union{Integer, Float64}, corr::Bool = true)
Calculate overall gene diversity with the adjustment/correction given by Nei:

Hₜ = 1 −sum(pbar²ᵢ + Hₛ/(ñ * np) − Het_obs/(2ñ*np))
- _ñ_ is the number of genotypes for a locus for a population
- _np_ is the number of genotypes of a locus across all populations
    - i.e. sum(_ñ_)
- _pbar²_ is the observed homozygosity of a locus for that population
- _Hₛ_ is the within population gene diversity given by:
    - Hₛ =  ñ/(ñ-1) * (1 - sum(pbar²ᵢ - Het_observed / 2ñ))

(Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press).
use `corr = false` to ignore sample-size correction `* n/(n-1)`.
"""
@inline function gene_diversity_nei87(het_exp::AbstractFloat, het_obs::AbstractFloat, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    if corr == true
        corr_val = n/(n-1)
    else
        corr_val = 1
    end
    return @fastmath (het_exp - (het_obs/n/2.0)) * corr_val
end

@inline function gene_diversity_nei87(het_exp::AbstractFloat, het_obs::Missing, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

@inline function gene_diversity_nei87(het_exp::Missing, het_obs::AbstractFloat, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

@inline function gene_diversity_nei87(het_exp::Missing, het_obs::Missing, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

"""
    hetero_o(data::T) where T <: GenoArray
Returns observed heterozygosity as a mean of the number of heterozygous genotypes, defined
as genotypes returning `true` for `ishet()`. This is numerically feasible because
`true` values are mathematically represented as `1`, whereas `false` are represented
as `0`.
"""
@inline function hetero_o(data::T) where T <: GenoArray
    adjusted_vector = ishet(data) |> skipmissing
    isempty(adjusted_vector) ? missing : mean(adjusted_vector)
end

"""
    hetero_e(allele_freqs::Vector{T}) where T <: GenoArray
Returns the expected heterozygosity of an array of genotypes,
calculated as 1 - sum of the squared allele frequencies.
"""
@inline function hetero_e(data::T) where T <: GenoArray
    if all(ismissing.(data)) == true
        return missing
    else
        1.0 - sum(@inbounds allele_freq_vec(data) .^ 2)
    end
end

"""
    heterozygosity(data::PopData, by::String = "locus")
Calculate observed and expected heterozygosity in a `PopData` object. For loci,
heterozygosity is calculated in the Nei fashion, such that heterozygosity is
calculated as the average over heterozygosity per locus per population.
### Modes
- `"locus"` or `"loci"` : heterozygosity per locus (default)
- `"sample"` or `"ind"` or `"individual"` : heterozygosity per individual/sample
- `"population"` or `"pop"` : heterozygosity per population
## Example
heterozygosity(@nancycats, "population" )
"""
function heterozygosity(data::PopData, by::String = "locus")
    if by ∈ ["locus", "loci"]
        tmp = DataFrames.combine(
                groupby(data.loci, [:locus, :population]),
                :genotype => nonmissing => :n_tmp,
                :genotype => hetero_o => :het_pop_obs,
                :genotype => hetero_e => :het_pop_exp
            )
        return DataFrames.combine(
                groupby(tmp, :locus),
                :n_tmp => sum => :n,
                :het_pop_obs => (h_o -> mean(skipmissing(h_o))) => :het_obs,
                :het_pop_exp => (h_e -> mean(skipmissing(h_e))) => :het_exp
            )

    elseif lowercase(by) ∈  ["sample", "ind", "individual"]
        return DataFrames.combine(
                groupby(data.loci, :name),
                :genotype => nonmissing => :n,
                :genotype => hetero_o => :het_obs
            )

    elseif lowercase(by) ∈  ["pop", "population"]
        return DataFrames.combine(
                groupby(data.loci, :population),
                :genotype => nonmissing => :n,
                :genotype => hetero_o => :het_obs,
                :genotype => hetero_e => :het_exp
            )
    else
        error("please specify by = \"locus\", \"sample\", or \"population\"")
    end
end

const het = heterozygosity

"""
    het_sample(data::PopData, individual::String)
Calculate the observed heterozygosity for an individual in a `PopData` object.
"""
@inline function het_sample(data::PopData, individual::String)
    data.loci[data.loci.name .== individual, :genotype] |> hetero_o
end


text:
mg!(st::SolveState;kw...) = iterate!(st,Vcycle!;mxiter=32,kw...)
function mg(A::AbstractMatrix,b::AbstractVector;kw...)
    x=zero(b)
    return x,mg!(mg_state(A,x,b);kw...)
end
mg_state(A,x,b) = fill_children!(SolveState(A,x,residual(A,x,b)))

@inline Vcycle!(st::SolveState;kw...) = Vcycle!(st,st.child;kw...)
@inline Vcycle!(fine::SolveState,coarse::Nothing;kw...) = GS!(fine;kw...)
function Vcycle!(fine::SolveState,coarse::SolveState;
                 precond!::Function=(st;kw...)->GS!(st,inner=0),smooth!::Function=GS!,kw...)
    # set up & solve coarse recursively
    precond!(fine;kw...)
    restrict!(coarse.r,fine.r)
    fill!(coarse.x,0.)
    Vcycle!(coarse;precond!,smooth!,kw...)
    # correct & solve fine
    prolongate!(fine.ϵ,coarse.x)
    increment!(fine)
    smooth!(fine;kw...)
end

restrict!(a,b) = @loop a[I] = sum(@inbounds(b[J]) for J ∈ up(I))
prolongate!(a,b) = @loop a[I] = b[down(I)]

@inline up(I::CartesianIndex{N},a=0) where N = (2I-2oneunit(I)):(2I-oneunit(I)-δ(a,N))
@inline down(I) = CartesianIndex((I+2oneunit(I)).I .÷2)

function fill_children!(st::SolveState)
    if isdivisible(st) && isnothing(st.child)
        st.child = create_child(st.A.L,eltype(st.x))
        fill_children!(st.child)
    end
    return st
end

function create_child(b::AbstractArray{T,N},xT::Type) where {T,N}
    dims = 1 .+ Base.front(size(b)) .÷ 2
    a = zeros(T,dims...,N-1)
    x = FieldVector(zeros(xT,dims))
    for i ∈ 1:N-1, I ∈ x.R
        a[I,i] = 0.5sum(@inbounds(b[J,i]) for J ∈ up(I,i))
    end
    SolveState(Poisson(a),x,zero(x))
end

@inline isdivisible(N::Int) = mod(N,2)==0 && N>2
@inline isdivisible(st::SolveState) = all(size(st.x.R) .|> isdivisible)


text:
###############################################################################
#
# Set titles for the axis
#
###############################################################################
"""
    set_titles!(axs::Array;
                labels::Array{String,1} = String[],
                paren::Bool = true,
                capital::Bool = false,
                dotsep::Bool = false,
                fontsize::Number = 16,
                loc::String = "center",
                usetex::Bool = true
    )

Set titles for the axes, given
- `axs` An array of axis
- `labels` Optional: labels after the panel title, e.g., `(a) label`
- `paren` Optional: if true, use format like `(a)`
- `capital` Optional: if true, use capital letters like `(A)`
- `dotsep` Optional: if true, add a dot after the label like `(a).`
- `fontsize` Optional: fontsize of the title
- `loc` Optional: location of the title
- `usetex` Optional: use latex render
"""
function set_titles!(
            axs::Array;
            labels::Array{String,1} = String[],
            paren::Bool = true,
            capital::Bool = false,
            dotsep::Bool = false,
            fontsize::Number = 16,
            loc::String = "center",
            usetex::Bool = true
)
    letters = ["a", "b", "c", "d", "e", "f", "g",
               "h", "i", "j", "k", "l", "m", "n",
               "o", "p", "q", "r", "s", "t",
               "u", "v", "w", "x", "y", "z"];

    # set titles for each panel
    for indx in eachindex(axs)
        _title  = "";
        _letter = letters[indx];
        _LETTER = uppercase(letters[indx]);

        usetex  ? _title *= "\\textbf{" : nothing;
        paren   ? _title *= "("         : nothing;
        capital ? _title *= _LETTER     : _title *= _letter;
        paren   ? _title *= ")"         : nothing;
        dotsep  ? _title *= "."         : nothing;
        usetex  ? _title *= "}"         : nothing;

        length(labels)>0 ? _title *= " " * labels[indx] : nothing;

        if usetex
            axs[indx].set_title(_title, fontsize=fontsize, loc=loc);
        else
            axs[indx].set_title(_title, fontsize=fontsize, loc=loc,
                                     fontweight="bold");
        end
    end

    return nothing
end


text:
@testset "index" begin
@testset "index_opens" begin
    @test Bgen(BGEN.datadir("example.15bits.bgen")).idx === nothing
    @test Bgen(BGEN.datadir("example.16bits.bgen")).idx !== nothing
end
@testset "index_region" begin
    chrom = "01"
    start = 5000
    stop = 50000

    b = Bgen(BGEN.datadir("example.16bits.bgen"))
    idx = b.idx
    @test length(select_region(b, chrom)) == length(gen_data)
    @test length(select_region(b, "02")) == 0
    @test length(select_region(b, chrom;
        start=start * 100, stop=stop * 100)) == 0

    after_pos_offsets = select_region(b, chrom; start=start)
    @test length(after_pos_offsets) ==
        length(filter(x -> x.pos >= start, gen_data))

    in_region_offsets = select_region(b, chrom; start=start, stop=stop)
    @test length(in_region_offsets) == length(filter(x -> start <= x.pos <=
        stop, gen_data))
    @test rsid(variant_by_rsid(b, "RSID_10")) == "RSID_10"
    @test rsid(variant_by_index(b, 4)) == "RSID_3"
end

end


text:
function UCST_mix(model::EoSModel,T;v0=nothing)
    if v0 == nothing
        v0 = x0_UCST_mix(model,T)
    end
    f! = (F,x) -> Obj_UCST_mix(model, F, x[2], exp10(x[1]), T)
    r  = Solvers.nlsolve(f!,v0[1:end],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    z_c = FractionVector(sol[2])
    V_c = exp10(sol[1])
    p_c = pressure(model, V_c, T, z_c)
    return (p_c, V_c, z_c)
end
function x0_UCST_mix(model::EoSModel,T)
    V  = x0_volume_liquid(model,T,[0.5,0.5])
    return [log10(V),0.5]
end

function Obj_UCST_mix(model::EoSModel,F,z,V,T)
    z    = FractionVector(z)
    f(x) = eos(model,V,T,x)
    H(x) = ForwardDiff.hessian(f,x)/8.134/T
    L(x) = det(H(x))
    dL(x) = ForwardDiff.gradient(L,x)
    M(x) = [H(x)[1:end-1,:];transpose(dL(x))]
    F[1] = L(z)
    F[2] = det(M(z))
    return F
end


text:
###############################################################################
#
# Parse h and v index from MODIS file name
#
###############################################################################
"""
    parse_HV(modis_name::String)

Return the H and V index, given
- `modis_name` MODIS file name
"""
function parse_HV(modis_name::String)
    i = findfirst(".h", modis_name)[1];
    h = parse(Int, modis_name[i+2:i+3]) + 1;
    v = parse(Int, modis_name[i+5:i+6]) + 1;

    return h,v
end


text:
# This file is adapted from
const CACHE_FLAG = Ref{Bool}(true)

"""
    enable_cache()

Enable command compile cache. See also [`disable_cache`](@ref).
"""
function enable_cache()
    CACHE_FLAG[] = true
    return
end

"""
    disable_cache()

Disable command compile cache. See also [`enable_cache`](@ref).
"""
function disable_cache()
    CACHE_FLAG[] = false
    return
end

"""
    @cast <function expr>
    @cast <module expr>
    @cast <function name>
    @cast <module name>

Cast a Julia object to a command object. `@cast` will
always execute the given expression, and only create and
register an command object via [`command`](@ref) after
analysing the given expression.

# Example

The most basic way is to use `@cast` on a function as your command
entry.

```julia
@cast function your_command(arg1, arg2::Int)
    # your processing code
end
```

This will create a [`LeafCommand`](@ref) object and register it
to the current module's `CASTED_COMMANDS` constant. You can access
this object via `MODULE_NAME.CASTED_COMMANDS["your_command"]`.

Note this will not create a functional CLI, to create a function
CLI you need to create an entry point, which can be declared by
[`@main`](@ref).
"""
macro cast(ex)
    if CACHE_FLAG[] && iscached()
        return esc(ex)
    else
        return esc(cast_m(__module__, QuoteNode(__source__), ex))
    end
end

function cast_m(m::Module, line::QuoteNode, ex)
    ret = Expr(:block)
    pushmaybe!(ret, create_casted_commands(m))

    if ex isa Symbol
        push!(ret.args, parse_module(m, line, ex))
        return ret
    end

    def = splitdef(ex; throw = false)
    if def === nothing
        push!(ret.args, parse_module(m, line, ex))
        return ret
    end

    push!(ret.args, parse_function(m, line, ex, def))
    return ret
end

# Entry
"""
    @main <function expr>
    @main [options...]

Create an `EntryCommand` and use it as the entry of the entire CLI.
If you only have one function to cast to a CLI, you can use `@main`
instead of `@cast` so it will create both the command object and
the entry.

If you have declared commands via `@cast`, you can create the entry
via `@main [options...]`, available options are:

- `name`: default is the current module name in lowercase.
- `version`: default is the current project version or `v"0.0.0"`. If
    it's `v"0.0.0"`, the version will not be printed.
- `doc`: a description of the entry command.
"""
macro main(xs...)
    return esc(main_m(__module__, QuoteNode(__source__), xs...))
end

function main_m(m::Module, line::QuoteNode, ex::Expr)
    if CACHE_FLAG[] && iscached()
        return quote
            Core.@__doc__ $ex
            include($(cachefile()[1]))
        end
    end

    ret = Expr(:block)
    def = splitdef(ex; throw = false)
    var_cmd, var_entry = gensym(:cmd), gensym(:entry)
    push!(ret.args, ex)

    if def === nothing
        ex.head === :(=) && return create_entry(m, line, ex)
        ex.head === :module ||
            throw(Meta.ParseError("invalid expression, can only cast functions or modules"))
        cmd = xcall(command, ex.args[2]; line = line)
        push!(ret.args, :($var_cmd = $cmd))
    else
        push!(ret.args, :(Core.@__doc__ $(def[:name])))
        cmd = xcall(command, def[:name], parse_args(def), parse_kwargs(def); line = line)
        push!(ret.args, :($var_cmd = $cmd))
    end

    push!(ret.args, :($var_entry = $(xcall(Types, :EntryCommand, var_cmd; version = get_version(m)))))
    push!(ret.args, precompile_or_exec(m, var_entry))
    return ret
end

function main_m(m::Module, line::QuoteNode, ex::Symbol)
    CACHE_FLAG[] && iscached() && return :(include($(cachefile()[1])))
    var_cmd, var_entry = gensym(:cmd), gensym(:entry)
    quote
        $var_cmd = $(xcall(command, ex; line = line))
        $var_entry = $(xcall(Types, :EntryCommand, var_cmd; line = line))
        $(precompile_or_exec(m, var_entry))
    end
end

function main_m(m::Module, line::QuoteNode, kwargs...)
    CACHE_FLAG[] && iscached() && return :(include($(cachefile()[1])))
    return create_entry(m, line, kwargs...)
end

function create_entry(m::Module, line::QuoteNode, kwargs...)
    configs = Dict{Symbol,Any}(:name => default_name(m), :version => get_version(m), :doc => "")
    for kw in kwargs
        for key in [:name, :version, :doc]
            if kw.args[1] === key
                configs[key] = kw.args[2]
            end
        end
    end

    ret = Expr(:block)

    var_cmd, var_entry = gensym(:cmd), gensym(:entry)
    cmd = xcall(
        Types,
        :NodeCommand,
        configs[:name],
        line,
        configs[:doc],
        :(collect(values($m.CASTED_COMMANDS))),
    )
    entry = xcall(Types, :EntryCommand, var_cmd, configs[:version], line)

    push!(ret.args, :($var_cmd = $cmd))
    push!(ret.args, :($var_entry = $entry))
    push!(ret.args, precompile_or_exec(m, var_entry))
    return ret
end

function precompile_or_exec(m::Module, entry)
    if m == Main && CACHE_FLAG[]
        return quote
            $create_cache($entry)
            include($(cachefile()[1]))
        end
    elseif m == Main
        return quote
            $(xcall(m, :eval, xcall(CodeGen, :codegen, entry)))
            poptart_main()
        end
    else
        quote
            $create_cache($entry)

            $(create_casted_commands(m))
            $(xcall(set_cmd!, casted_commands(m), entry, "main"))
            $(xcall(m, :eval, xcall(CodeGen, :codegen, entry)))

            precompile(Tuple{typeof($m.poptart_main),Array{String,1}})
        end
    end
end