text:
"""
    MPO

A finite size matrix product operator type.
Keeps track of the orthogonality center.
"""
mutable struct MPO <: AbstractMPS
  data::Vector{ITensor}
  llim::Int
  rlim::Int
end

function MPO(A::Vector{<:ITensor}; ortho_lims::UnitRange=1:length(A))
  return MPO(A, first(ortho_lims) - 1, last(ortho_lims) + 1)
end

set_data(A::MPO, data::Vector{ITensor}) = MPO(data, A.llim, A.rlim)

MPO() = MPO(ITensor[], 0, 0)

function convert(::Type{MPS}, M::MPO)
  return MPS(data(M); ortho_lims=ortho_lims(M))
end

function convert(::Type{MPO}, M::MPS)
  return MPO(data(M); ortho_lims=ortho_lims(M))
end

function MPO(::Type{ElT}, sites::Vector{<:Index}) where {ElT<:Number}
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  if N == 0
    return MPO()
  elseif N == 1
    v[1] = emptyITensor(ElT, dag(sites[1]), sites[1]')
    return MPO(v)
  end
  space_ii = all(hasqns, sites) ? [QN() => 1] : 1
  l = [Index(space_ii, "Link,l=$ii") for ii in 1:(N - 1)]
  for ii in eachindex(sites)
    s = sites[ii]
    if ii == 1
      v[ii] = emptyITensor(ElT, dag(s), s', l[ii])
    elseif ii == N
      v[ii] = emptyITensor(ElT, dag(l[ii - 1]), dag(s), s')
    else
      v[ii] = emptyITensor(ElT, dag(l[ii - 1]), dag(s), s', l[ii])
    end
  end
  return MPO(v)
end

MPO(sites::Vector{<:Index}) = MPO(Float64, sites)

"""
    MPO(N::Int)

Make an MPO of length `N` filled with default ITensors.
"""
MPO(N::Int) = MPO(Vector{ITensor}(undef, N))

"""
    MPO([::Type{ElT} = Float64}, ]sites, ops::Vector{String})

Make an MPO with pairs of sites `s[i]` and `s[i]'`
and operators `ops` on each site.
"""
function MPO(::Type{ElT}, sites::Vector{<:Index}, ops::Vector{String}) where {ElT<:Number}
  N = length(sites)
  ampo = OpSum() + [ops[n] => n for n in 1:N]
  M = MPO(ampo, sites)

  # Currently, OpSum does not output the optimally truncated
  # MPO (see https://github.com/ITensor/ITensors.jl/issues/526)
  # So here, we need to first normalize, then truncate, then
  # return the normalization.
  lognormM = lognorm(M)
  M ./= exp(lognormM / N)
  truncate!(M; cutoff=1e-15)
  M .*= exp(lognormM / N)
  return M
end

function MPO(::Type{ElT}, sites::Vector{<:Index}, fops::Function) where {ElT<:Number}
  ops = [fops(n) for n in 1:length(sites)]
  return MPO(ElT, sites, ops)
end

MPO(sites::Vector{<:Index}, ops) = MPO(Float64, sites, ops)

"""
    MPO([::Type{ElT} = Float64, ]sites, op::String)

Make an MPO with pairs of sites `s[i]` and `s[i]'`
and operator `op` on every site.
"""
function MPO(::Type{ElT}, sites::Vector{<:Index}, op::String) where {ElT<:Number}
  return MPO(ElT, sites, fill(op, length(sites)))
end

MPO(sites::Vector{<:Index}, op::String) = MPO(Float64, sites, op)

function randomMPO(sites::Vector{<:Index}, m::Int=1)
  M = MPO(sites, "Id")
  for i in eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  m > 1 && throw(ArgumentError("randomMPO: currently only m==1 supported"))
  return M
end

function MPO(A::ITensor, sites::Vector{<:Index}; kwargs...)
  return MPO(A, IndexSet.(prime.(sites), dag.(sites)); kwargs...)
end

function outer_mps_mps_deprecation_warning()
  return "Calling `outer(ψ::MPS, ϕ::MPS)` for MPS `ψ` and `ϕ` with shared indices is deprecated. Currently, we automatically prime `ψ` to make sure the site indices don't clash, but that will no longer be the case in ITensors v0.4. To upgrade your code, call `outer(ψ', ϕ)`. Although the new interface seems less convenient, it will allow `outer` to accept more general outer products going forward, such as outer products where some indices are shared (a batched outer product) or outer products of MPS between site indices that aren't just related by a single prime level."
end

function deprecate_make_inds_unmatch(::typeof(outer), ψ::MPS, ϕ::MPS; kw...)
  if hassameinds(siteinds, ψ, ϕ)
    warn_once(outer_mps_mps_deprecation_warning(), :outer_mps_mps)
    ψ = ψ'
  end
  return ψ, ϕ
end

"""
    outer(x::MPS, y::MPS; <keyword argument>) -> MPO

Compute the outer product of `MPS` `x` and `MPS` `y`,
returning an `MPO` approximation. Note that `y` will be conjugated.

In Dirac notation, this is the operation `|x⟩⟨y|`.

If you want an outer product of an MPS with itself, you should
call `outer(x', x; kwargs...)` so that the resulting MPO
has site indices with indices coming in pairs of prime levels
of 1 and 0. If not, the site indices won't be unique which would
not be an outer product.

For example:
```julia
s = siteinds("S=1/2", 5)
x = randomMPS(s)
y = randomMPS(s)
outer(x, y) # Incorrect! Site indices must be unique.
outer(x', y) # Results in an MPO with pairs of primed and unprimed indices.
```
This allows for more general outer products, such as more general
MPO outputs which don't have pairs of primed and unprimed indices,
or outer products where the input MPS are vectorizations of MPOs.

For example:
```julia
s = siteinds("S=1/2", 5)
X = MPO(s, "Id")
Y = MPO(s, "Id")
x = convert(MPS, X)
y = convert(MPS, Y)
outer(x, y) # Incorrect! Site indices must be unique.
outer(x', y) # Incorrect! Site indices must be unique.
outer(addtags(x, "Out"), addtags(y, "In")) # This performs a proper outer product.
```

The keyword arguments determine the truncation, and accept
the same arguments as `contract(::MPO, ::MPO; kwargs...)`.

See also [`apply`](@ref), [`contract`](@ref).
"""
function outer(ψ::MPS, ϕ::MPS; kw...)
  ψ, ϕ = deprecate_make_inds_unmatch(outer, ψ, ϕ; kw...)

  ψmat = convert(MPO, ψ)
  ϕmat = convert(MPO, dag(ϕ))
  return contract(ψmat, ϕmat; kw...)
end

"""
    projector(x::MPS; <keyword argument>) -> MPO

Computes the projector onto the state `x`. In Dirac notation, this is the operation `|x⟩⟨x|/|⟨x|x⟩|²`.

Use keyword arguments to control the level of truncation, which are
the same as those accepted by `contract(::MPO, ::MPO; kw...)`.

# Keywords
- `normalize::Bool=true`: whether or not to normalize the input MPS before forming the projector. If `normalize==false` and the input MPS is not already normalized, this function will not output a proper project, and simply outputs `outer(x, x) = |x⟩⟨x|`, i.e. the projector scaled by `norm(x)^2`.
- truncation keyword arguments accepted by `contract(::MPO, ::MPO; kw...)`.

See also [`outer`](@ref), [`contract`](@ref).
"""
function projector(ψ::MPS; normalize::Bool=true, kw...)
  ψψᴴ = outer(ψ', ψ; kw...)
  if normalize
    normalize!(ψψᴴ[orthocenter(ψψᴴ)])
  end
  return ψψᴴ
end

# XXX: rename originalsiteind?
"""
    siteind(M::MPO, j::Int; plev = 0, kwargs...)

Get the first site Index of the MPO found, by
default with prime level 0.
"""
siteind(M::MPO, j::Int; kwargs...) = siteind(first, M, j; plev=0, kwargs...)

# TODO: make this return the site indices that would have
# been used to create the MPO? I.e.:
# [dag(siteinds(M, j; plev = 0, kwargs...)) for j in 1:length(M)]
"""
    siteinds(M::MPO; kwargs...)

Get a Vector of IndexSets of all the site indices of M.
"""
siteinds(M::MPO; kwargs...) = siteinds(all, M; kwargs...)

function siteinds(Mψ::Tuple{MPO,MPS}, n::Int; kwargs...)
  return siteinds(uniqueinds, Mψ[1], Mψ[2], n; kwargs...)
end

function nsites(Mψ::Tuple{MPO,MPS})
  M, ψ = Mψ
  N = length(M)
  @assert N == length(ψ)
  return N
end

siteinds(Mψ::Tuple{MPO,MPS}; kwargs...) = [siteinds(Mψ, n; kwargs...) for n in 1:nsites(Mψ)]

# XXX: rename originalsiteinds?
"""
    firstsiteinds(M::MPO; kwargs...)

Get a Vector of the first site Index found on each site of M.

By default, it finds the first site Index with prime level 0.
"""
firstsiteinds(M::MPO; kwargs...) = siteinds(first, M; plev=0, kwargs...)

function hassameinds(::typeof(siteinds), ψ::MPS, Hϕ::Tuple{MPO,MPS})
  N = length(ψ)
  @assert N == length(Hϕ[1]) == length(Hϕ[1])
  for n in 1:N
    !hassameinds(siteinds(Hϕ, n), siteinds(ψ, n)) && return false
  end
  return true
end

function inner_mps_mpo_mps_deprecation_warning()
  return """
 Calling `inner(x::MPS, A::MPO, y::MPS)` where the site indices of the `MPS` `x` and the `MPS` resulting from contracting `MPO` `A` with `MPS` `y` don't match is deprecated as of ITensors v0.3 and will result in an error in ITensors v0.4. The most common cause of this is something like the following:
 ```julia
 s = siteinds("S=1/2")
 psi = randomMPS(s)
 H = MPO(s, "Id")
 inner(psi, H, psi)
 ```
 `psi` has the Index structure `-s-(psi)` and `H` has the Index structure `-s'-(H)-s-`, so the Index structure of would be `(dag(psi)-s- -s'-(H)-s-(psi)` unless the prime levels were fixed. Previously we tried fixing the prime level in situations like this, but we will no longer be doing that going forward.

 There are a few ways to fix this. You can simply change:
 ```julia
 inner(psi, H, psi)
 ```
 to:
 ```julia
 inner(psi', H, psi)
 ```
 in which case the Index structure will be `(dag(psi)-s'-(H)-s-(psi)`.

 Alternatively, you can use the `Apply` function:
 ```julia
 inner(psi, Apply(H, psi))
 ```
 In this case, `Apply(H, psi)` represents the "lazy" evaluation of `apply(H, psi)`. The function `apply(H, psi)` performs the contraction of `H` with `psi` and then unprimes the results, so this versions ensures that the prime levels of the inner product will match.

 Although the new behavior seems less convenient, it makes it easier to generalize `inner(::MPS, ::MPO, ::MPS)` to other types of inputs, like `MPS` and `MPO` with different tag and prime conventions, multiple sites per tensor, `ITensor` inputs, etc.
 """
end

function deprecate_make_inds_match!(
  ::typeof(dot), ydag::MPS, A::MPO, x::MPS; make_inds_match::Bool=true
)
  N = length(x)
  if !hassameinds(siteinds, ydag, (A, x))
    sAx = siteinds((A, x))
    if any(s -> length(s) > 1, sAx)
      n = findfirst(n -> !hassameinds(siteinds(ydag, n), siteinds((A, x), n)), 1:N)
      error(
        """Calling `dot(ϕ::MPS, H::MPO, ψ::MPS)` with multiple site indices per MPO/MPS tensor but the site indices don't match. Even with `make_inds_match = true`, the case of multiple site indices per MPO/MPS is not handled automatically. The sites with unmatched site indices are:

            inds(ϕ[$n]) = $(inds(ydag[n]))

            inds(H[$n]) = $(inds(A[n]))

            inds(ψ[$n]) = $(inds(x[n]))

        Make sure the site indices of your MPO/MPS match. You may need to prime one of the MPS, such as `dot(ϕ', H, ψ)`.""",
      )
    end
    if !hassameinds(siteinds, ydag, (A, x)) && make_inds_match
      warn_once(inner_mps_mpo_mps_deprecation_warning(), :inner_mps_mpo_mps)
      replace_siteinds!(ydag, sAx)
    end
  end
  return ydag, A, x
end

"""
    dot(y::MPS, A::MPO, x::MPS)

Same as [`inner`](@ref).
"""
function dot(y::MPS, A::MPO, x::MPS; make_inds_match::Bool=true, kwargs...)::Number
  N = length(A)
  check_hascommoninds(siteinds, A, x)
  ydag = dag(y)
  sim!(linkinds, ydag)
  ydag, A, x = deprecate_make_inds_match!(dot, ydag, A, x; make_inds_match)
  check_hascommoninds(siteinds, A, y)
  O = ydag[1] * A[1] * x[1]
  for j in 2:N
    O = O * ydag[j] * A[j] * x[j]
  end
  return O[]
end

"""
    inner(y::MPS, A::MPO, x::MPS)

Compute `⟨y|A|x⟩ = ⟨y|Ax⟩` efficiently and exactly without making any intermediate
MPOs. In general it is more efficient and accurate than `inner(y, apply(A, x))`.

This is helpful for computing the expectation value of an operator `A`, which would be:
```julia
inner(x, A, x)
```
assuming `x` is normalized.

If you want to compute `⟨By|Ax⟩` you can use `inner(B::MPO, y::MPS, A::MPO, x::MPS)`.

This is helpful for computing the variance of an operator `A`, which would be:
```julia
inner(A, x, A, x) - inner(x, A, x) ^ 2
```
assuming `x` is normalized.

$(make_inds_match_docstring_warning())

Same as [`dot`](@ref).
"""
inner(y::MPS, A::MPO, x::MPS; kwargs...) = dot(y, A, x; kwargs...)

function inner(y::MPS, Ax::Apply{Tuple{MPO,MPS}})
  return inner(y', Ax.args[1], Ax.args[2])
end

"""
    dot(B::MPO, y::MPS, A::MPO, x::MPS)

Same as [`inner`](@ref).
"""
function dot(B::MPO, y::MPS, A::MPO, x::MPS; make_inds_match::Bool=true, kwargs...)::Number
  !make_inds_match && error(
    "make_inds_match = false not currently supported in dot(::MPO, ::MPS, ::MPO, ::MPS)"
  )
  N = length(B)
  if length(y) != N || length(x) != N || length(A) != N
    throw(
      DimensionMismatch(
        "inner: mismatched lengths $N and $(length(x)) or $(length(y)) or $(length(A))"
      ),
    )
  end
  ydag = dag(y)
  prime!(ydag, 2)
  Bdag = dag(B)
  prime!(Bdag)
  # Swap prime levels 1 -> 2 and 2 -> 1.
  for j in eachindex(Bdag)
    Axcommon = commonind(A[j], x[j])
    ABcommon = uniqueind(filterinds(A[j]; tags="Site"), IndexSet(Axcommon))
    swapprime!(Bdag[j], 2, 3)
    swapprime!(Bdag[j], 1, 2)
    swapprime!(Bdag[j], 3, 1)
    noprime!(Bdag[j], prime(ABcommon, 2))
  end
  yB = ydag[1] * Bdag[1]
  Ax = A[1] * x[1]
  O = yB * Ax
  for j in 2:N
    yB = ydag[j] * Bdag[j]
    Ax = A[j] * x[j]
    yB *= O
    O = yB * Ax
  end
  return O[]
end

# TODO: maybe make these into tuple inputs?
# Also can generalize to:
# inner((β, B, y), (α, A, x))
"""
    inner(B::MPO, y::MPS, A::MPO, x::MPS)

Compute `⟨By|A|x⟩ = ⟨By|Ax⟩` efficiently and exactly without making any intermediate
MPOs. In general it is more efficient and accurate than `inner(apply(B, y), apply(A, x))`.

This is helpful for computing the variance of an operator `A`, which would be:
```julia
inner(A, x, A, x) - inner(x, A, x) ^ 2
```

$(make_inds_match_docstring_warning())

Same as [`dot`](@ref).
"""
inner(B::MPO, y::MPS, A::MPO, x::MPS) = dot(B, y, A, x)

function dot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...)
  if make_inds_match
    error("In dot(::MPO, ::MPO), make_inds_match is not currently supported")
  end
  return _log_or_not_dot(M1, M2, false; make_inds_match=make_inds_match)
end

# TODO: implement by combining the MPO indices and converting
# to MPS
function logdot(M1::MPO, M2::MPO; make_inds_match::Bool=false, kwargs...)
  if make_inds_match
    error("In dot(::MPO, ::MPO), make_inds_match is not currently supported")
  end
  return _log_or_not_dot(M1, M2, true; make_inds_match=make_inds_match)
end

function tr(M::MPO; plev::Pair{Int,Int}=0 => 1, tags::Pair=ts"" => ts"")
  N = length(M)
  #
  # TODO: choose whether to contract or trace
  # first depending on the bond dimension. The scaling is:
  #
  # 1. Trace last:  O(χ²d²) + O(χd²)
  # 2. Trace first: O(χ²d²) + O(χ²)
  #
  # So tracing first is better if d > √χ.
  #
  L = tr(M[1]; plev=plev, tags=tags)
  for j in 2:N
    L *= M[j]
    L = tr(L; plev=plev, tags=tags)
  end
  return L
end

"""
    error_contract(y::MPS, A::MPO, x::MPS;
                   make_inds_match::Bool = true)
    error_contract(y::MPS, x::MPS, x::MPO;
                   make_inds_match::Bool = true)

Compute the distance between A|x> and an approximation MPS y:
`| |y> - A|x> |/| A|x> | = √(1 + (<y|y> - 2*real(<y|A|x>))/<Ax|A|x>)`.

If `make_inds_match = true`, the function attempts match the site
indices of `y` with the site indices of `A` that are not common
with `x`.
"""
function error_contract(y::MPS, A::MPO, x::MPS; kwargs...)
  N = length(A)
  if length(y) != N || length(x) != N
    throw(
      DimensionMismatch("inner: mismatched lengths $N and $(length(x)) or $(length(y))")
    )
  end
  iyy = dot(y, y; kwargs...)
  iyax = dot(y', A, x; kwargs...)
  iaxax = dot(A, x, A, x; kwargs...)
  return sqrt(abs(1.0 + (iyy - 2 * real(iyax)) / iaxax))
end

error_contract(y::MPS, x::MPS, A::MPO) = error_contract(y, A, x)

"""
    apply(A::MPO, x::MPS; kwargs...)

Contract the `MPO` `A` with the `MPS` `x` and then map the prime level of the resulting
MPS back to 0.

Equivalent to `replaceprime(contract(A, x; kwargs...), 2 => 1)`.

See also [`contract`](@ref) for details about the arguments available.
"""
function apply(A::MPO, ψ::MPS; kwargs...)
  Aψ = contract(A, ψ; kwargs...)
  return replaceprime(Aψ, 1 => 0)
end

(A::MPO)(ψ::MPS; kwargs...) = apply(A, ψ; kwargs...)

function contract(A::MPO, ψ::MPS; alg="densitymatrix", kwargs...)
  if haskey(kwargs, :method)
    # Backwards compatibility, use `method`.
    alg = get(kwargs, :method, "densitymatrix")
  end

  # Keyword argument deprecations
  if alg == "DensityMatrix"
    @warn "In contract, method DensityMatrix is deprecated in favor of densitymatrix"
    alg = "densitymatrix"
  end
  if alg == "Naive"
    @warn "In contract, `alg=\"Naive\"` is deprecated in favor of `alg=\"naive\"`"
    alg = "naive"
  end

  return contract(Algorithm(alg), A, ψ; kwargs...)
end

contract_mpo_mps_doc = """
    contract(ψ::MPS, A::MPO; kwargs...) -> MPS
    *(::MPS, ::MPO; kwargs...) -> MPS

    contract(A::MPO, ψ::MPS; kwargs...) -> MPS
    *(::MPO, ::MPS; kwargs...) -> MPS

Contract the `MPO` `A` with the `MPS` `ψ`, returning an `MPS` with the unique
site indices of the `MPO`.

For example, for an MPO with site indices with prime levels of 1 and 0, such as
`-s'-A-s-`, and an MPS with site indices with prime levels of 0, such as
`-s-x`, the result is an MPS `y` with site indices with prime levels of 1,
`-s'-y = -s'-A-s-x`.

Since it is common to contract an MPO with prime levels of 1 and 0 with an MPS with
prime level of 0 and want a resulting MPS with prime levels of 0, we provide a
convenience function `apply`:
```julia
apply(A, x; kwargs...) = replaceprime(contract(A, x; kwargs...), 2 => 1)`.
```

Choose the method with the `method` keyword, for example
`"densitymatrix"` and `"naive"`.

# Keywords
- `cutoff::Float64=1e-13`: the cutoff value for truncating the density matrix eigenvalues. Note that the default is somewhat arbitrary and subject to change, in general you should set a `cutoff` value.
- `maxdim::Int=maxlinkdim(A) * maxlinkdim(ψ))`: the maximal bond dimension of the results MPS.
- `mindim::Int=1`: the minimal bond dimension of the resulting MPS.
- `normalize::Bool=false`: whether or not to normalize the resulting MPS.
- `method::String="densitymatrix"`: the algorithm to use for the contraction. Currently the options are "densitymatrix", where the network formed by the MPO and MPS is squared and contracted down to a density matrix which is diagonalized iteratively at each site, and "naive", where the MPO and MPS tensor are contracted exactly at each site and then a truncation of the resulting MPS is performed.

See also [`apply`](@ref).
"""

@doc """
$contract_mpo_mps_doc
""" contract(::MPO, ::MPS)

contract(ψ::MPS, A::MPO; kwargs...) = contract(A, ψ; kwargs...)

*(A::MPO, B::MPS; kwargs...) = contract(A, B; kwargs...)
*(A::MPS, B::MPO; kwargs...) = contract(A, B; kwargs...)

# TODO: try this to copy the docstring
# Causing an error in Revise
#@doc """
#$contract_mpo_mps_doc
#""" *(::MPO, ::MPS)

#@doc (@doc contract(::MPO, ::MPS)) *(::MPO, ::MPS)

function contract(::Algorithm"densitymatrix", A::MPO, ψ::MPS; kwargs...)::MPS
  n = length(A)
  n != length(ψ) &&
    throw(DimensionMismatch("lengths of MPO ($n) and MPS ($(length(ψ))) do not match"))
  if n == 1
    return MPS([A[1] * ψ[1]])
  end

  ψ_out = similar(ψ)
  cutoff::Float64 = get(kwargs, :cutoff, 1e-13)
  requested_maxdim::Int = get(kwargs, :maxdim, maxlinkdim(A) * maxlinkdim(ψ))
  mindim::Int = max(get(kwargs, :mindim, 1), 1)
  normalize::Bool = get(kwargs, :normalize, false)

  any(i -> isempty(i), siteinds(commoninds, A, ψ)) &&
    error("In `contract(A::MPO, x::MPS)`, `A` and `x` must share a set of site indices")

  # In case A and ψ have the same link indices
  A = sim(linkinds, A)

  ψ_c = dag(ψ)
  A_c = dag(A)

  # To not clash with the link indices of A and ψ
  sim!(linkinds, A_c)
  sim!(linkinds, ψ_c)
  sim!(siteinds, commoninds, A_c, ψ_c)

  # A version helpful for making the density matrix
  simA_c = sim(siteinds, uniqueinds, A_c, ψ_c)

  # Store the left environment tensors
  E = Vector{ITensor}(undef, n - 1)

  E[1] = ψ[1] * A[1] * A_c[1] * ψ_c[1]
  for j in 2:(n - 1)
    E[j] = E[j - 1] * ψ[j] * A[j] * A_c[j] * ψ_c[j]
  end
  R = ψ[n] * A[n]
  simR_c = ψ_c[n] * simA_c[n]
  ρ = E[n - 1] * R * simR_c
  l = linkind(ψ, n - 1)
  ts = isnothing(l) ? "" : tags(l)
  Lis = siteinds(uniqueinds, A, ψ, n)
  Ris = siteinds(uniqueinds, simA_c, ψ_c, n)
  F = eigen(ρ, Lis, Ris; ishermitian=true, tags=ts, kwargs...)
  D, U, Ut = F.D, F.V, F.Vt
  l_renorm, r_renorm = F.l, F.r
  ψ_out[n] = Ut
  R = R * dag(Ut) * ψ[n - 1] * A[n - 1]
  simR_c = simR_c * U * ψ_c[n - 1] * simA_c[n - 1]
  for j in reverse(2:(n - 1))
    # Determine smallest maxdim to use
    cip = commoninds(ψ[j], E[j - 1])
    ciA = commoninds(A[j], E[j - 1])
    prod_dims = dim(cip) * dim(ciA)
    maxdim = min(prod_dims, requested_maxdim)

    s = siteinds(uniqueinds, A, ψ, j)
    s̃ = siteinds(uniqueinds, simA_c, ψ_c, j)
    ρ = E[j - 1] * R * simR_c
    l = linkind(ψ, j - 1)
    ts = isnothing(l) ? "" : tags(l)
    Lis = IndexSet(s..., l_renorm)
    Ris = IndexSet(s̃..., r_renorm)
    F = eigen(ρ, Lis, Ris; ishermitian=true, maxdim=maxdim, tags=ts, kwargs...)
    D, U, Ut = F.D, F.V, F.Vt
    l_renorm, r_renorm = F.l, F.r
    ψ_out[j] = Ut
    R = R * dag(Ut) * ψ[j - 1] * A[j - 1]
    simR_c = simR_c * U * ψ_c[j - 1] * simA_c[j - 1]
  end
  if normalize
    R ./= norm(R)
  end
  ψ_out[1] = R
  setleftlim!(ψ_out, 0)
  setrightlim!(ψ_out, 2)
  return ψ_out
end

function contract(::Algorithm"naive", A::MPO, ψ::MPS; kwargs...)::MPS
  truncate = get(kwargs, :truncate, true)

  N = length(A)
  if N != length(ψ)
    throw(DimensionMismatch("lengths of MPO ($N) and MPS ($(length(ψ))) do not match"))
  end

  ψ_out = MPS(N)
  for j in 1:N
    ψ_out[j] = A[j] * ψ[j]
  end

  for b in 1:(N - 1)
    Al = commonind(A[b], A[b + 1])
    pl = commonind(ψ[b], ψ[b + 1])
    C = combiner(Al, pl)
    ψ_out[b] *= C
    ψ_out[b + 1] *= dag(C)
  end

  if truncate
    truncate!(ψ_out; kwargs...)
  end

  return ψ_out
end

function contract(A::MPO, B::MPO; kwargs...)
  if hassameinds(siteinds, A, B)
    error(
      "In `contract(A::MPO, B::MPO)`, MPOs A and B have the same site indices. The indices of the MPOs in the contraction are taken literally, and therefore they should only share on site index per site so the contraction results in an MPO. You may want to use `replaceprime(contract(A', B), 2 => 1)` or `apply(A, B)` which automatically adjusts the prime levels assuming the input MPOs have pairs of primed and unprimed indices.",
    )
  end
  cutoff::Float64 = get(kwargs, :cutoff, 1e-14)
  resp_degen::Bool = get(kwargs, :respect_degenerate, true)
  maxdim::Int = get(kwargs, :maxdim, maxlinkdim(A) * maxlinkdim(B))
  mindim::Int = max(get(kwargs, :mindim, 1), 1)
  N = length(A)
  N != length(B) &&
    throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
  # Special case for a single site
  N == 1 && return MPO([A[1] * B[1]])
  A = orthogonalize(A, 1)
  B = orthogonalize(B, 1)
  A = sim(linkinds, A)
  sA = siteinds(uniqueinds, A, B)
  sB = siteinds(uniqueinds, B, A)
  C = MPO(N)
  lCᵢ = Index[]
  R = ITensor(1)
  for i in 1:(N - 2)
    RABᵢ = R * A[i] * B[i]
    left_inds = [sA[i]..., sB[i]..., lCᵢ...]
    C[i], R = factorize(
      RABᵢ,
      left_inds;
      ortho="left",
      tags=commontags(linkinds(A, i)),
      cutoff=cutoff,
      maxdim=maxdim,
      mindim=mindim,
      kwargs...,
    )
    lCᵢ = dag(commoninds(C[i], R))
  end
  i = N - 1
  RABᵢ = R * A[i] * B[i] * A[i + 1] * B[i + 1]
  left_inds = [sA[i]..., sB[i]..., lCᵢ...]
  C[N - 1], C[N] = factorize(
    RABᵢ,
    left_inds;
    ortho="right",
    tags=commontags(linkinds(A, i)),
    cutoff=cutoff,
    maxdim=maxdim,
    mindim=mindim,
    kwargs...,
  )
  truncate!(C; kwargs...)
  return C
end

"""
    apply(A::MPO, B::MPO; kwargs...)

Contract the `MPO` `A'` with the `MPO` `B` and then map the prime level of the resulting
MPO back to having pairs of indices with prime levels of 1 and 0.

Equivalent to `replaceprime(contract(A', B; kwargs...), 2 => 1)`.

See also [`contract`](@ref) for details about the arguments available.
"""
function apply(A::MPO, B::MPO; kwargs...)
  AB = contract(A', B; kwargs...)
  return replaceprime(AB, 2 => 1)
end

(A::MPO)(B::MPO; kwargs...) = apply(A, B; kwargs...)

contract_mpo_mpo_doc = """
    contract(A::MPO, B::MPO; kwargs...) -> MPO
    *(::MPO, ::MPO; kwargs...) -> MPO

Contract the `MPO` `A` with the `MPO` `B`, returning an `MPO` with the
site indices that are not shared between `A` and `B`.

If you are contracting two MPOs with the same sets of indices, likely you
want to call something like:
```julia
C = contract(A', B; cutoff=1e-12)
C = replaceprime(C, 2 => 1)
```
That is because if MPO `A` has the index structure `-s'-A-s-` and MPO `B`
has the Index structure `-s'-B-s-`, if we only want to contract over
on set of the indices, we would do `(-s'-A-s-)'-s'-B-s- = -s''-A-s'-s'-B-s- = -s''-C-s-`,
and then map the prime levels back to pairs of primed and unprimed indices with:
`replaceprime(-s''-C-s-, 2 => 1) = -s'-C-s-`.

Since this is a common use case, you can use the convenience function:
```julia
C = apply(A, B; cutoff=1e-12)
```
which is the same as the code above.

# Keywords
- `cutoff::Float64=1e-13`: the cutoff value for truncating the density matrix eigenvalues. Note that the default is somewhat arbitrary and subject to change, in general you should set a `cutoff` value.
- `maxdim::Int=maxlinkdim(A) * maxlinkdim(B))`: the maximal bond dimension of the results MPS.
- `mindim::Int=1`: the minimal bond dimension of the resulting MPS.

See also [`apply`](@ref) for details about the arguments available.
"""

@doc """
$contract_mpo_mpo_doc
""" contract(::MPO, ::MPO)

*(A::MPO, B::MPO; kwargs...) = contract(A, B; kwargs...)

# TODO: try this to copy the docstring
# Causing an error in Revise
#@doc """
#$contract_mpo_mpo_doc
#""" *(::MPO, ::MPO)

#@doc (@doc contract(::MPO, ::MPO)) *(::MPO, ::MPO)

"""
    sample(M::MPO)

Given a normalized MPO `M`,
returns a `Vector{Int}` of `length(M)`
corresponding to one sample of the
probability distribution defined by the MPO,
treating the MPO as a density matrix.

The MPO `M` should have an (approximately)
positive spectrum.
"""
function sample(M::MPO)
  N = length(M)
  s = siteinds(M)
  R = Vector{ITensor}(undef, N)
  R[N] = M[N] * δ(dag(s[N]))
  for n in reverse(1:(N - 1))
    R[n] = M[n] * δ(dag(s[n])) * R[n + 1]
  end

  if abs(1.0 - R[1][]) > 1E-8
    error("sample: MPO is not normalized, norm=$(norm(M[1]))")
  end

  result = zeros(Int, N)
  ρj = M[1] * R[2]
  Lj = ITensor()

  for j in 1:N
    s = siteind(M, j)
    d = dim(s)
    # Compute the probability of each state
    # one-by-one and stop when the random
    # number r is below the total prob so far
    pdisc = 0.0
    r = rand()
    # Will need n, An, and pn below
    n = 1
    projn = ITensor()
    pn = 0.0
    while n <= d
      projn = ITensor(s)
      projn[s => n] = 1.0
      pnc = (ρj * projn * prime(projn))[]
      if imag(pnc) > 1e-8
        @warn "In sample, probability $pnc is complex."
      end
      pn = real(pnc)
      pdisc += pn
      (r < pdisc) && break
      n += 1
    end
    result[j] = n
    if j < N
      if j == 1
        Lj = M[j] * projn * prime(projn)
      elseif j > 1
        Lj = Lj * M[j] * projn * prime(projn)
      end
      if j == N - 1
        ρj = Lj * M[j + 1]
      else
        ρj = Lj * M[j + 1] * R[j + 2]
      end
      s = siteind(M, j + 1)
      normj = (ρj * δ(s', s))[]
      ρj ./= normj
    end
  end
  return result
end

function HDF5.write(parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, M::MPO)
  g = create_group(parent, name)
  attributes(g)["type"] = "MPO"
  attributes(g)["version"] = 1
  N = length(M)
  write(g, "rlim", M.rlim)
  write(g, "llim", M.llim)
  write(g, "length", N)
  for n in 1:N
    write(g, "MPO[$(n)]", M[n])
  end
end

function HDF5.read(parent::Union{HDF5.File,HDF5.Group}, name::AbstractString, ::Type{MPO})
  g = open_group(parent, name)
  if read(attributes(g)["type"]) != "MPO"
    error("HDF5 group or file does not contain MPO data")
  end
  N = read(g, "length")
  rlim = read(g, "rlim")
  llim = read(g, "llim")
  v = [read(g, "MPO[$(i)]", ITensor) for i in 1:N]
  return MPO(v, llim, rlim)
end


text:
using Test

# this line because the const objects are not being exported
include(joinpath(dirname(@__FILE__), "..", "src/.timeseriesrc.jl"))

@testset "timeseriesrc" begin

@testset "const values are set the package defaults" begin
    @testset "DECIMALS" begin
        @test DECIMALS == 4
    end

    @testset "MISSING" begin
        @test MISSING == NAN
    end
end

@testset "const values are correct" begin
    @testset "NAN" begin
        @test NAN == "NaN"
    end

    @testset "NA" begin
        @test NA == "NA"
    end

    @testset "BLACKHOLE" begin
        @test BLACKHOLE == "\u2B24"
    end

    @testset "DOTCIRCLE" begin
        @test DOTCIRCLE == "\u25CC"
    end

    @testset "QUESTION" begin
        @test QUESTION == "\u003F"
    end
end

end  # @testset "timeseriesrc"

ext: jl


text:
# From https://github.com/vchuravy/GPUifyLoops.jl/blob/b4d43198ec4aee75e6bcfdda7165ed5e469d2df5/src/loopinfo.jl
module LoopInfo

const HAS_LOOPINFO_EXPR = VERSION >= v"1.2.0-DEV.462"
export @unroll

##
# Uses the loopinfo expr node to attach LLVM loopinfo to loops
# the full list of supported metadata nodes is available at
# https://llvm.org/docs/LangRef.html#llvm-loop
# TODO: Figure out how to deal with compile-time constants in `@unroll(N, expr)`
#       so constants that come from `Val{N}` but are not parse time constant.
#       Most likely will require changes to base Julia.
##

module MD
    unroll_count(n) = (Symbol("llvm.loop.unroll.count"), convert(Int, n))
    unroll_disable() = (Symbol("llvm.loop.unroll.disable"), 1)
    unroll_enable() = (Symbol("llvm.loop.unroll.enable"), 1)
    unroll_full() = (Symbol("llvm.loop.unroll.full"), 1)
end

function loopinfo(expr, nodes...)
    if expr.head != :for
        error("Syntax error: loopinfo needs a for loop")
    end
    if HAS_LOOPINFO_EXPR
        push!(expr.args[2].args, Expr(:loopinfo, nodes...))
    end
    return expr
end

"""
   @unroll expr
Takes a for loop as `expr` and informs the LLVM unroller to fully unroll it, if
it is safe to do so and the loop count is known.
"""
macro unroll(expr)
    expr = loopinfo(expr, MD.unroll_full())
    return esc(expr)
end

"""
    @unroll N expr
Takes a for loop as `expr` and informs the LLVM unroller to unroll it `N` times,
if it is safe to do so.
"""
macro unroll(N, expr)
    if !(N isa Integer)
        error("Syntax error: `@unroll N expr` needs a constant integer N")
    end
    expr = loopinfo(expr, MD.unroll_count(N))
    return esc(expr)
end

end #module


text:
################################################################################
# Maptables
################################################################################

"""
TODO

API
---
* [`location_type`](@ref)
* [`getlocations`](@ref)
* [`isvalid`](@ref)

Implementations
---------------
* [`MapTable`](@ref)
"""
abstract type AbstractMapTable end

#######
# API #
#######

"""
Return the stored location type for a `MapTable`.

Method List
-----------
$(METHODLIST)
"""
function location_type end

"""
    getlocations(maptable, class::Int)

Return a vector of locations that nodes of type `class` can occupy.

Method List
-----------
$(METHODLIST)
"""
function getlocations end

"""
    isvalid(maptable, class, address :: Address)

Return `true` if nodes of type `class` can occupy `address`. In other words,
there is some component at `adddress` that class can be mapped to.

Method List
-----------
$(METHODLIST)
"""
function isvalid end

################################################################################
# MapTable
################################################################################

"""
Default implementation of [`AbstractMapTable`](@ref)

Parameters:
* `D` - The dimensionality of the `Addresses` in the table.
"""
struct MapTable{D} <: AbstractMapTable
    """
    Bit mask of whether a node class may be mapped to an address.

    Accessing strategy: Look up the node class to index the outer vector. Index the inner
    array with an address.

    A `true` entry means the task class can be mapped.
    A `false` entry means the task calss cannot be mapped.
    """
    mask::Vector{Array{Bool,D}}
end

function MapTable(
    toplevel::TopLevel, ruleset::RuleSet, equivalence_classes, pathtable::PathTable{D}
) where {D}

    # For each node class C, create an array the same size as the
    # pathtable. For each address A, record the indices of the paths at
    # componant_table[A] that C can be mapped to.
    mask = map(equivalence_classes.reps) do node
        map(pathtable) do path
            return path !== nothing && canmap(ruleset, node, toplevel[path])
        end
    end
    return MapTable(mask)
end
location_type(::MapTable{D}) where {D} = Address{D}

# Makes repeated calls to the 3 argument version of getlocations starting at
# all base addresses to generate ALL locations for this class.
function getlocations(maptable::MapTable{D}, class::Int) where {D}
    # Allocation an empty vector of the appropriate location type.
    addresses = location_type(maptable)[]
    table = maptable.mask[class]

    for address in CartesianIndices(table)
        isvalid(maptable, class, address) && push!(addresses, address)
    end
    return addresses
end

function isvalid(maptable::MapTable, class::Integer, address::Address)
    return maptable.mask[class][address]
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ c9cd29ea-906d-11eb-390f-87f6b15c1e04
begin
using Images, TestImages
using ImageQualityIndexes

img = testimage("cameraman") .|> float64
noisy_img = img .+ 0.1 .* randn(size(img))
assess_ssim(noisy_img, img) # 0.24112
assess_psnr(noisy_img, img) # 19.9697

kernel = ones(3, 3)./9 # mean filter
denoised_img = imfilter(noisy_img, kernel)
assess_psnr(denoised_img, img) # 28.4249
assess_ssim(denoised_img, img) # 0.6390
assess_msssim(denoised_img, img) # 0.8460

img = testimage("fabio");
colorfulness(img) # 68.5530
end

# ╔═╡ 803452c2-9063-11eb-1aa6-43489226ffa0
begin
# using ImageFeatures, TestImages, Images, ImageDraw, CoordinateTransformations,Rotations
#BRIEF
# img = testimage("lena_gray_512");
# img1 = Gray.(img);
# trans = Translation(-100, -200)
# 	img2 = warp(img1, trans, axes(img1));
# 	keypoints_1 = Keypoints(fastcorners(img1, 12, 0.4))
# 	keypoints_2 = Keypoints(fastcorners(img2, 12, 0.4))
# 	brief_params = BRIEF(size = 256, window = 10, seed = 123)

# desc_1, ret_keypoints_1 = create_descriptor(img1, keypoints_1, brief_params);
# desc_2, ret_keypoints_2 = create_descriptor(img2, keypoints_2, brief_params);

# matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.1)

# grid = hcat(img1, img2)
# offset = CartesianIndex(0, size(img1, 2))
# data=map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)
end

# ╔═╡ 5268c78c-9064-11eb-12d5-6d8c1b779872
# #ORB
# begin
# img = testimage("lighthouse")
# img1 = Gray.(img)
# rot = recenter(RotMatrix(5pi/6), [size(img1)...] .÷ 2)  # a rotation around the center
# tform = rot ∘ Translation(-50, -40)
# img2 = warp(img1, tform, axes(img1))
# orb_params = ORB(num_keypoints = 1000)
# desc_1, ret_keypoints_1 = create_descriptor(img1, orb_params)
# desc_2, ret_keypoints_2 = create_descriptor(img2, orb_params)
# matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.2)
# nothing # hide
# grid = hcat(img1, img2)
# offset = CartesianIndex(0, size(img1, 2))
# map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)
# end

# ╔═╡ 05b8e934-9065-11eb-02c4-652584ee0c8c
# begin
# #BRISK
# img = testimage("lighthouse")
# img1 = Gray.(img)
# rot = recenter(RotMatrix(5pi/6), [size(img1)...] .÷ 2)  # a rotation around the center
# tform = rot ∘ Translation(-50, -40)
# img2 = warp(img1, tform, axes(img1))

# features_1 = Features(fastcorners(img1, 12, 0.35))
# features_2 = Features(fastcorners(img2, 12, 0.35))

# brisk_params = BRISK()

# desc_1, ret_features_1 = create_descriptor(img1, features_1, brisk_params)
# desc_2, ret_features_2 = create_descriptor(img2, features_2, brisk_params)

# matches = match_keypoints(Keypoints(ret_features_1), Keypoints(ret_features_2), desc_1, desc_2, 0.1)

# grid = hcat(img1, img2)
# offset = CartesianIndex(0, size(img1, 2))
# map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)
# end

# ╔═╡ b7c4ec72-9065-11eb-3f25-b5008729692a
# begin
# 	#freak
# img = testimage("lighthouse")
# img1 = Gray.(img)
# rot = recenter(RotMatrix(5pi/6), [size(img1)...] .÷ 2)  # a rotation around the center
# tform = rot ∘ Translation(-50, -40)
# img2 = warp(img1, tform, axes(img1))

# keypoints_1 = Keypoints(fastcorners(img1, 12, 0.35))
# keypoints_2 = Keypoints(fastcorners(img2, 12, 0.35))

# freak_params = FREAK()

# desc_1, ret_keypoints_1 = create_descriptor(img1, keypoints_1, freak_params)
# desc_2, ret_keypoints_2 = create_descriptor(img2, keypoints_2, freak_params)

# matches = match_keypoints(ret_keypoints_1, ret_keypoints_2, desc_1, desc_2, 0.1)

# grid = hcat(img1, img2)
# offset = CartesianIndex(0, size(img1, 2))
# map(m -> draw!(grid, LineSegment(m[1], m[2] + offset)), matches)
# end

# ╔═╡ 99a6d316-9067-11eb-3b83-4f467da9f6c9
# begin
# #GLCM
# img = testimage("lighthouse")
# img = Gray.(img)
# glcm1 = glcm(img, 5, pi)
# end

# ╔═╡ 9532fc2e-9067-11eb-2f3f-9bb43fd1d6bb
# Gray.(glcm1/255)

# ╔═╡ 3b2289d8-9068-11eb-17ce-7f5b27f7bb99
# prop = glcm_prop(glcm1, 500, 500, glcm_entropy(zeros(5,5)))

# ╔═╡ Cell order:
# ╠═803452c2-9063-11eb-1aa6-43489226ffa0
# ╠═5268c78c-9064-11eb-12d5-6d8c1b779872
# ╠═05b8e934-9065-11eb-02c4-652584ee0c8c
# ╠═b7c4ec72-9065-11eb-3f25-b5008729692a
# ╠═99a6d316-9067-11eb-3b83-4f467da9f6c9
# ╠═9532fc2e-9067-11eb-2f3f-9bb43fd1d6bb
# ╠═3b2289d8-9068-11eb-17ce-7f5b27f7bb99
# ╠═c9cd29ea-906d-11eb-390f-87f6b15c1e04

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
abstract type AbstractGrid{C <: AbstractCell, N <: Tuple} end

floattype(::Type{<:AbstractGrid{C}}) where {C} = floattype(C)
arraytype(::Type{<:AbstractGrid{C}}) where {C} = arraytype(C)
celltype(::Type{<:AbstractGrid{C}}) where {C} = C
Base.ndims(::Type{<:AbstractGrid{C}}) where {C} = ndims(C)
Base.size(::Type{<:AbstractGrid{C, N}}) where {C, N} = size_to_tuple(N)
Base.length(::Type{<:AbstractGrid{C, N}}) where {C, N} = tuple_prod(N)

floattype(grid::AbstractGrid) = floattype(typeof(grid))
arraytype(grid::AbstractGrid) = arraytype(typeof(grid))
celltype(grid::AbstractGrid) = celltype(typeof(grid))
Base.ndims(grid::AbstractGrid)= Base.ndims(typeof(grid))
Base.size(grid::AbstractGrid) = Base.size(typeof(grid))
Base.length(grid::AbstractGrid) = Base.length(typeof(grid))

function WriteVTK.vtk_grid(filename::AbstractString, grid::AbstractGrid,
                           args...; kwargs...)
    vtk = vtk_grid(filename, points_vtk(grid), cells_vtk(grid), args...;
                   kwargs...)
    data_vtk!(vtk, grid)

    return vtk
end

struct NodalGrid{C <: AbstractCell, N <: Tuple, V, Y, P, Q, F, G, H, B} <: AbstractGrid{C, N}
    referencecell::C
    vertices::V
    connectivity::Y
    points::P
    metrics::Q
    faces::F
    faceindices::G
    facemetrics::H
    boundaryfaces::B
    type::Symbol
end

function Adapt.adapt_structure(to, grid::NodalGrid{C1, N}) where {C1, N}
    names = fieldnames(NodalGrid)
    args = ntuple(j->adapt(to, getfield(grid, names[j])), length(names))
    C2 = typeof(args[1])
    NodalGrid{C2, N, typeof.(args[2:end-1])...}(args...)
end

function NodalGrid(warp::Function, referencecell, vertices, connectivity, type=:unknown;
                   faces=nothing, boundaryfaces=nothing)
    C = typeof(referencecell)
    N = size(connectivity)
    V = typeof(vertices)
    Y = typeof(connectivity)

    unwarpedbrick = warp === identity && type === :brick

    points = materializepoints(referencecell, vertices, connectivity)
    points = warp.(points)
    P = typeof(points)

    metrics, facemetrics = materializemetrics(referencecell, points, unwarpedbrick)
    Q = typeof(metrics)
    H = typeof(facemetrics)

    if isnothing(faces)
        faces = materializefaces(referencecell, connectivity)
    end
    F = typeof(faces)

    faceindices = materializefaceindices(referencecell, faces)
    G = typeof(faceindices)

    if isnothing(boundaryfaces)
        boundaryfaces = materializeboundaryfaces(referencecell, faces)
    end
    B = typeof(boundaryfaces)

    types = (V, Y, P, Q, F, G, H, B)
    return NodalGrid{C, Tuple{N...}, types...}(referencecell, vertices,
                                               connectivity, points, metrics,
                                               faces, faceindices, facemetrics,
                                               boundaryfaces, type)
end

function NodalGrid(referencecell, vertices, connectivity, type=:unknown;
                   faces=nothing, boundaryfaces=nothing)
    return NodalGrid(identity, referencecell, vertices, connectivity, type;
                     faces=faces, boundaryfaces=boundaryfaces)
end

referencecell(grid::NodalGrid) = grid.referencecell
vertices(grid::NodalGrid) = grid.vertices
connectivity(grid::NodalGrid) = grid.connectivity
points(grid::NodalGrid) = grid.points
metrics(grid::NodalGrid) = grid.metrics
faces(grid::NodalGrid) = grid.faces
faceindices(grid::NodalGrid) = grid.faceindices
facemetrics(grid::NodalGrid) = grid.facemetrics
boundaryfaces(grid::NodalGrid) = grid.boundaryfaces

function points_vtk(grid::NodalGrid)
    P = toequallyspaced(referencecell(grid))
    x = P * points(grid)

    return reinterpret(reshape, floattype(grid), vec(adapt(Array, x)))
end

function cells_vtk(grid::NodalGrid)
    type = celltype_vtk(referencecell(grid))
    connectivity = connectivity_vtk(referencecell(grid))

    cells = [MeshCell(type, e * length(connectivity) .+ connectivity)
             for e = 0:length(grid)-1]

    return cells
end

function data_vtk!(vtk, grid::NodalGrid)
    higherorderdegrees = zeros(Int, 3, length(grid))
    ds = [degrees(referencecell(grid))...]
    higherorderdegrees[1:length(ds), :] .= repeat(ds, 1, length(grid))

    vtk["HigherOrderDegrees", VTKCellData()] = higherorderdegrees

    return
end

function materializefaces(referencecell::AbstractCell, connectivity)
    cellfaces = materializefaces(referencecell)[2:end]
    return ntuple(i->connect(cellfaces[i], connectivity), length(cellfaces))
end

function connect(cellfaces, connectivity)
    A = arraytype(connectivity)
    # TODO Should we move this calculation to the GPU?
    connectivity = adapt(Array, connectivity)
    connectivity = reinterpret(reshape, Int, vec(connectivity))

    faces = Array{Int}(undef, size(cellfaces)..., size(connectivity, 2))
    @tullio faces[i, j, k] = connectivity[cellfaces[i, j], k]
    faces = vec(reinterpret(reshape, SVector{size(cellfaces, 1), Int}, faces))

    localfacenumbers = 1:length(faces)
    globalfacenumbers = numbercontiguous(faces; by=sort)
    facepermutations = tuplesortpermutation.(Tuple.(faces))

    M = sparse(localfacenumbers, globalfacenumbers, facepermutations)
    if A <: CuArray
        M = adapt(A, GeneralSparseMatrixCSC(M))
    end
    return M
end

struct FaceConnectionException <: Exception end

@kernel function materializefaceindices!(referencecell::AbstractCell,
                                         faceindices⁺, faceindices⁻,
                                         @Const(globaltolocalfaces))
    @uniform begin
        localfaces = rowvals(globaltolocalfaces)
        vertexpermutations = nonzeros(globaltolocalfaces)
        num_faces = number_of_faces(referencecell)[2]
        num_dof_per_cell = length(referencecell)
        conn = connectivity(referencecell)[2]
        offsets = connectivityoffsets(referencecell, Val(2))
    end

    i, j = @index(Global, NTuple)
    r = nzrange(globaltolocalfaces, j)

    r₁ = first(r)
    r₂ = last(r)

    if r₂ - r₁ ≥ 2 || r₂ - r₁ < 0
        throw(FaceConnectionException())
    end

    e₁, f₁ = divrem(localfaces[r₁] - 1, num_faces) .+ 1
    e₂, f₂ = divrem(localfaces[r₂] - 1, num_faces) .+ 1

    p₂₁ = vertexpermutations[r₁] ∘ inv(vertexpermutations[r₂])
    p₁₂ = vertexpermutations[r₂] ∘ inv(vertexpermutations[r₁])

    if i ≤ length(conn[f₁])
        i₁ = i + offsets[f₁]
        i₂ = i + offsets[f₂]

        j₁ = conn[f₁][i] + (e₁ - 1) * num_dof_per_cell
        j₂ = conn[f₂][i] + (e₂ - 1) * num_dof_per_cell

        j₂₁ = getpermutedindex(conn[f₁], p₂₁, i) .+ (e₁ - 1) * num_dof_per_cell
        j₁₂ = getpermutedindex(conn[f₂], p₁₂, i) .+ (e₂ - 1) * num_dof_per_cell

        faceindices⁻[i₁, e₁] = j₁
        faceindices⁻[i₂, e₂] = j₂

        faceindices⁺[i₁, e₁] = j₁₂
        faceindices⁺[i₂, e₂] = j₂₁
    end
end

function materializefaceindices(referencecell::AbstractCell, faces)
    C = typeof(referencecell)
    A = arraytype(C)
    globaltolocalfaces = first(faces)
    num_localfaces, num_globalfaces = size(globaltolocalfaces)
    num_cells = div(num_localfaces, number_of_faces(C)[2])
    faceconn = connectivity(referencecell)[2]
    max_num_indices_per_face = maximum(length.(faceconn))
    num_facesindices = sum(length.(faceconn))

    faceindices⁺ = fill!(A{Int, 2}(undef, num_facesindices, num_cells), 0)
    faceindices⁻ = fill!(A{Int, 2}(undef, num_facesindices, num_cells), 0)

    kernel = materializefaceindices!(device(A), (max_num_indices_per_face, 5))
    event = Event(device(A))
    event = kernel(referencecell, faceindices⁺, faceindices⁻,
                   globaltolocalfaces;
                   ndrange = (max_num_indices_per_face, num_globalfaces),
                   dependencies = (event,))
    wait(event)

    inidices = (faceindices⁻, faceindices⁺)
    return inidices
end

@kernel function materializeboundaryfaces!(referencecell::AbstractCell,
                                           boundaryfaces,
                                           @Const(globaltolocalfaces))
    @uniform begin
        localfaces = rowvals(globaltolocalfaces)
        vertexpermutations = nonzeros(globaltolocalfaces)
        num_faces = number_of_faces(referencecell)[2]
        num_dof_per_cell = length(referencecell)
        conn = connectivity(referencecell)[2]
        offsets = connectivityoffsets(referencecell, Val(2))
    end

    j = @index(Global)
    r = nzrange(globaltolocalfaces, j)

    r₁ = first(r)
    r₂ = last(r)

    if r₂ - r₁ ≥ 2 || r₂ - r₁ < 0
        throw(FaceConnectionException())
    end

    e₁, f₁ = divrem(localfaces[r₁] - 1, num_faces) .+ 1
    e₂, f₂ = divrem(localfaces[r₂] - 1, num_faces) .+ 1

    if r₁ == r₂
        boundaryfaces[f₁, e₁] = 1
    else
        boundaryfaces[f₁, e₁] = 0
        boundaryfaces[f₂, e₂] = 0
    end
end

function materializeboundaryfaces(referencecell, faces)
    A = arraytype(referencecell)
    globaltolocalfaces = first(faces)
    num_localfaces, num_globalfaces = size(globaltolocalfaces)
    num_cells = div(num_localfaces, number_of_faces(referencecell)[2])
    num_faces = number_of_faces(referencecell)[2]

    boundaryfaces = fill!(A{Int, 2}(undef, num_faces, num_cells), 0)

    kernel = materializeboundaryfaces!(device(A), (256,))
    event = Event(device(A))
    event = kernel(referencecell, boundaryfaces, globaltolocalfaces;
                   ndrange = (num_globalfaces,), dependencies = (event,))
    wait(event)

    return boundaryfaces
end

function faceviews(referencecell, A::AbstractMatrix)
    num_faces = number_of_faces(referencecell)[2]
    offsets = connectivityoffsets(referencecell, Val(2))
    facesizes = size.(connectivity(referencecell)[2])

    if last(offsets) != size(A, 1)
        throw(ArgumentError("The first dimension of A needs to contain the face degrees of freedom."))
    end

    return ntuple(Val(num_faces)) do f
        reshape(view(A, (1+offsets[f]):offsets[f+1], :), facesizes[f]..., :)
    end
end

function min_node_distance(grid::NodalGrid; dims = 1:ndims(grid))
    @assert maximum(dims) <= ndims(grid)

    A = arraytype(grid)
    T = floattype(grid)
    cell = referencecell(grid)

    min_neighbour_distance = similar(points(grid), T)

    Np = length(cell)
    event = min_neighbour_distance_kernel(device(A), 256)(
        min_neighbour_distance, points(grid), Val(strides(cell)), Val(Np), Val(dims);
        ndrange = length(grid) * length(cell))
    wait(event)

    minimum(min_neighbour_distance)
end

@kernel function min_neighbour_distance_kernel(min_neighbour_distance, points,
                                               ::Val{S}, ::Val{Np}, ::Val{dims}) where {S, Np, dims}
    I = @index(Global, Linear)

    @inbounds begin
        e = (I - 1) ÷ Np + 1
        ijk = (I - 1) % Np + 1

        md = typemax(eltype(min_neighbour_distance))
        x⃗ = points[ijk, e]
        for d in dims
            for m in (-1, 1)
                ijknb = ijk + S[d] * m
                    if 1 <= ijknb <= Np
                        x⃗nb = points[ijknb, e]
                        md = min(norm(x⃗ - x⃗nb), md)
                    end
            end
        end
        min_neighbour_distance[ijk, e] = md
    end
end


text:
import Elemental
const El = Elemental
using LinearAlgebra: mul!

n0 = 50
n1 = 50
testNative = true
display = false
worldRank = El.MPI.commRank(El.MPI.CommWorld[])
worldSize = El.MPI.commSize(El.MPI.CommWorld[])
if worldRank == 0
  print("worldSize=$worldSize\n")
end

using DistributedArrays

function randDArray(m,n,sparsity=0.01)
    DA = DistributedArrays.DArray((m,n)) do I
        sprandn(length(I[1]), length(I[2]), sparsity)
    end
    A = El.DistSparseMatrix{Float64}(DA)
    return A
end

# Stack two 2D finite-difference matrices on top of each other
# and make the last column dense
function stackedFD2D(n0, n1)
    height = 2*n0*n1
    width = n0*n1
    A = El.DistSparseMatrix(Float64, height, width)
    localHeight = El.localHeight(A)
    El.reserve(A, 6*localHeight)

    for sLoc in 1:localHeight
        s = El.globalRow(A, sLoc)
        if s <= n0*n1
            x0 = ((s-1) % n0) + 1
            x1 = div((s-1), n0) + 1
            El.queueUpdate(A, s, s, 11.0)
            if x0 > 1
                El.queueUpdate(A, s, s - 1, -10.0)
            end
            if x0 < n0
                El.queueUpdate(A, s, s + 1, 20.0)
            end
            if x1 > 1
                El.queueUpdate(A, s, s - n0, -30.0)
            end
            if x1 < n1
                El.queueUpdate(A, s, s + n0, 40.0)
            end
        else
            sRel = s - n0*n1
            x0 = ((sRel-1) % n0) + 1
            x1 = div(sRel-1, n0) + 1
            El.queueUpdate(A, s, sRel, -20.0)
            if x0 > 1
                El.queueUpdate(A, s, sRel - 1, -1.0)
            end
            if x0 < n0
                El.queueUpdate(A, s, sRel + 1, -2.0)
            end
            if x1 > 1
                El.queueUpdate(A, s, sRel - n0, -3.0)
            end
            if x1 < n1
                El.queueUpdate(A, s, sRel + n0, 3.0)
            end
        end

        # The dense last column
        El.queueUpdate(A, s, width, floor(-10/height))
    end
    El.processQueues(A)
    return A
end

if testNative
    A = stackedFD2D(n0, n1)
else
    A = randDArray(2*n0*n1,n0*n1)
end

b = El.DistMultiVec(Float64)
El.gaussian!(b, 2*n0*n1, 1)

AHeight = El.height(A)
AWidth = El.width(A)
bHeight = El.height(b)
bWidth = El.width(b)

if display
    El.print(A, "A")
    # El.print(b, "b")
end

ctrl = El.LPAffineCtrl(Float64,
    mehrotraCtrl=El.MehrotraCtrl(Float64,
        solveCtrl=El.RegSolveCtrl(Float64,
                                  progress=true),
        print=true,
        outerEquil=true,
        time=true)
    )

# elapsedLAV = @elapsed x = El.lav(A, b);#Elemental.print(A, "Matrix A")
elapsedLAV = @elapsed x = El.lav(A, b, ctrl)

if El.MPI.commRank(El.MPI.CommWorld[]) == 0
    println("LAV time: $elapsedLAV seconds")
end
bTwoNorm = El.nrm2(b)
bInfNorm = El.maxNorm(b)

r = copy(b)
mul!(r, A, x, -1.0, 1.0)

rTwoNorm = El.nrm2(r)
rOneNorm = El.entrywiseNorm(r, 1)

if El.MPI.commRank(El.MPI.CommWorld[]) == 0
    println("|| b ||_2       = $bTwoNorm")
    println("|| b ||_oo      = $bInfNorm")
    println("|| A x - b ||_2 = $rTwoNorm")
    println("|| A x - b ||_1 = $rOneNorm")
end

elapsedLS = @elapsed xLS = El.leastSquares(A, b)

if El.MPI.commRank(El.MPI.CommWorld[]) == 0
    println("LS time: $elapsedLAV seconds")
end

rLS = copy(b)
mul!(rLS, A, xLS, -1.0, 1.0)
if display
    El.print( rLS, "A x_{LS} - b" )
end

rLSTwoNorm = El.nrm2(rLS)
rLSOneNorm = El.entrywiseNorm(rLS, 1)
if El.MPI.commRank(El.MPI.CommWorld[]) == 0
    println("|| A x_{LS} - b ||_2 = $rLSTwoNorm")
    println("|| A x_{LS} - b ||_1 = $rLSOneNorm")
end

elapsedRidge = @elapsed xRidge = El.ridge(A, b, 0.0)

if El.MPI.commRank(El.MPI.CommWorld[]) == 0
    println("Ridge time: $elapsedRidge seconds")
end

rRidge = copy(b)
mul!(rRidge, A, xRidge, -1.0, 1.0)
if display
    El.print( rRidge, "A x_{Ridge} - b" )
end

rRidgeTwoNorm = El.nrm2(rRidge)
rRidgeOneNorm = El.entrywiseNorm(rRidge, 1)
if El.MPI.commRank(El.MPI.CommWorld[]) == 0
    println("|| A x_{Ridge} - b ||_2 = $rRidgeTwoNorm")
    println("|| A x_{Ridge} - b ||_1 = $rRidgeOneNorm")
end

# Require the user to press a button before the figures are closed
# commSize = El.mpi.Size( El.mpi.COMM_WORLD() )
# El.Finalize()


text:
export add_meta!, locations, locations!, loci, genotypes, get_genotypes, get_genotype, populations, population, populations!, population!, exclude, remove, omit, exclude!, remove!, omit!, samples

"""
    add_meta!(popdata::PopData, metadata::T; name::String, loci::Bool = true, categorical::Bool = true) where T <: AbstractVector
Add an additional metadata information to a `PopData` object. Mutates `PopData` in place. Metadata
must be in the same order as the samples in `PopData.meta`.

#### Arguments
- `popdata` : The `PopData` object to add information to
- `metadata` : A `Vector` with the metadata you wish to add to the `PopData`, in the same order as the names appear in `PopData.meta`

#### Keyword Arguments
- `name` : String of the name of this new column
- `loci` : Boolean of whether to also add this information to `PopData.loci` (default: `true`)
- `categorical` : Boolean of whether the metadata being added is categorical aka "factors" (default: `true`)
"""
function add_meta!(popdata::PopData, metadata::T; name::String, loci::Bool = true, categorical::Bool = true) where T <: AbstractVector
    length(metadata) != length(popdata.meta.name) && error("Provided metadata vector (n = $length(metadata)) and samples in PopData (n = $length(popdata.meta.name)) have different lengths")
    @info "Adding $Symbol(name) column to .meta dataframe"
    # add to meta
    insertcols!(popdata.meta, Symbol(name) => metadata)

    # add to loci
    if loci
        @info "Adding $Symbol(name) column to .meta and .loci dataframes"
        tmp = DataFrame(:name => popdata.meta.name, Symbol(name) => metadata)
        popdata.loci = outerjoin(popdata.loci, tmp, on = :name)
        if categorical == true
            popdata.loci[name] = PooledArray(popdata.loci[name])
        end
    end
    return
end

"""
    add_meta!(popdata::PopData, samples::Vector{String}, metadata::T; name::String, loci::Bool = true, categorical::Bool = true) where T <: AbstractVector
Add an additional metadata information to a `PopData` object. Mutates `PopData` in place. Takes a vector of
sample names if the metadata is not in the same order as samples appear in `PopData.meta`.

#### Arguments
- `popdata` : The `PopData` object to add information to
- `sample` : A `Vector{String}` of sample names corresponding to the order of the `metadata`
- `metadata` : A `Vector` with the metadata you wish to add to the `PopData`, in the same order as the names appear in `PopData.meta`

#### Keyword Arguments
- `name` : String of the name of this new column
- `loci` : Boolean of whether to also add this information to `PopData.loci` (default: `true`)
- `categorical` : Boolean of whether the metadata being added is categorical aka "factors" (default: `true`)
"""
function add_meta!(popdata::PopData, samples::Vector{String}, metadata::T; name::String, loci::Bool = true) where T <: AbstractVector
    length(samples) != length(popdata.meta.name) && error("Provided sample vector (n = $length(samples)) and samples in PopData (n = $length(popdata.meta.name)) have different lengths")
    length(samples) != length(metadata) && error("Sample names (n = $length(samples)) and metadata vectors (n = $length(metadata)) have different lengths")
    sort(samples) != sort(popdata.meta.name) && error("Sample names are not identical")
    @info "Adding $Symbol(name) column to .meta dataframe"

    tmp = DataFrame(:name => samples, Symbol(name) => metadata)
    popdata.meta = outerjoin(popdata.meta, tmp, on = :name)

    if loci
        @info "Adding $Symbol(name) column to .meta and .loci dataframes"
        popdata.loci = outerjoin(popdata.loci, tmp, on = :name)
        if categorical == true
            popdata.loci[name] = PooledArray(popdata.loci[name])
        end
    end
    return
end

"""
    locations(data::PopData)
View the longitude and latitude data in a `PopData` object. Returns a table
derived from the PopData. Changes made to this table will not alter the source
`PopData` object.

Use `locations!` to add spatial data to a `PopData` object.
"""
function locations(data::PopData)
    @view data.meta[!, [:longitude, :latitude]]
end

"""
    locations!(data::PopData; long::Vector{Float64}, lat::Vector{Float64})
Replaces existing `PopData` location data (longitude `long`, latitude `lat`).
Takes **decimal degrees** as a `Vector` of any `AbstractFloat`.
## Formatting requirements
- Decimal Degrees format: `-11.431`
- **Must** use negative sign `-` instead of cardinal directions
- Location data must be in the order that samples appear in your `PopData`
- Missing data should be coded as `missing` values of type `Missing` (can be accomplished with `replace!()`)
### Example
```
ncats = @nancycats ;
x = rand(237) ; y = rand(237)
locations!(ncats, long = x, lat = y)
```
"""
function locations!(data::PopData, long::Vector{Union{Missing,T}}, lat::Vector{Union{Missing,T}}) where T <: AbstractFloat
    long_len = length(long)
    lat_len = length(lat)
    long_len != lat_len && error("latitude ($lat_len) and longitude ($long_len) arrays not equal in length")
    long_len != length(data.meta.name) && error("lat/long array length ($long_len) and number of samples in PopData ($long_len) are not equal")

    data.meta.longitude .= long
    data.meta.latitude .= lat
    return
end

function locations!(data::PopData, long::Vector{T}, lat::Vector{T}) where T <: AbstractFloat
    # convert to the right type and use locations!()
    lat_adjust = lat |> Vector{Union{Missing, Float32}}
    long_adjust = long |> Vector{Union{Missing, Float32}}
    locations!(data, long = long_adjust, lat = lat_adjust)
end

"""
    locations!(data::PopData; long::Vector{String}, lat::Vector{String})
Replaces existing `PopData` location data (longitude `long`, latitude `lat`). Takes
**decimal minutes** or **degrees minutes seconds** format as a `Vector` of `String`. Recommended to use `CSV.read`
from `CSV.jl` to import your spatial coordinates from a text file.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as the string `"missing"` (can be accomplished with `replace!()`)
- Can mix colons and spaces (although it's bad practice)
### NOTE
If you read in the coordinate data as 4 vectors (longitude degrees, longitude minutes, latitude degrees, latitude minutes),
then the easiest course of action would be to merge them into two vectors of strings
(one for longitude, one for latitude):
```
long_string = string.(lat_deg, " ", lat_min)
lat_string = string.(long_deg, " ", long_min)
```
and use these as inputs into `locations!`

### Example
```
ncats = @nancycats;
x = fill("11 22.33W", 237) ; y = fill("-41 31.52", 237)
locations!(ncats, long = x, lat = y)
```
"""
function locations!(data::PopData, long::Vector{String}, lat::Vector{String})
    long_len = length(long)
    lat_len = length(lat)
    lat_len != long_len && error("latitude ($lat_len) and longitude ($long_len) arrays not equal in length")
    lat_len != length(data.meta.name) && error("lat/long array length ($lat_len) and number of samples in PopData ($long_len) are not equal")
    @info "Converting decimal minutes to decimal degrees"
    # convert coordinates to decimal degrees
    data.meta.longitude .= convert_coord.(long)
    data.meta.latitude .= convert_coord.(lat)
    return
end

function locations!(data::PopData; kwargs...)
    kwargs = Dict(kwargs)
    # check for matching lat and long keywords
    if all([haskey(kwargs, :lat), haskey(kwargs, :long)])
        locations!(data, kwargs[:long], kwargs[:lat])
    else
        error("keyword arguments \"lat\" and \"long\" must be supplied together")
    end
end

"""
    loci(data::PopData)
Returns an array of strings of the loci names in a `PopData` object.
"""
function loci(data::PopData)
    unique(data.loci.locus)
end

"""
    get_genotype(data::PopObj; sample::String, locus::String)
Return the genotype of one sample at one locus in a `PopData` object.
### Example
```
cats = @nancycats;
get_genotype(cats, sample = "N115", locus = "fca8")
```
"""
function get_genotype(data::PopData; sample::String, locus::String)
    @views data.loci[(data.loci.name .== sample) .& (data.loci.locus .== locus), :genotype][1]
end

"""
    get_genotypes(data::PopData, sample::String)
Return a vector of all the genotypes of a sample in a `PopData` object. To return a
single genotype at a locus, see `get_genotype`.
```
cats = @nancycats
get_genotypes(cats, "N115")
```
"""
function get_genotypes(data::PopObj, sample::String)
    data.loci[data.loci.name .== sample, :genotype]
end

"""
    get_genotypes(data::PopObj; name::Union{String, Vector{String}}, locus::Union{String, Vector{String}})
Return a table of the genotype(s) of one or more `samples` for one or more
specific `loci` (both as keywords) in a `PopData` object.
### Examples
```
cats = @nancycats;
get_genotypes(cats, name = "N115" , locus = "fca8")
get_genotypes(cats, name = ["N115", "N7"] , locus = "fca8")
get_genotypes(cats, name = "N115" , locus = ["fca8", "fca37"])
get_genotypes(cats, name = ["N1", "N2"] , locus = ["fca8", "fca37"])
```
"""
function get_genotypes(data::PopData; name::Union{String, Vector{String}}, locus::Union{String, Vector{String}})
    if typeof(name) == String
        name = [name]
    end
    if typeof(locus) == String
        locus = [locus]
    end
    @view data.loci[(data.loci.name .∈ Ref(name)) .& (data.loci.locus .∈ Ref(locus)), :]
end

"""
    genotypes(data::PopData, locus::Union{String, Symbol})
Convenience wrapper to return a vector of all the genotypes of a single locus

### Example
```
genotypes(@gulfsharks, "contig_475")
```
"""
function genotypes(data::PopData, locus::String)
    @view data.loci[data.loci.locus .== locus, :genotype]
end

"""
    populations(data::PopData; listall::Bool = false)
View unique population ID's and their counts in a `PopData`.

- `listall = true` displays all samples and their `population` instead (default = `false`)
"""
@inline function populations(data::PopData; listall::Bool = false)
    if all(ismissing.(data.meta.population)) == true
        @info "no population data present in PopData"
        return
    elseif listall == true
        return @view data.meta[:, [:name, :population]]
    else
        pops = countmap(data.meta.population)
        return DataFrame(:population => collect(keys(pops)), :count => collect(values(pops)))
    end
end
"""
```
# Replace by matching
populations!(data::PopData, rename::Dict)
populations!(data::PopData, rename::Vector{String})
populations!(data::PopData, samples::Vector{String}, populations::Vector{String})
```
Multiple methods to rename or reassign population names to a `PopData`.

## Rename using a Dictionary
Rename existing population ID's of `PopData` using a `Dict` of
`population_name => replacement`
\n**Example**
```
potatopops = Dict("1" => "Idaho", "2" => "Russet")
populations!(potatoes, potatopops)
```

## Rename using a Vector of Strings
`Vector` of new unique population names in the order that they appear in the PopData.meta
\n**Example**
```
potatopops = ["Idaho", "Russet"]
populations!(potatoes, potatopops)
```
## Reassign using samples and new population assignments
Completely reassign populations for each individual. Takes two vectors of strings
as input: one of the sample names, and the other with their new corresponding
population name. This can be useful to change population names for only some individuals.

\n**Example**
```
populations!(potatoes, ["potato_1", "potato_2"], ["north_russet", "south_russet"])
```

"""
function populations!(data::PopData, rename::Dict)
    msg = ""
    @inbounds for key in keys(rename)
        if key ∉ unique(data.meta.population)
            msg *= "  Population \"$key\" not found in PopData\n"
        else
            replace!(data.meta.population, key => rename[key])
            replace!(data.loci.population.pool, key => rename[key])
        end
    end
    msg != "" && printstyled("Warnings:", color = :yellow) ; print("\n"*msg)
    return
end

function populations!(data::PopData, rename::Vector{String})
    current_popnames = unique(data.meta.population)
    rn_dict = Dict{String, String}()
    [rn_dict[string(i)] = j for (i,j) in zip(current_popnames, rename)]
    populations!(data, rn_dict)
    return
end

function populations!(data::PopData, samples::Vector{String}, populations::Vector{String})
    meta_df = groupby(data.meta, :name)
    loci_df = groupby(data.loci, :name)
    for (sample, new_pop) in zip(samples, populations)
        meta_df[(name = sample,)].population .= new_pop
        loci_df[(name = sample,)].population .= new_pop
    end
    # drop old levels
    data.loci.population = data.loci.population |> Array |> PooledArray
    return
end

const population = populations
const population! = populations!
const popnames! = populations!

##### Exclusion #####
"""
    exclude!(data::PopData, kwargs...)
Edit a `PopData` object in-place by excluding all occurences of the specified information.
The keywords can be used in any combination. Synonymous with `omit!` and `remove!`.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to remove from the `PopData`.
The keyword `loci` also works.

#### `population`
A `String` or `Vector{String}` of populations you want to remove from the `PopData`
The keyword `populations` also works.

#### `name`
A `String` or `Vector{String}` of samples you want to remove from the `PopData`
The keywords `names`, `sample`, and `samples` also work.

**Examples**
```
cats = @nancycats;
exclude!(cats, name = "N100", population = ["1", "15"])
exclude!(cats, samples = ["N100", "N102", "N211"], locus = ["fca8", "fca23"])
exclude!(cats, names = "N102", loci = "fca8", population = "3")
```
"""
function exclude!(data::PopData; kwargs...)
    filter_by = Dict(kwargs...)
    tmp = data
    filter_params = keys(filter_by) |> collect
    notices = ""
    # populations
    # check for keywords
    filt_pop = get.(Ref(filter_by), [:population, :populations], nothing)
    filt_pop = filt_pop[filt_pop .!== nothing]
    if length(filt_pop) != 0
        filt_pop = filt_pop[begin]
        if typeof(filt_pop) == String
            filt_pop = [filt_pop]
        end
        err = filt_pop[filt_pop .∉ Ref(unique(tmp.meta.population))]
        if length(err) > 0
            [notices *= "\n  population \"$i\" not found" for i in err]
        end
        # choose the cheaper method
        all_pops = unique(tmp.meta.population)
        if length(filt_pop) < length(all_pops)/2
            filter!(:population => x -> x ∉ filt_pop, tmp.loci)
            filter!(:population => x -> x ∉ filt_pop, tmp.meta)
        else
            keep = all_pops[all_pops .∉ Ref(filt_pop)]
            filter!(:population => x -> x in filt_pop, tmp.loci)
            filter!(:population => x -> x in filt_pop, tmp.meta)
        end
        tmp.loci.population = tmp.loci.population |> Array |> PooledArray
    end

    # samples
    # check for keywords
    filt_name = get.(Ref(filter_by), [:name, :names, :sample, :samples], nothing)
    filt_name = filt_name[filt_name .!== nothing]
    if length(filt_name) != 0
        filt_name = filt_name[begin]
        if typeof(filt_name) == String
            filt_name = [filt_name]
        end
        err = filt_name[filt_name .∉ Ref(tmp.meta.name)]
        if length(err) > 0
            [notices *= "\n  sample \"$i\" not found" for i in err]
        end
        # choose the cheaper method
        all_samples = tmp.meta.name
        if length(filt_name) < length(all_samples)/2
            filter!(:name => x -> x ∉ filt_name, tmp.loci)
            filter!(:name => x -> x ∉ filt_name, tmp.meta)
        else
            keep = all_samples[all_samples .∉ Ref(filt_name)]
            filter!(:name => x -> x in keep, tmp.loci)
            filter!(:name => x -> x in keep, tmp.meta)
        end
        tmp.loci.name = tmp.loci.name |> Array |> PooledArray
    end

    # loci
    # check for keywords
    filt_loci = get.(Ref(filter_by), [:locus, :loci], nothing)
    filt_loci = filt_loci[filt_loci .!= nothing]
    if length(filt_loci) != 0
        filt_loci = filt_loci[begin]
        if typeof(filt_loci) == String
            filt_loci = [filt_loci]
        end
        err = filt_loci[filt_loci .∉ Ref(loci(tmp))]
        if length(err) > 0
            [notices *= "\n  locus \"$i\" not found" for i in err]
        end
        # choose the cheaper method
        all_loci = loci(tmp)
        if length(filt_loci) < length(all_loci)/2
            filter!(:locus => x -> x ∉ filt_loci, tmp.loci)
        else
            keep = all_loci[all_loci .∉ Ref(filt_loci)]
            filter!(:locus => x -> x in keep, tmp.loci)
        end
        tmp.loci.locus = tmp.loci.locus |> Array |> PooledArray
    end

    # print the notices, if any
    if notices != ""
        printstyled("Notices:", bold = true, color = :blue)
        print(notices, "\n\n")
    end
    return tmp
end

const omit! = exclude!
const remove! = exclude!

"""
    exclude(data::PopData, kwargs...)
Returns a new `PopData` object excluding all occurrences of the specified keywords.
The keywords can be used in any combination. Synonymous with `omit` and `remove`.

### Keyword Arguments
#### `locus`
A `String` or `Vector{String}` of loci you want to remove from the `PopData`.
The keyword `loci` also works.

#### `population`
A `String` or `Vector{String}` of populations you want to remove from the `PopData`.
The keyword `populations` also works.

#### `name`
A `String` or `Vector{String}` of samples you want to remove from the `PopData`.
The keywords `names`, `sample`, and `samples` also work.

**Examples**
```
cats = @nancycats;
exclude(cats, name = "N100", population = ["1", "15"])
exclude(cats, samples = ["N100", "N102", "N211"], locus = ["fca8", "fca23"])
exclude(cats, names = "N102", loci = "fca8", population = "3")
```
"""
function exclude(data::PopData; kwargs...)
    tmp = copy(data)
    exclude!(tmp; kwargs...)
    return tmp
end

const omit = exclude
const remove = exclude

"""
    samples(data::PopData)
View individual/sample names in a `PopData`
"""
function samples(data::PopData)
    @view data.meta[!, :name]
end


text:
# %% Interactive GALI psos for henon heiles
# color coded with lyapunov exponent
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))

using DynamicalSystems, InteractiveDynamics
using OrdinaryDiffEq
import GLMakie

diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)

hh = Systems.henonheiles()

potential(x, y) = 0.5(x^2 + y^2) + (x^2*y - (y^3)/3)
energy(x,y,px,py) = 0.5(px^2 + py^2) + potential(x,y)
E = energy(get_state(hh)...)

function complete(y, py, x)
    V = potential(x, y)
    Ky = 0.5*(py^2)
    Ky + V ≥ E && error("Point has more energy!")
    px = sqrt(2(E - V - Ky))
    ic = [x, y, px, py]
    return ic
end

plane = (1, 0.0) # first variable crossing 0

cmap = cgrad(:cividis)
function λcolor(u)
    λ = lyapunov(hh, 10000; u0 = u, diffeq...)
    λmax = 0.06
    v = clamp(λ/λmax, 0, 1)
    return cmap.colors[round(Int, v*255 + 1)]
    # RGBf0(v/2, 0, v)
end

fig, state = interactive_poincaresos(
    hh, plane, (2, 4), complete;
    labels = ("q₂" , "p₂"), color = λcolor, diffeq...
)

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 10)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Distributions, CmdStan, Random, Parameters
ProjDir = @__DIR__
cd(ProjDir)
include("LNR.jl")

stream = open("LNR_Model.stan","r")
Model = read(stream,String)
close(stream)
seed = 343

stanmodel = Stanmodel(Sample(save_warmup=false, num_warmup=1000,
  num_samples=1000, thin=1), nchains=1, name="LNR", model=Model,
  printsummary=false, output_format=:mcmcchains, random=CmdStan.Random(seed))

Random.seed!(seed)
Nreps = 30
Nr = 3
N = 30

for i in 1:Nreps
    μ = rand(Normal(0 ,3), Nr)
    σ = rand(Uniform(.5, 2))
    dist = LNR(μ=μ, σ=σ, ϕ=0.0)
    data = rand(dist, N)
    choice = map(x->x[1], data)
    rt = map(x->x[2], data)
    stanInput = Dict("choice"=>choice, "rt"=>rt, "Nr"=>Nr, "N"=>N)
    rc, chain, cnames = stan(stanmodel, stanInput, ProjDir)
end


text:
function swsetupmatrix(s1, s2)
    setupmatrix = zeros(Int, length(s1)+1, length(s2)+1)
    return setupmatrix
end

function swscorematrix(seq1, seq2; match=1, mismatch=-1, gap=-1)
    scoremat = swsetupmatrix(seq1, seq2)
    for i in 2:size(scoremat, 1) # iterate through row indices
        for j in 2:size(scoremat, 2) # iterate through column indices
            above = nwscore(seq1[i-1], nothing; match, mismatch, gap) + scoremat[i-1,j]
            left = nwscore(nothing, seq2[j-1]; match, mismatch, gap) + scoremat[i,j-1]
            diagonal = nwscore(seq1[i-1], seq2[j-1]; match, mismatch, gap) + scoremat[i-1,j-1]
            scoremat[i,j] = max(above, left, diagonal)
            # or `scoremat[i,j] = max(above, left, diagonal, 0)`

            if scoremat[i,j] < 0
                scoremat[i,j] = 0
            end
        end
    end
    return scoremat
end


text:
"""
A conditional linear Gaussian CPD, always returns a Normal{Float64}

    This is a combination of the CategoricalCPD and the LinearGaussianCPD.
    For a variable with N discrete parents and M continuous parents, it will construct
    a linear gaussian distribution for all M parents for each discrete instantiation.

                      { Normal(μ=a₁×continuous_parents(x) + b₁, σ₁) for discrete instantiation 1
	P(x|parents(x)) = { Normal(μ=a₂×continuous_parents(x) + b₂, σ₂) for discrete instantiation 2
                      { ...
"""
struct ConditionalLinearGaussianCPD <: CPD{Normal}
    target::NodeName
    parents::NodeNames # list of all parents

    parents_disc::NodeNames # list of discrete parents
    parental_ncategories::Vector{Int} # list of instantiation counts for each discrete parent
    linear_gaussians::Vector{LinearGaussianCPD} # set of linear gaussian CPDs over the continuous parents
end

name(cpd::ConditionalLinearGaussianCPD) = cpd.target
parents(cpd::ConditionalLinearGaussianCPD) = cpd.parents
nparams(cpd::ConditionalLinearGaussianCPD) = sum(d->nparams(d), cpd.linear_gaussians)
function (cpd::ConditionalLinearGaussianCPD)(a::Assignment)

    idx = 1
    if !isempty(cpd.parents_disc)

        # get the index in cpd.distributions
        N = length(cpd.parents_disc)
        idx = a[cpd.parents_disc[N]] - 1
        for i in N-1:-1:1
            idx = (a[cpd.parents_disc[i]] - 1 + cpd.parental_ncategories[i]*idx)
        end
        idx += 1
    end

    lingaussian = cpd.linear_gaussians[idx]
    lingaussian(a)
end
(cpd::ConditionalLinearGaussianCPD)() = (cpd)(Assignment()) # cpd()
(cpd::ConditionalLinearGaussianCPD)(pair::Pair{NodeName}...) = (cpd)(Assignment(pair)) # cpd(:A=>1)

function Distributions.fit(::Type{ConditionalLinearGaussianCPD},
    data::DataFrame,
    target::NodeName;
    min_stdev::Float64=0.0, # an optional minimum on the standard deviation
    )

    # no parents

    arr = data[!,target]
    eltype(arr) <: Real || error("fit ConditionalLinearGaussianCPD requrires target to be numeric")

    lingaussian = fit(LinearGaussianCPD, data, target, min_stdev=min_stdev)

    ConditionalLinearGaussianCPD(target, NodeName[], NodeName[], Int[], [lingaussian])
end
function Distributions.fit(::Type{ConditionalLinearGaussianCPD},
    data::DataFrame,
    target::NodeName,
    parents::NodeNames;
    min_stdev::Float64=0.0, # an optional minimum on the standard deviation
    )

    if isempty(parents)
        return fit(ConditionalLinearGaussianCPD, data, target, min_stdev=min_stdev)
    end

    # ---------------------
    # identify discrete and continuous parents

    parents_disc = filter(p->eltype(data[!,p]) <: Int, parents)
    parents_cont = filter(p->eltype(data[!,p]) <: AbstractFloat, parents)

   # ---------------------
    # pull discrete dataset
    # 1st row is all of the data for the 1st parent
    # 2nd row is all of the data for the 2nd parent, etc.
    # calc parent_instantiation_counts

    nparents_disc = length(parents_disc)

    if nparents_disc != 0

        parental_ncategories = Array{Int}(undef, nparents_disc)
        dims = Array{UnitRange{Int64}}(undef, nparents_disc)
        for (i,p) in enumerate(parents_disc)
            parental_ncategories[i] = infer_number_of_instantiations(data[!,p])
            dims[i] = 1:parental_ncategories[i]
        end

        # ---------------------
        # fit linear gaussians

        linear_gaussians = Array{LinearGaussianCPD}(undef, prod(parental_ncategories))
        for (q, parent_instantiation) in enumerate(Iterators.product(dims...))
            indeces = Int[]
            for i in 1 : nrow(data)
                if all(j->data[i,parents_disc[j]]==parent_instantiation[j], 1:nparents_disc) # parental instantiation matches
                    push!(indeces, i)
                end
            end
            linear_gaussians[q] = fit(LinearGaussianCPD, data[indeces, :], target, parents_cont, min_stdev=min_stdev)
        end
        ConditionalLinearGaussianCPD(target, parents, parents_disc, parental_ncategories, linear_gaussians)

    else # no discrete parents
        lingaussian = fit(LinearGaussianCPD, data, target, parents, min_stdev=min_stdev)
        ConditionalLinearGaussianCPD(target, parents, NodeName[], Int[], [lingaussian])
    end
end


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

using Cxx

cxxinclude("TH1D.h")
cxxinclude("TH1F.h")
cxxinclude("TH1I.h")
cxxinclude("TH2D.h")
cxxinclude("TH2F.h")
cxxinclude("TH2I.h")
cxxinclude("TH3D.h")
cxxinclude("TH3F.h")
cxxinclude("TH3I.h")

import Base: show, display, size, push!
import StatsBase: fit

export edge, edges
export getname
export gettitle
export draw

export TH1I, TH1IPtr, TH1IInst, new_TH1I
const TH1I = cxxt"TH1I"
const TH1IPtr = pcpp"TH1I"
const TH1IInst = Union{TH1I, TH1IPtr}

export TH1F, TH1FPtr, TH1FInst, new_TH1F
const TH1F = cxxt"TH1F"
const TH1FPtr = pcpp"TH1F"
const TH1FInst = Union{TH1F, TH1FPtr}

export TH1D, TH1DPtr, TH1DInst, new_TH1D
const TH1D = cxxt"TH1D"
const TH1DPtr = pcpp"TH1D"
const TH1DInst = Union{TH1D, TH1DPtr}

export AnyTH1, AnyTH1Ptr, AnyTH1Inst
const AnyTH1 = Union{TH1I, TH1F, TH1D}
const AnyTH1Ptr = Union{TH1IPtr, TH1FPtr, TH1DPtr}
const AnyTH1Inst = Union{AnyTH1, AnyTH1Ptr}

export TH2I, TH2IPtr, TH2IInst, new_TH2I
const TH2I = cxxt"TH2I"
const TH2IPtr = pcpp"TH2I"
const TH2IInst = Union{TH2I, TH2IPtr}

export TH2F, TH2FPtr, TH2FInst, new_TH2F
const TH2F = cxxt"TH2F"
const TH2FPtr = pcpp"TH2F"
const TH2FInst = Union{TH2F, TH2FPtr}

export TH2D, TH2DPtr, TH2DInst, new_TH2D
const TH2D = cxxt"TH2D"
const TH2DPtr = pcpp"TH2D"
const TH2DInst = Union{TH2D, TH2DPtr}

export AnyTH2, AnyTH2Ptr, AnyTH2Inst
const AnyTH2 = Union{TH2I, TH2F, TH2D}
const AnyTH2Ptr = Union{TH2IPtr, TH2FPtr, TH2DPtr}
const AnyTH2Inst = Union{AnyTH2, AnyTH2Ptr}

export TH3I, TH3IPtr, TH3IInst, new_TH3I
const TH3I = cxxt"TH3I"
const TH3IPtr = pcpp"TH3I"
const TH3IInst = Union{TH3I, TH3IPtr}

export TH3F, TH3FPtr, TH3FInst, new_TH3F
const TH3F = cxxt"TH3F"
const TH3FPtr = pcpp"TH3F"
const TH3FInst = Union{TH3F, TH3FPtr}

export TH3D, TH3DPtr, TH3DInst, new_TH3D
const TH3D = cxxt"TH3D"
const TH3DPtr = pcpp"TH3D"
const TH3DInst = Union{TH3D, TH3DPtr}

export AnyTH3, AnyTH3Ptr, AnyTH3Inst
const AnyTH3 = Union{TH3I, TH3F, TH3D}
const AnyTH3Ptr = Union{TH3IPtr, TH3FPtr, TH3DPtr}
const AnyTH3Inst = Union{AnyTH3, AnyTH3Ptr}

const AnyTH = Union{AnyTH1, AnyTH2, AnyTH3}
const AnyTHPtr = Union{AnyTH1Ptr, AnyTH2Ptr, AnyTH3Ptr}
const AnyTHInst = Union{AnyTH1Inst, AnyTH2Inst, AnyTH3Inst}

export THxx, new_THxx

init_thxx(hist::AnyTHInst, edges::NTuple{N,AbstractRange}) where N = begin
    @cxx hist->SetBins(axis_spec(edges)...)
    hist
end

init_thxx(hist::AnyTHInst, edges::NTuple{N, AbstractRange}, name::AbstractString, title::AbstractString) where N = begin
    init_thxx(hist, edges)
    @cxx hist->SetNameTitle(pointer(name), pointer(title))
    hist
end

size(hist::AnyTH1Inst) = ((@cxx hist->GetNbinsX()),)
size(hist::AnyTH2Inst) = ((@cxx hist->GetNbinsX()), (@cxx hist->GetNbinsY()))
size(hist::AnyTH3Inst) = ((@cxx hist->GetNbinsX()), (@cxx hist->GetNbinsY()), (@cxx hist->GetNbinsZ()))

xaxis(hist::AnyTHInst) = @cxx hist->GetXaxis()
yaxis(hist::Union{AnyTH2Inst, AnyTH3Inst}) = @cxx hist->GetYaxis()
zaxis(hist::AnyTH3Inst) = @cxx hist->GetZaxis()

edge(hist::AnyTH1Inst) = range(xaxis(hist))

edges(hist::AnyTH1Inst) = (range(xaxis(hist)),)
edges(hist::AnyTH2Inst) = range(xaxis(hist)), range(yaxis(hist))
edges(hist::AnyTH3Inst) = (range(xaxis(hist)), range(yaxis(hist)), range(zaxis(hist)))

THxx(::Type{Int32}, edges::NTuple{1, AbstractRange}, args...) = init_thxx((@cxx TH1I()), edges, args...)
new_THxx(::Type{Int32}, edges::NTuple{1, AbstractRange}, args...) = init_thxx(icxx""" new TH1I; """, edges, args...)

THxx(::Type{Float32}, edges::NTuple{1, AbstractRange}, args...) = init_thxx((@cxx TH1F()), edges, args...)
new_THxx(::Type{Float32}, edges::NTuple{1, AbstractRange}, args...) = init_thxx(icxx""" new TH1F; """, edges, args...)

THxx(::Type{Float64}, edges::NTuple{1, AbstractRange}, args...) = init_thxx((@cxx TH1D()), edges, args...)
new_THxx(::Type{Float64}, edges::NTuple{1, AbstractRange}, args...) = init_thxx(icxx""" new TH1D; """, edges, args...)

THxx(::Type{T}, edge::AbstractRange, args...) where {T<:Real} = THxx(T, (edge,), args...)
new_THxx(::Type{T}, edge::AbstractRange, args...) where {T<:Real} = new_THxx(T, (edge,), args...)

THxx(::Type{Int32}, edges::NTuple{2, AbstractRange}, args...) = init_thxx((@cxx TH2I()), edges, args...)
new_THxx(::Type{Int32}, edges::NTuple{2, AbstractRange}, args...) = init_thxx(icxx""" new TH2I; """, edges, args...)

THxx(::Type{Float32}, edges::NTuple{2, AbstractRange}, args...) = init_thxx((@cxx TH2F()), edges, args...)
new_THxx(::Type{Float32}, edges::NTuple{2, AbstractRange}, args...) = init_thxx(icxx""" new TH2F; """, edges, args...)

THxx(::Type{Float64}, edges::NTuple{2, AbstractRange}, args...) = init_thxx((@cxx TH2D()), edges, args...)
new_THxx(::Type{Float64}, edges::NTuple{2, AbstractRange}, args...) = init_thxx(icxx""" new TH2D; """, edges, args...)

THxx(::Type{Int32}, edges::NTuple{3, AbstractRange}, args...) = init_thxx((@cxx TH3I()), edges, args...)
new_THxx(::Type{Int32}, edges::NTuple{3, AbstractRange}, args...) = init_thxx(icxx""" new TH3I; """, edges, args...)

THxx(::Type{Float32}, edges::NTuple{3, AbstractRange}, args...) = init_thxx((@cxx TH3F()), edges, args...)
new_THxx(::Type{Float32}, edges::NTuple{3, AbstractRange}, args...) = init_thxx(icxx""" new TH3F; """, edges, args...)

THxx(::Type{Float64}, edges::NTuple{3, AbstractRange}, args...) = init_thxx((@cxx TH3D()), edges, args...)
new_THxx(::Type{Float64}, edges::NTuple{3, AbstractRange}, args...) = init_thxx(icxx""" new TH3D; """, edges, args...)

THxx(edges::NTuple{N, AbstractRange}, args...) where N = THxx(Float64, edges, args...)
new_THxx(edges::NTuple{N, AbstractRange}, args...) where N = new_THxx(Float64, edges, args...)

THxx(edge::AbstractRange, args...) = THxx(Float64, (edge, ), args...)
new_THxx(edge::AbstractRange, args...) = new_THxx(Float64, (edge, ), args...)

show(io::IO, hist::AnyTH) = begin
    tpname = string(typeof(hist).parameters[1].parameters[1].parameters[1])
    print(io, "$tpname($(size(hist)), \"$(escape_string(getname(hist)))\", \"$(escape_string(gettitle(hist)))\")")
end

getname(obj::AnyTHInst) = unsafe_string(@cxx obj->GetName())

gettitle(obj::AnyTHInst) = unsafe_string(@cxx obj->GetTitle())

push!(hist::AnyTH1Inst, x::T) where {T <: Union{Integer, AbstractFloat}} = @cxx hist->Fill(x)
push!(hist::AnyTH1Inst, x::NTuple{1, T}) where {T <: Union{Integer, AbstractFloat}} = @cxx hist->Fill(x[1])
push!(hist::AnyTH2Inst, x::NTuple{2, T}) where {T <: Union{Integer, AbstractFloat}} = @cxx hist->Fill(x[1], x[2])
push!(hist::AnyTH3Inst, x::NTuple{3, T}) where {T <: Union{Integer, AbstractFloat}} = @cxx hist->Fill(x[1], x[2], x[3])

draw(hist::AnyTHInst) = @cxx hist->Draw()
draw(hist::AnyTHInst, option) = @cxx hist->Draw(pointer(option))

fit(hist::AnyTHInst, f::AbstractString) = @cxx hist->Fit(pointer(f))


text:
###########################
# Ladder Operator Methods #
###########################
    # n specifies which particle (a.k.a tensor product
    # factor) the operator acts on
    raiseop(b::AbstractBasis, n=1) = QuArray(raisematrix(size(b), n), b, b)
    raiseop(lens, n=1) = raiseop(FiniteBasis(lens), n)

    lowerop(b::AbstractBasis, n=1) = QuArray(lowermatrix(size(b), n), b, b)
    lowerop(lens, n=1) = lowerop(FiniteBasis(lens), n)

    ##########################
    # Helper Functions/Types #
    ##########################
    raisematrix(lens, particle) = eye_sandwich(lens, particle, raise_single(lens[particle]))
    lowermatrix(lens, particle) = eye_sandwich(lens, particle, lower_single(lens[particle]))

    laddercoeffs(n) = sqrt(linspace(1, n, n))
    lower_single(n) = sparse([1:n-1;], [2:n;], laddercoeffs(n-1), n, n)
    raise_single(n) = sparse([2:n;], [1:n-1;], laddercoeffs(n-1), n, n)

    before_eye(lens, pivot) = speye(prod(lens[1:pivot-1]))
    after_eye(lens, pivot) = speye(prod(lens[pivot+1:length(lens)]))
    function eye_sandwich(lens, pivot, op)
        return kron(before_eye(lens, pivot),
                    op,
                    after_eye(lens, pivot))
    end

function positionop(n::Int)
    cop = raiseop(n)
    return scale!(1/sqrt(2.),cop+cop')
end

function momentumop(n::Int)
    cop = raiseop(n)
    return scale(im/sqrt(2.), cop-cop')
end

squeeze_construct(a::AbstractQuMatrix,b::AbstractQuMatrix, z::Number) = scale!(0.5,(scale!(z',a*b)-scale!(z,a'*b')))
squeezingop(a::AbstractQuMatrix, b::AbstractQuMatrix, z::Number) = expm(squeeze_construct(a,b,z))
squeezingop(n::Int,z::Number) = squeezingop(lowerop(n),lowerop(n),z)

displaceop(n::Int,alpha::Number) = expm(scale(alpha,lowerop(n)')-scale(alpha',lowerop(n)))

export raiseop,
    lowerop,
    positionop,
    momentumop,
    squeezingop,
    displaceop


text:
getindex(a::Number, ::Block{0}) = a

function _sym_axes(A)
    ax = axes(parent(A),2)
    (ax, ax)
end

# Symmetric and Triangular should inherit blocks from parent
axes(A::HermOrSym{<:Any,<:AbstractBlockMatrix}) = _sym_axes(A)
axes(A::HermOrSym{<:Any,<:SubArray{<:Any,2,<:AbstractBlockMatrix}}) = _sym_axes(A)
axes(A::AbstractTriangular{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))
axes(A::AbstractTriangular{<:Any,<:SubArray{<:Any,2,<:AbstractBlockMatrix}}) = axes(parent(A))
function axes(D::Diagonal{<:Any,<:AbstractBlockVector})
    a = axes(parent(D),1)
    (a,a)
end
ArrayLayouts.diagonallayout(lay::AbstractBlockLayout) = DiagonalLayout{typeof(lay)}()

blocksize(A::AbstractTriangular) = blocksize(parent(A))
blocksize(A::AbstractTriangular, i::Int) = blocksize(parent(A), i)
blockaxes(A::AbstractTriangular) = blockaxes(parent(A))

hasmatchingblocks(A::AbstractTriangular) = hasmatchingblocks(parent(A))
hasmatchingblocks(A::HermOrSym) = true


text:
module ConstraintProgrammingExtensions
  using MathOptInterface

  """
      AllDifferent(dimension::Int)
  The set corresponding to an all-different constraint.
  All expressions of a vector-valued function are enforced to take distinct
  values in the solution: for all pairs of expressions, their values must
  differ.
  This constraint is sometimes called `distinct`.
  ## Example
      [x, y, z] in AllDifferent(3)
      # enforces `x != y` AND `x != z` AND `y != z`.
  """
  struct AllDifferent <: AbstractVectorSet
      dimension::Int
  end

  """
      Domain{T <: Number}(values::Set{T})
  The set corresponding to an enumeration of constant values.
  The value of a scalar function is enforced to take a value from this set of
  values.
  This constraint is sometimes called `in`.
  ## Example
      x in Domain(Set(1, 2, 3))
      # enforces `x == 1` OR `x == 2` OR `x == 3`.
  """
  struct Domain{T <: Number} <: AbstractScalarSet
      values::Set{T}
  end

  function Base.copy(set::Domain{T}) where T
      return Domain(copy(set.values))
  end

  """
      Membership(dimension)
  The first element of a function of dimension `dimension` must equal at least
  one of the following `dimension - 1` elements of the function.
  This constraint is sometimes called `in_set`.
  ## Example
      [x, y, z] in Membership(dimension)
      # enforces `x == y` OR `x == z`.
  """
  struct Membership <: AbstractScalarSet
      values::Set{AbstractScalarFunction}
  end

  function Base.copy(set::Membership)
      return Domain(copy(set.values))
  end

  """
      DifferentFrom(dimension)
  The first element of a function of dimension `dimension` cannot equal any
  of the following `dimension - 1` elements of the function.
  ## Example
      [x, y, z] in DifferentFrom(dimension)
      # enforces `x != y` AND `x != z`.
  """
  struct DifferentFrom <: AbstractScalarSet
      dimension::Int
  end

  """
      Count{T <: Real}(value::T, dimension::Int)
  ``\\{(y, x) \\in \\mathbb{N} \\times \\mathbb{R}^n : y = |\\{i | x_i = value\\}|\\}``
  """
  struct Count{T <: Real} <: AbstractVectorSet
      value::T
      dimension::Int
  end

  dimension(set::Count{T}) where T = set.dimension + 1

  function Base.copy(set::Count)
      return Count(copy(set.value), value)
  end

  """
      CountDistinct(dimension::Int)
  The first variable in the set is forced to be the number of distinct values in the rest of the expressions.
  Also called `nvalues`. This is a relaxed version of `AllDifferent`; it encodes an `AllDifferent` constraint
  when the first variable is the number of variables in the set.
  """
  struct CountDistinct <: AbstractVectorSet
      dimension::Int
  end

  dimension(set::CountDistinct) where T = set.dimension + 1

  """
      Strictly{S <: Union{LessThan, GreaterThan}}
  Converts an inequality set to a set with the same inequality made strict. For instance, while LessThan(1)
  corresponds to the inequality `<= 1`, Strictly(LessThan(1)) corresponds to the inequality `< 1`, i.e.
  the value 1 is no more allowed.
  """
  struct Strictly{S <: Union{LessThan, GreaterThan}} <: AbstractScalarSet
      set::S
  end

  function Base.copy(set::Strictly{S}) where S
      return Count(copy(set.set))
  end

  """
      Element{T <: Real}(values::Vector{T})
  ``\\{(x, i) \\in \\mathbb{R}^d \\times \\mathbb{N}^d | x_j = values[i_j]\\}``
  Less formally, the j-th element constrained in this set will take the value of `values` at the index given by
  the (j + dimension)-th element.
  ## Examples
      [x, 3] in Element([4, 5, 6])
      # Enforces that x = 6, because 6 is the 3rd element from the array.
      [x, 3, y, j] in Element([4, 5, 6])
      # Enforces that x = 6, because 6 is the 3rd element from the array.
      # Enforces that y = array[j], depending on the value of j (an integer between 1 and 3).
  """
  struct Element{T <: Real} <: AbstractVectorSet
      values::Vector{T}
      dimension::Int
  end

  dimension(set::Element{T}) where T = 2 * set.dimension

  function Base.copy(set::Element{T}) where T
      return Element(copy(set.values), set.dimension)
  end

  """
      Sort(dimension::Int)
  Ensures that the first `dimension` elements is a sorted copy of the next `dimension` elements.
  ## Example
      [a, b, c, d] in Sort(2)
      # Enforces that:
      # - the first part is sorted: a <= b
      # - the first part corresponds to the second one:
      #     - either a = c and b = d
      #     - or a = d and b = c
  """
  struct Sort <: AbstractVectorSet
      dimension::Int
  end

  dimension(set::Sort) where T = 2 * set.dimension

  """
      SortPermutation(dimension::Int)
  Ensures that the first `dimension` elements is a sorted copy of the next `dimension` elements.
  The last `dimension` elemenst give a permutation to get from the original array to its sorted version.
  ## Example
      [a, b, c, d, i, j] in SortPermutation(2)
      # Enforces that:
      # - the first part is sorted: a <= b
      # - the first part corresponds to the second one:
      #     - either a = c and b = d: in this case, i = 1 and j = 2
      #     - or a = d and b = c: in this case, i = 2 and j = 1
  """
  struct SortPermutation <: AbstractVectorSet
      dimension::Int
  end

  dimension(set::SortPermutation) where T = 3 * set.dimension

  """
      BinPacking(n_bins::Int, n_items::Int)
  Implements an uncapacitated version of the bin packing. The first `n_bins` variables give the load in each bin,
  the next `n_items` give the number of the bin to which the item is assigned to, and the last `n_items` ones
  give the size of each item. The load of a bin is defined as the sum of the sizes of the items put in that bin.
  Also called `pack`.
  ## Example
      [a, b, c, d, e] in BinPacking(1, 2)
      # As there is only one bin, the only solution is to put all the items in that bin.
      # Enforces that:
      # - the bin load is the sum of the objects in that bin: a = d + e
      # - the bin number of the two items is 1: b = c = 1
  """
  struct BinPacking <: AbstractVectorSet
      n_bins::Int
      n_items::Int
  end

  dimension(set::BinPacking) where T = set.n_bins + 2 * set.n_items

  """
      CapacitatedBinPacking(n_bins::Int, n_items::Int)
  Implements an uncapacitated version of the bin packing. The first `n_bins` variables give the load in each bin,
  the next `n_bins` are the capacity of each bin, the followin `n_items` give the number of the bin to which the item
  is assigned to, and the last `n_items` ones give the size of each item. The load of a bin is defined as the sum
  of the sizes of the items put in that bin.
  This constraint is equivalent to `BinPacking` with inequality constraints on the loads of the bins. However, there
  are more efficient propagators for the combined constraint (bin packing with maximum load).
  Also called `bin_packing_capa`.
  ## Example
      [a, 2, b, c, d, e] in CapacitatedBinPacking(1, 2)
      # As there is only one bin, the only solution is to put all the items in that bin if its capacity is large enough.
      # Enforces that:
      # - the bin load is the sum of the objects in that bin: a = d + e
      # - the bin load is at most its capacity: a <= 2
      # - the bin number of the two items is 1: b = c = 1
  """
  struct CapacitatedBinPacking <: AbstractVectorSet
      n_bins::Int
      n_items::Int
  end

  dimension(set::CapacitatedBinPacking) where T = 2 * set.n_bins + 2 * set.n_items

  """
      ReificationSet{S <: AbstractScalarSet}(set::S)
  ``\\{(y, x) \\in \\{0, 1\\} \\times \\mathbb{R}^n | y = 1 \\iff x \\in set, y = 0 otherwise\\}``.
  `S` has to be a sub-type of `AbstractScalarSet`.
  This set serves to find out whether a given constraint is satisfied. The only possible values are
  0 and 1.
  """
  struct ReificationSet{S <: AbstractScalarSet} <: AbstractVectorSet
      set::S
      dimension::Int
  end

  dimension(set::ReificationSet{T}) where T = set.dimension + 1

  function Base.copy(set::ReificationSet{S}) where S
      return ReificationSet(copy(set.set), set.dimension)
  end

  # isbits types, nothing to copy
  function Base.copy(
      set::Union{AllDifferent, DifferentFrom, CountDistinct, Element,
          BinPacking, CapacitatedBinPacking, DifferentFrom}
  )
      return set
  end
end


text:
# This is the canonical kernel: 1D Gaussian.
function evalkernel(x, z::T, θ::GaussianKernelODEType)::T where T <: Real

    # a is ϵ^2 in "Stable eval of Gaussian radial basis func interpolants", 2012.
    out::T = θ.s[1]*exp(-θ.ϵ_sq[1]*(x-z)^2)

    @assert isfinite(out)
    return out
end

"""
evalDkernel(x, z::T, θ)::T

evaluates ∂k(a,b)_∂a at (a,b) = (x,z). See Mathematica script for derivation.
"""
function evalDkernel(x, z::T, θ::GaussianKernelODEType)::T where T <: Real

    τ::T = x-z
    out::T = -2*θ.ϵ_sq[1]*τ*exp(-θ.ϵ_sq[1]*τ^2)
    out = θ.s[1]*out

    @assert isfinite(out)
    return out
end

"""
evalkernel𝐷(x, z::T, θ)::T

evaluates ∂k(a,b)_∂b at (a,b) = (x,z).
See Mathematica script for derivation.

Test:

import RKHSRegularization
using LinearAlgebra
θ = RKHSRegularization.GaussianKernelODEType(1.23, 0.321)
s = randn()
z = randn()
p1 = RKHSRegularization.evalDkernel(s,z,θ)
p2 = RKHSRegularization.evalkernel𝐷(z,s,θ)
println("discrepancy = ", norm(p1-p2))

"""
function evalkernel𝐷(x, z::T, θ::GaussianKernelODEType)::T where T <: Real

    τ::T = x-z
    out::T = 2*θ.ϵ_sq[1]*τ*exp(-θ.ϵ_sq[1]*τ^2)
    out = θ.s[1]*out

    @assert isfinite(out)
    return out
end

"""
evalDkernel𝐷(x, z::T, θ)::T

evaluates ∂k(a,b)_∂a_∂b at (a,b) = (x,z). See Mathematica script for derivation.
"""
function evalDkernel𝐷(x, z::T, θ::GaussianKernelODEType)::T where T <: Real

    τ_sq::T = (x-z)^2
    out::T = 2*θ.ϵ_sq[1] *exp(-θ.ϵ_sq[1] *τ_sq) *(1 - 2 *θ.ϵ_sq[1] *τ_sq)
    out = θ.s[1]*out

    @assert isfinite(out)
    return out
end

##########

max_stars_repo_name: RoyCCWang/RKHSRegularization

max_issues_repo_name: RoyCCWang/RKHSRegularization

max_forks_repo_name: RoyCCWang/RKHSRegularization


text:
# Extra constants, expanding PhysicalConstants
@constant(
    MuonMass, m_μ, "Muon mass",
    1.883_531_627e-28,
    BigFloat(1_883_531_627)/BigFloat(10_000_000_000_000_000_000_000_000_000_000_000_000),
    kg,
    4.2e-36,
    BigFloat(42)/BigFloat(10_000_000_000_000_000_000_000_000_000_000_000_000),
    "CODATA 2018"
)

const muon_mass_au = AtomicMiscellany.MuonMass / CODATA2018.ElectronMass

const α = real(CODATA2018.FineStructureConstant)


text:
function rotate_char(num::Int64, char::Char)
    if char in 'a':'z'
        return 'a' + (char - 'a' + num) % 26
    elseif char in 'A':'Z'
        return 'A' + (char - 'A' + num) % 26
    else
        return char
    end
end

function rotate(num::Int64, text::String)
    map(char->rotate_char(num, char), text)
end

rotate(num::Int64, char::Char)  = rotate_char(num, char)

for i=0:26
    @eval macro $(Symbol('R', i, "_str"))(text::String)
        rotate($i, text)
    end
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism