text:
abstract type AbstractMPS end

"""
    length(::MPS/MPO)

The number of sites of an MPS/MPO.
"""
length(m::AbstractMPS) = length(data(m))

"""
    size(::MPS/MPO)

The number of sites of an MPS/MPO, in a tuple.
"""
size(m::AbstractMPS) = size(data(m))

ndims(m::AbstractMPS) = ndims(data(m))

"""
    promote_itensor_eltype(m::MPS)
    promote_itensor_eltype(m::MPO)

Return the promotion of the elements type of the
tensors in the MPS or MPO. For example,
if all tensors have type `Float64` then
return `Float64`. But if one or more tensors
have type `ComplexF64`, return `ComplexF64`.
"""
promote_itensor_eltype(m::AbstractMPS) = promote_itensor_eltype(data(m))

function promote_itensor_eltype(m::Vector{ITensor})
  T = isassigned(m, 1) ? eltype(m[1]) : Number
  for n in 2:length(m)
    Tn = isassigned(m, n) ? eltype(m[n]) : Number
    T = promote_type(T, Tn)
  end
  return T
end

"""
    eltype(m::MPS)
    eltype(m::MPO)

The element type of the MPS/MPO. Always returns `ITensor`.

For the element type of the ITensors of the MPS/MPO,
use `promote_itensor_eltype`.
"""
eltype(::AbstractMPS) = ITensor

"""
    ITensors.data(::MPS/MPO)

Returns a view of the Vector storage of an MPS/MPO.

This is not exported and mostly for internal usage, please let us
know if there is functionality not available for MPS/MPO you would like.
"""
data(m::AbstractMPS) = m.data

leftlim(m::AbstractMPS) = m.llim

rightlim(m::AbstractMPS) = m.rlim

function setleftlim!(m::AbstractMPS, new_ll::Integer)
  return m.llim = new_ll
end

function setrightlim!(m::AbstractMPS, new_rl::Integer)
  return m.rlim = new_rl
end

"""
    ortho_lims(::MPS/MPO)

Returns the range of sites of the orthogonality center of the MPS/MPO.

# Examples

```julia
s = siteinds("S=½", 5)
ψ = randomMPS(s)
orthogonalize!(ψ, 3)

# ortho_lims(ψ) = 3:3
@show ortho_lims(ψ)

ψ[2] = randomITensor(inds(ψ[2]))

# ortho_lims(ψ) = 2:3
@show ortho_lims(ψ)
```
"""
function ortho_lims(ψ::AbstractMPS)
  return (leftlim(ψ) + 1):(rightlim(ψ) - 1)
end

"""
    ITensors.set_ortho_lims!(::MPS/MPO, r::UnitRange{Int})

Sets the range of sites of the orthogonality center of the MPS/MPO.

This is not exported and is an advanced feature that should be used with
care. Setting the orthogonality limits wrong can lead to incorrect results
when using ITensor MPS/MPO functions.

If you are modifying an MPS/MPO and want the orthogonality limits to be
preserved, please see the `@preserve_ortho` macro.
"""
function set_ortho_lims!(ψ::AbstractMPS, r::UnitRange{Int})
  setleftlim!(ψ, first(r) - 1)
  setrightlim!(ψ, last(r) + 1)
  return ψ
end

reset_ortho_lims!(ψ::AbstractMPS) = set_ortho_lims!(ψ, 1:length(ψ))

isortho(m::AbstractMPS) = leftlim(m) + 1 == rightlim(m) - 1

# Could also define as `only(ortho_lims)`
function orthocenter(m::AbstractMPS)
  !isortho(m) && error("$(typeof(m)) has no well-defined orthogonality center")
  return leftlim(m) + 1
end

getindex(M::AbstractMPS, n) = getindex(data(M), n)

isassigned(M::AbstractMPS, n) = isassigned(data(M), n)

lastindex(M::AbstractMPS) = lastindex(data(M))

"""
    @preserve_ortho

Specify that a block of code preserves the orthogonality limits of
an MPS/MPO that is being modified in-place. The first input is either
a single MPS/MPO or a tuple of the MPS/MPO whose orthogonality limits
should be preserved.

# Examples

```julia
s = siteinds("S=1/2", 4)

# Make random MPS with bond dimension 2
ψ₁ = randomMPS(s, "↑", 2)
ψ₂ = randomMPS(s, "↑", 2)
orthogonalize!(ψ₁, 1)
orthogonalize!(ψ₂, 1)

# ortho_lims(ψ₁) = 1:1
@show ortho_lims(ψ₁)

# ortho_lims(ψ₂) = 1:1
@show ortho_lims(ψ₂)

@preserve_ortho (ψ₁, ψ₂) begin
  ψ₁ .= addtags.(ψ₁, "x₁"; tags = "Link")
  ψ₂ .= addtags.(ψ₂, "x₂"; tags = "Link")
end

# ortho_lims(ψ₁) = 1:1
@show ortho_lims(ψ₁)

# ortho_lims(ψ₂) = 1:1
@show ortho_lims(ψ₂)
```
"""
macro preserve_ortho(ψ, block)
  quote
    if $(esc(ψ)) isa AbstractMPS
      local ortho_limsψ = ortho_lims($(esc(ψ)))
    else
      local ortho_limsψ = ortho_lims.($(esc(ψ)))
    end
    r = $(esc(block))
    if $(esc(ψ)) isa AbstractMPS
      set_ortho_lims!($(esc(ψ)), ortho_limsψ)
    else
      set_ortho_lims!.($(esc(ψ)), ortho_limsψ)
    end
    r
  end
end

function setindex!(M::AbstractMPS, T::ITensor, n::Integer; set_limits::Bool=true)
  if set_limits
    (n <= leftlim(M)) && setleftlim!(M, n - 1)
    (n >= rightlim(M)) && setrightlim!(M, n + 1)
  end
  data(M)[n] = T
  return M
end

function setindex!(M::MPST, v::MPST, ::Colon) where {MPST<:AbstractMPS}
  setleftlim!(M, leftlim(v))
  setrightlim!(M, rightlim(v))
  data(M)[:] = data(v)
  return M
end

setindex!(M::AbstractMPS, v::Vector{<:ITensor}, ::Colon) = setindex!(M, MPS(v), :)

"""
    copy(::MPS)
    copy(::MPO)

Make a shallow copy of an MPS or MPO. By shallow copy, it means that a new MPS/MPO
is returned, but the data of the tensors are still shared between the returned MPS/MPO
and the original MPS/MPO.

Therefore, replacing an entire tensor of the returned MPS/MPO will not modify the input MPS/MPO,
but modifying the data of the returned MPS/MPO will modify the input MPS/MPO.

Use [`deepcopy`](@ref) for an alternative that copies the ITensors as well.

# Examples
```julia
julia> using ITensors

julia> s = siteinds("S=1/2", 3);

julia> M1 = randomMPS(s; linkdims=3);

julia> norm(M1)
0.9999999999999999

julia> M2 = copy(M1);

julia> M2[1] *= 2;

julia> norm(M1)
0.9999999999999999

julia> norm(M2)
1.9999999999999998

julia> M3 = copy(M1);

julia> M3[1] .*= 3; # Modifies the tensor data

julia> norm(M1)
3.0000000000000004

julia> norm(M3)
3.0000000000000004
```
"""
copy(m::AbstractMPS) = typeof(m)(copy(data(m)), leftlim(m), rightlim(m))

similar(m::AbstractMPS) = typeof(m)(similar(data(m)), 0, length(m))

"""
    deepcopy(::MPS)
    deepcopy(::MPO)

Make a deep copy of an MPS or MPO. By deep copy, it means that a new MPS/MPO
is returned that doesn't share any data with the input MPS/MPO.

Therefore, modifying the resulting MPS/MPO will note modify the original MPS/MPO.

Use [`copy`](@ref) for an alternative that performs a shallow copy that avoids
copying the ITensor data.

# Examples
```julia
julia> using ITensors

julia> s = siteinds("S=1/2", 3);

julia> M1 = randomMPS(s; linkdims=3);

julia> norm(M1)
1.0

julia> M2 = deepcopy(M1);

julia> M2[1] .*= 2; # Modifies the tensor data

julia> norm(M1)
1.0

julia> norm(M2)
2.0

julia> M3 = copy(M1);

julia> M3[1] .*= 3; # Modifies the tensor data

julia> norm(M1)
3.0

julia> norm(M3)
3.0
```
"""
deepcopy(m::AbstractMPS) = typeof(m)(copy.(data(m)), leftlim(m), rightlim(m))

eachindex(m::AbstractMPS) = 1:length(m)

iterate(M::AbstractMPS) = iterate(data(M))

iterate(M::AbstractMPS, state) = iterate(data(M), state)

"""
    linkind(M::MPS, j::Integer)
    linkind(M::MPO, j::Integer)

Get the link or bond Index connecting the
MPS or MPO tensor on site j to site j+1.

If there is no link Index, return `nothing`.
"""
function linkind(M::AbstractMPS, j::Integer)
  N = length(M)
  (j ≥ length(M) || j < 1) && return nothing
  return commonind(M[j], M[j + 1])
end

"""
    linkinds(M::MPS, j::Integer)
    linkinds(M::MPO, j::Integer)

Get all of the link or bond Indices connecting the
MPS or MPO tensor on site j to site j+1.
"""
function linkinds(M::AbstractMPS, j::Integer)
  N = length(M)
  (j ≥ length(M) || j < 1) && return IndexSet()
  return commoninds(M[j], M[j + 1])
end

linkinds(ψ::AbstractMPS) = [linkind(ψ, b) for b in 1:(length(ψ) - 1)]

function linkinds(::typeof(all), ψ::AbstractMPS)
  return IndexSet[linkinds(ψ, b) for b in 1:(length(ψ) - 1)]
end

#
# Internal tools for checking for default link tags
#

"""
    ITensors.defaultlinktags(b::Integer)

Default link tags for link index connecting sites `b` to `b+1`.
"""
defaultlinktags(b::Integer) = TagSet("Link,l=$b")

"""
    ITensors.hasdefaultlinktags(ψ::MPS/MPO)

Return true if the MPS/MPO has default link tags.
"""
function hasdefaultlinktags(ψ::AbstractMPS)
  ls = linkinds(all, ψ)
  for (b, lb) in enumerate(ls)
    if length(lb) ≠ 1 || tags(only(lb)) ≠ defaultlinktags(b)
      return false
    end
  end
  return true
end

"""
    ITensors.eachlinkinds(ψ::MPS/MPO)

Return an iterator over each of the sets of link indices of the MPS/MPO.
"""
eachlinkinds(ψ::AbstractMPS) = (linkinds(ψ, n) for n in eachindex(ψ)[1:(end - 1)])

"""
    ITensors.eachsiteinds(ψ::MPS/MPO)

Return an iterator over each of the sets of site indices of the MPS/MPO.
"""
eachsiteinds(ψ::AbstractMPS) = (siteinds(ψ, n) for n in eachindex(ψ))

"""
    ITensors.hasnolinkinds(ψ::MPS/MPO)

Return true if the MPS/MPO has no link indices.
"""
function hasnolinkinds(ψ::AbstractMPS)
  for l in eachlinkinds(ψ)
    if length(l) > 0
      return false
    end
  end
  return true
end

"""
    ITensors.insertlinkinds(ψ::MPS/MPO)

If any link indices are missing, insert default ones.
"""
function insertlinkinds(ψ::AbstractMPS)
  ψ = copy(ψ)
  space = hasqns(ψ) ? [QN() => 1] : 1
  linkind(b::Integer) = Index(space; tags=defaultlinktags(b))
  for b in 1:(length(ψ) - 1)
    if length(linkinds(ψ, b)) == 0
      lb = ITensor(1, linkind(b))
      @preserve_ortho ψ begin
        ψ[b] = ψ[b] * lb
        ψ[b + 1] = ψ[b + 1] * dag(lb)
      end
    end
  end
  return ψ
end

"""
    dense(::MPS/MPO)

Given an MPS (or MPO), return a new MPS (or MPO)
having called `dense` on each ITensor to convert each
tensor to use dense storage and remove any QN or other
sparse structure information, if it is not dense already.
"""
function dense(ψ::AbstractMPS)
  ψ = copy(ψ)
  @preserve_ortho ψ ψ .= dense.(ψ)
  return ψ
end

"""
    siteinds(uniqueinds, A::MPO, B::MPS, j::Integer; kwargs...)
    siteinds(uniqueind, A::MPO, B::MPS, j::Integer; kwargs...)

Get the site index (or indices) of MPO `A` that is unique to `A` (not shared with MPS/MPO `B`).
"""
function siteinds(
  f::Union{typeof(uniqueinds),typeof(uniqueind)},
  A::AbstractMPS,
  B::AbstractMPS,
  j::Integer;
  kwargs...,
)
  N = length(A)
  N == 1 && return f(A[j], B[j]; kwargs...)
  j == 1 && return f(A[j], A[j + 1], B[j]; kwargs...)
  j == N && return f(A[j], A[j - 1], B[j]; kwargs...)
  return f(A[j], A[j - 1], A[j + 1], B[j]; kwargs...)
end

"""
    siteinds(uniqueinds, A::MPO, B::MPS)
    siteinds(uniqueind, A::MPO, B::MPO)

Get the site indices of MPO `A` that are unique to `A` (not shared with MPS/MPO `B`), as a `Vector{<:Index}`.
"""
function siteinds(
  f::Union{typeof(uniqueinds),typeof(uniqueind)}, A::AbstractMPS, B::AbstractMPS; kwargs...
)
  return [siteinds(f, A, B, j; kwargs...) for j in eachindex(A)]
end

"""
    siteinds(commoninds, A::MPO, B::MPS, j::Integer; kwargs...)
    siteinds(commonind, A::MPO, B::MPO, j::Integer; kwargs...)

Get the site index (or indices) of  the `j`th MPO tensor of `A` that is shared with MPS/MPO `B`.
"""
function siteinds(
  f::Union{typeof(commoninds),typeof(commonind)},
  A::AbstractMPS,
  B::AbstractMPS,
  j::Integer;
  kwargs...,
)
  return f(A[j], B[j]; kwargs...)
end

"""
    siteinds(commoninds, A::MPO, B::MPS; kwargs...)
    siteinds(commonind, A::MPO, B::MPO; kwargs...)

Get a vector of the site index (or indices) of MPO `A` that is shared with MPS/MPO `B`.
"""
function siteinds(
  f::Union{typeof(commoninds),typeof(commonind)}, A::AbstractMPS, B::AbstractMPS; kwargs...
)
  return [siteinds(f, A, B, j) for j in eachindex(A)]
end

keys(ψ::AbstractMPS) = keys(data(ψ))

#
# Find sites of an MPS or MPO
#

# TODO: accept a keyword argument sitedict that
# is a dictionary from the site indices to the site.
"""
    findsite(M::Union{MPS, MPO}, is)

Return the first site of the MPS or MPO that has at least one
Index in common with the Index or collection of indices `is`.

To find all sites with common indices with `is`, use the
`findsites` function.

# Examples
```julia
s = siteinds("S=1/2", 5)
ψ = randomMPS(s)
findsite(ψ, s[3]) == 3
findsite(ψ, (s[3], s[4])) == 3

M = MPO(s)
findsite(M, s[4]) == 4
findsite(M, s[4]') == 4
findsite(M, (s[4]', s[4])) == 4
findsite(M, (s[4]', s[3])) == 3
```
"""
findsite(ψ::AbstractMPS, is) = findfirst(hascommoninds(is), ψ)

findsite(ψ::AbstractMPS, s::Index) = findsite(ψ, IndexSet(s))

#
# TODO: Maybe make:
# findall(f::Function, siteindsM::Tuple{typeof(siteinds), ::AbstractMPS})
# findall(siteindsM::Tuple{typeof(siteinds), <:AbstractMPS}, is) =
#   findall(hassameinds(is), siteindsM)
#
"""
    findsites(M::Union{MPS, MPO}, is)

Return the sites of the MPS or MPO that have
indices in common with the collection of site indices
`is`.

# Examples
```julia
s = siteinds("S=1/2", 5)
ψ = randomMPS(s)
findsites(ψ, s[3]) == [3]
findsites(ψ, (s[4], s[1])) == [1, 4]

M = MPO(s)
findsites(M, s[4]) == [4]
findsites(M, s[4]') == [4]
findsites(M, (s[4]', s[4])) == [4]
findsites(M, (s[4]', s[3])) == [3, 4]
```
"""
findsites(ψ::AbstractMPS, is) = findall(hascommoninds(is), ψ)

findsites(ψ::AbstractMPS, s::Index) = findsites(ψ, IndexSet(s))

#
# TODO: Maybe make:
# findfirst(f::Function, siteindsM::Tuple{typeof(siteinds), ::AbstractMPS})
# findfirst(siteindsM::Tuple{typeof(siteinds), <:AbstractMPS}, is) =
#   findfirst(hassameinds(is), siteindsM)
#
"""
    findfirstsiteind(M::MPS, i::Index)
    findfirstsiteind(M::MPO, i::Index)

Return the first site of the MPS or MPO that has the
site index `i`.
"""
findfirstsiteind(ψ::AbstractMPS, s::Index) = findfirst(hasind(s), ψ)

# TODO: depracate in favor of findsite.
"""
    findfirstsiteinds(M::MPS, is)
    findfirstsiteinds(M::MPO, is)

Return the first site of the MPS or MPO that has the
site indices `is`.
"""
findfirstsiteinds(ψ::AbstractMPS, s) = findfirst(hasinds(s), ψ)

"""
    siteind(::typeof(first), M::Union{MPS,MPO}, j::Integer; kwargs...)

Return the first site Index found on the MPS or MPO
(the first Index unique to the `j`th MPS/MPO tensor).

You can choose different filters, like prime level
and tags, with the `kwargs`.
"""
function siteind(::typeof(first), M::AbstractMPS, j::Integer; kwargs...)
  N = length(M)
  (N == 1) && return firstind(M[1]; kwargs...)
  if j == 1
    si = uniqueind(M[j], M[j + 1]; kwargs...)
  elseif j == N
    si = uniqueind(M[j], M[j - 1]; kwargs...)
  else
    si = uniqueind(M[j], M[j - 1], M[j + 1]; kwargs...)
  end
  return si
end

"""
    siteinds(M::Union{MPS, MPO}}, j::Integer; kwargs...)

Return the site Indices found of the MPO or MPO
at the site `j` as an IndexSet.

Optionally filter prime tags and prime levels with
keyword arguments like `plev` and `tags`.
"""
function siteinds(M::AbstractMPS, j::Integer; kwargs...)
  N = length(M)
  (N == 1) && return inds(M[1]; kwargs...)
  if j == 1
    si = uniqueinds(M[j], M[j + 1]; kwargs...)
  elseif j == N
    si = uniqueinds(M[j], M[j - 1]; kwargs...)
  else
    si = uniqueinds(M[j], M[j - 1], M[j + 1]; kwargs...)
  end
  return si
end

function siteinds(::typeof(all), ψ::AbstractMPS, n::Integer; kwargs...)
  return siteinds(ψ, n; kwargs...)
end

function siteinds(::typeof(first), ψ::AbstractMPS; kwargs...)
  return [siteind(first, ψ, j; kwargs...) for j in 1:length(ψ)]
end

function siteinds(::typeof(only), ψ::AbstractMPS; kwargs...)
  return [siteind(only, ψ, j; kwargs...) for j in 1:length(ψ)]
end

function siteinds(::typeof(all), ψ::AbstractMPS; kwargs...)
  return [siteinds(ψ, j; kwargs...) for j in 1:length(ψ)]
end

function replaceinds!(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index})
  for i in eachindex(M)[1:(end - 1)]
    l = linkind(M, i)
    l̃ = l̃s[i]
    if !isnothing(l)
      @preserve_ortho M begin
        M[i] = replaceinds(M[i], l => l̃)
        M[i + 1] = replaceinds(M[i + 1], l => l̃)
      end
    end
  end
  return M
end

function replaceinds(::typeof(linkinds), M::AbstractMPS, l̃s::Vector{<:Index})
  return replaceinds!(linkinds, copy(M), l̃s)
end

# TODO: change kwarg from `set_limits` to `preserve_ortho`
function map!(f::Function, M::AbstractMPS; set_limits::Bool=true)
  for i in eachindex(M)
    M[i, set_limits=set_limits] = f(M[i])
  end
  return M
end

# TODO: change kwarg from `set_limits` to `preserve_ortho`
function map(f::Function, M::AbstractMPS; set_limits::Bool=true)
  return map!(f, copy(M); set_limits=set_limits)
end

for fname in (
  :dag,
  :prime,
  :setprime,
  :noprime,
  :swapprime,
  :replaceprime,
  :addtags,
  :removetags,
  :replacetags,
  :settags,
)
  fname! = Symbol(fname, :!)

  @eval begin
    """
        $($fname)[!](M::MPS, args...; kwargs...)
        $($fname)[!](M::MPO, args...; kwargs...)

    Apply $($fname) to all ITensors of an MPS/MPO, returning a new MPS/MPO.

    The ITensors of the MPS/MPO will be a view of the storage of the original ITensors. Alternatively apply the function in-place.
    """
    function $fname(M::AbstractMPS, args...; set_limits::Bool=false, kwargs...)
      return map(m -> $fname(m, args...; kwargs...), M; set_limits=set_limits)
    end

    function $(fname!)(M::AbstractMPS, args...; set_limits::Bool=false, kwargs...)
      return map!(m -> $fname(m, args...; kwargs...), M; set_limits=set_limits)
    end
  end
end

adjoint(M::AbstractMPS) = prime(M)

function hascommoninds(::typeof(siteinds), A::AbstractMPS, B::AbstractMPS)
  N = length(A)
  for n in 1:N
    !hascommoninds(siteinds(A, n), siteinds(B, n)) && return false
  end
  return true
end

function check_hascommoninds(::typeof(siteinds), A::AbstractMPS, B::AbstractMPS)
  N = length(A)
  if length(B) ≠ N
    throw(
      DimensionMismatch(
        "$(typeof(A)) and $(typeof(B)) have mismatched lengths $N and $(length(B))."
      ),
    )
  end
  for n in 1:N
    !hascommoninds(siteinds(A, n), siteinds(B, n)) && error(
      "$(typeof(A)) A and $(typeof(B)) B must share site indices. On site $n, A has site indices $(siteinds(A, n)) while B has site indices $(siteinds(B, n)).",
    )
  end
  return nothing
end

function map!(f::Function, ::typeof(linkinds), M::AbstractMPS)
  for i in eachindex(M)[1:(end - 1)]
    l = linkinds(M, i)
    if !isempty(l)
      l̃ = f(l)
      @preserve_ortho M begin
        M[i] = replaceinds(M[i], l, l̃)
        M[i + 1] = replaceinds(M[i + 1], l, l̃)
      end
    end
  end
  return M
end

map(f::Function, ::typeof(linkinds), M::AbstractMPS) = map!(f, linkinds, copy(M))

function map!(f::Function, ::typeof(siteinds), M::AbstractMPS)
  for i in eachindex(M)
    s = siteinds(M, i)
    if !isempty(s)
      @preserve_ortho M begin
        M[i] = replaceinds(M[i], s, f(s))
      end
    end
  end
  return M
end

map(f::Function, ::typeof(siteinds), M::AbstractMPS) = map!(f, siteinds, copy(M))

function map!(
  f::Function, ::typeof(siteinds), ::typeof(commoninds), M1::AbstractMPS, M2::AbstractMPS
)
  length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
  for i in eachindex(M1)
    s = siteinds(commoninds, M1, M2, i)
    if !isempty(s)
      s̃ = f(s)
      @preserve_ortho (M1, M2) begin
        M1[i] = replaceinds(M1[i], s .=> s̃)
        M2[i] = replaceinds(M2[i], s .=> s̃)
      end
    end
  end
  return M1, M2
end

function map!(
  f::Function, ::typeof(siteinds), ::typeof(uniqueinds), M1::AbstractMPS, M2::AbstractMPS
)
  length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
  for i in eachindex(M1)
    s = siteinds(uniqueinds, M1, M2, i)
    if !isempty(s)
      @preserve_ortho M1 begin
        M1[i] = replaceinds(M1[i], s .=> f(s))
      end
    end
  end
  return M1
end

function map(
  f::Function,
  ffilter::typeof(siteinds),
  fsubset::Union{typeof(commoninds),typeof(uniqueinds)},
  M1::AbstractMPS,
  M2::AbstractMPS,
)
  return map!(f, ffilter, fsubset, copy(M1), copy(M2))
end

function hassameinds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS)
  length(M1) ≠ length(M2) && return false
  for n in 1:length(M1)
    !hassameinds(siteinds(all, M1, n), siteinds(all, M2, n)) && return false
  end
  return true
end

function hassamenuminds(::typeof(siteinds), M1::AbstractMPS, M2::AbstractMPS)
  length(M1) ≠ length(M2) && return false
  for n in 1:length(M1)
    length(siteinds(M1, n)) ≠ length(siteinds(M2, n)) && return false
  end
  return true
end

for fname in
    (:sim, :prime, :setprime, :noprime, :addtags, :removetags, :replacetags, :settags)
  fname! = Symbol(fname, :!)
  @eval begin
    """
        $($fname)[!](linkinds, M::MPS, args...; kwargs...)
        $($fname)[!](linkinds, M::MPO, args...; kwargs...)

    Apply $($fname) to all link indices of an MPS/MPO, returning a new MPS/MPO.

    The ITensors of the MPS/MPO will be a view of the storage of the original ITensors.
    """
    function $fname(ffilter::typeof(linkinds), M::AbstractMPS, args...; kwargs...)
      return map(i -> $fname(i, args...; kwargs...), ffilter, M)
    end

    function $(fname!)(ffilter::typeof(linkinds), M::AbstractMPS, args...; kwargs...)
      return map!(i -> $fname(i, args...; kwargs...), ffilter, M)
    end

    """
        $($fname)[!](siteinds, M::MPS, args...; kwargs...)
        $($fname)[!](siteinds, M::MPO, args...; kwargs...)

    Apply $($fname) to all site indices of an MPS/MPO, returning a new MPS/MPO.

    The ITensors of the MPS/MPO will be a view of the storage of the original ITensors.
    """
    function $fname(ffilter::typeof(siteinds), M::AbstractMPS, args...; kwargs...)
      return map(i -> $fname(i, args...; kwargs...), ffilter, M)
    end

    function $(fname!)(ffilter::typeof(siteinds), M::AbstractMPS, args...; kwargs...)
      return map!(i -> $fname(i, args...; kwargs...), ffilter, M)
    end

    """
        $($fname)[!](siteinds, commoninds, M1::MPO, M2::MPS, args...; kwargs...)
        $($fname)[!](siteinds, commoninds, M1::MPO, M2::MPO, args...; kwargs...)

    Apply $($fname) to the site indices that are shared by `M1` and `M2`.

    Returns new MPSs/MPOs. The ITensors of the MPSs/MPOs will be a view of the storage of the original ITensors.
    """
    function $fname(
      ffilter::typeof(siteinds),
      fsubset::typeof(commoninds),
      M1::AbstractMPS,
      M2::AbstractMPS,
      args...;
      kwargs...,
    )
      return map(i -> $fname(i, args...; kwargs...), ffilter, fsubset, M1, M2)
    end

    function $(fname!)(
      ffilter::typeof(siteinds),
      fsubset::typeof(commoninds),
      M1::AbstractMPS,
      M2::AbstractMPS,
      args...;
      kwargs...,
    )
      return map!(i -> $fname(i, args...; kwargs...), ffilter, fsubset, M1, M2)
    end

    """
        $($fname)[!](siteinds, uniqueinds, M1::MPO, M2::MPS, args...; kwargs...)

    Apply $($fname) to the site indices of `M1` that are not shared with `M2`. Returns new MPSs/MPOs.

    The ITensors of the MPSs/MPOs will be a view of the storage of the original ITensors.
    """
    function $fname(
      ffilter::typeof(siteinds),
      fsubset::typeof(uniqueinds),
      M1::AbstractMPS,
      M2::AbstractMPS,
      args...;
      kwargs...,
    )
      return map(i -> $fname(i, args...; kwargs...), ffilter, fsubset, M1, M2)
    end

    function $(fname!)(
      ffilter::typeof(siteinds),
      fsubset::typeof(uniqueinds),
      M1::AbstractMPS,
      M2::AbstractMPS,
      args...;
      kwargs...,
    )
      return map!(i -> $fname(i, args...; kwargs...), ffilter, fsubset, M1, M2)
    end
  end
end

"""
    maxlinkdim(M::MPS)
    maxlinkdim(M::MPO)

Get the maximum link dimension of the MPS or MPO.

The minimum this will return is `1`, even if there
are no link indices.
"""
function maxlinkdim(M::AbstractMPS)
  md = 1
  for b in eachindex(M)[1:(end - 1)]
    l = linkind(M, b)
    linkdim = isnothing(l) ? 1 : dim(l)
    md = max(md, linkdim)
  end
  return md
end

"""
    linkdim(M::MPS, j::Integer)
    linkdim(M::MPO, j::Integer)

Get the dimension of the link or bond connecting the
MPS or MPO tensor on site j to site j+1.

If there is no link Index, return `nothing`.
"""
function linkdim(ψ::AbstractMPS, b::Integer)
  l = linkind(ψ, b)
  isnothing(l) && return nothing
  return dim(l)
end

linkdims(ψ::AbstractMPS) = [linkdim(ψ, b) for b in 1:(length(ψ) - 1)]

function _log_or_not_dot(
  M1::MPST, M2::MPST, loginner::Bool; make_inds_match::Bool=true
)::Number where {MPST<:AbstractMPS}
  N = length(M1)
  if length(M2) != N
    throw(DimensionMismatch("inner: mismatched lengths $N and $(length(M2))"))
  end
  M1dag = dag(M1)
  sim!(linkinds, M1dag)
  siteindsM1dag = siteinds(all, M1dag)
  siteindsM2 = siteinds(all, M2)
  if any(n -> length(n) > 1, siteindsM1dag) ||
    any(n -> length(n) > 1, siteindsM2) ||
    !hassamenuminds(siteinds, M1, M2)
    # If the MPS have more than one site Indices on any site or they don't have
    # the same number of site indices on each site, don't try to make the
    # indices match
    if !hassameinds(siteinds, M1, M2)
      n = findfirst(n -> !hassameinds(siteinds(M1, n), siteinds(M2, n)), 1:N)
      error(
        """Calling `dot(ϕ::MPS/MPO, ψ::MPS/MPO)` with multiple site indices per MPS/MPO tensor but the site indices don't match. Even with `make_inds_match = true`, the case of multiple site indices per MPS/MPO is not handled automatically. The sites with unmatched site indices are:

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

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

        Make sure the site indices of your MPO/MPS match. You may need to prime one of the MPS, such as `dot(ϕ', ψ)`.""",
      )
    end
    make_inds_match = false
  end
  if make_inds_match
    replace_siteinds!(M1dag, siteindsM2)
  end
  O = M1dag[1] * M2[1]

  if loginner
    normO = norm(O)
    log_inner_tot = log(normO)
    O ./= normO
  end

  for j in eachindex(M1)[2:end]
    O = (O * M1dag[j]) * M2[j]

    if loginner
      normO = norm(O)
      log_inner_tot += log(normO)
      O ./= normO
    end
  end

  if loginner
    return log_inner_tot
  end

  return O[]
end

"""
    dot(A::MPS, B::MPS; make_inds_match = true)
    inner(A::MPS, B::MPS; make_inds_match = true)

    dot(A::MPO, B::MPO)
    inner(A::MPO, B::MPO)

Compute the inner product `<A|B>`. If `A` and `B` are MPOs, computes the Frobenius inner product.

If `make_inds_match = true`, the function attempts to make
the site indices match before contracting (so for example, the
inputs can have different site indices, as long as they
have the same dimensions or QN blocks).

For now, `make_inds_match` is only supported for MPSs.

See also `logdot`/`loginner`.
"""
function dot(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
  return _log_or_not_dot(M1, M2, false; kwargs...)
end

"""
    logdot(A::MPS, B::MPS; make_inds_match = true)
    loginner(A::MPS, B::MPS; make_inds_match = true)

    logdot(A::MPO, B::MPO)
    loginner(A::MPO, B::MPO)

Compute the logarithm of the inner product `<A|B>`. If `A` and `B` are MPOs, computes the logarithm of the Frobenius inner product.

This is useful for larger MPS/MPO, where in the limit of large numbers of sites the inner product can diverge or approach zero.

If `make_inds_match = true`, the function attempts to make
the site indices match before contracting (so for example, the
inputs can have different site indices, as long as they
have the same dimensions or QN blocks).

For now, `make_inds_match` is only supported for MPSs.
"""
function logdot(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
  return _log_or_not_dot(M1, M2, true; kwargs...)
end

inner(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS} = dot(M1, M2; kwargs...)

function loginner(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
  return logdot(M1, M2; kwargs...)
end

"""
    norm(A::MPS)
    norm(A::MPO)

Compute the norm of the MPS or MPO.

If the MPS or MPO has a well defined orthogonality center, this reduces to the norm of the orthogonality center tensor. Otherwise, it computes the norm with the full inner product of the MPS/MPO with itself.

See also [`lognorm`](@ref).
"""
function norm(M::AbstractMPS)
  if isortho(M)
    return norm(M[orthocenter(M)])
  end
  return sqrt(dot(M, M))
end

"""
    lognorm(A::MPS)
    lognorm(A::MPO)

Compute the logarithm of the norm of the MPS or MPO.

This is useful for larger MPS/MPO that are not gauged, where in the limit of large numbers of sites the norm can diverge or approach zero.

See also [`norm`](@ref), [`logdot`](@ref).
"""
function lognorm(M::AbstractMPS)
  if isortho(M)
    return log(norm(M[orthocenter(M)]))
  end
  return 0.5 * logdot(M, M)
end

"""
    dist(A::MPS, B::MPS)
    dist(A::MPO, B::MPO)

Compute the Euclidean distance between to MPS/MPO. Equivalent to `norm(A - B)` but done more efficiently as:
`sqrt(abs(inner(A, A) + inner(B, B) - 2 * real(inner(A, B))))`.

Note that if the MPS/MPO are not normalized, the normalizations may diverge and this may not be accurate. For those cases, likely it is best to use `norm(A - B)` directly (or `lognorm(A - B)` if you expect the result may be very large).
"""
function dist(A::AbstractMPS, B::AbstractMPS)
  return sqrt(abs(inner(A, A) + inner(B, B) - 2 * real(inner(A, B))))
end

function site_combiners(ψ::AbstractMPS)
  N = length(ψ)
  Cs = Vector{ITensor}(undef, N)
  for n in 1:N
    s = siteinds(all, ψ, n)
    Cs[n] = combiner(s; tags=commontags(s))
  end
  return Cs
end

# The maximum link dimensions when adding MPS/MPO
function _add_maxlinkdims(ψ⃗::AbstractMPS...)
  N = length(ψ⃗[1])
  maxdims = Vector{Int}(undef, N - 1)
  for b in 1:(N - 1)
    maxdims[b] = sum(ψ -> linkdim(ψ, b), ψ⃗)
  end
  return maxdims
end

"""
    +(A::MPS/MPO...; kwargs...)
    add(A::MPS/MPO...; kwargs...)

Add arbitrary numbers of MPS/MPO with each other, with some optional
truncation.

A cutoff of 1e-15 is used by default, and in general users should set their own cutoff for their particular application.

In the future we will give an interface for returning the truncation error.

# Examples

```julia
N = 10

s = siteinds("S=1/2", N; conserve_qns = true)

state = n -> isodd(n) ? "↑" : "↓"
ψ₁ = randomMPS(s, state, 2)
ψ₂ = randomMPS(s, state, 2)
ψ₃ = randomMPS(s, state, 2)

ψ = +(ψ₁, ψ₂; cutoff = 1e-8)

# Can use:
#
# ψ = ψ₁ + ψ₂
#
# but generally you want to set a custom `cutoff` and `maxdim`.

println()
@show inner(ψ, ψ)
@show inner(ψ₁, ψ₂) + inner(ψ₁, ψ₂) + inner(ψ₂, ψ₁) + inner(ψ₂, ψ₂)

# Computes ψ₁ + 2ψ₂
ψ = ψ₁ + 2ψ₂

println()
@show inner(ψ, ψ)
@show inner(ψ₁, ψ₁) + 2 * inner(ψ₁, ψ₂) + 2 * inner(ψ₂, ψ₁) + 4 * inner(ψ₂, ψ₂)

# Computes ψ₁ + 2ψ₂ + ψ₃
ψ = ψ₁ + 2ψ₂ + ψ₃

println()
@show inner(ψ, ψ)
@show inner(ψ₁, ψ₁) + 2 * inner(ψ₁, ψ₂) + inner(ψ₁, ψ₃) +
      2 * inner(ψ₂, ψ₁) + 4 * inner(ψ₂, ψ₂) + 2 * inner(ψ₂, ψ₃) +
      inner(ψ₃, ψ₁) + 2 * inner(ψ₃, ψ₂) + inner(ψ₃, ψ₃)
```
"""
function +(ψ⃗::MPST...; cutoff=1e-15, kwargs...) where {MPST<:AbstractMPS}
  # TODO: Check that the inputs have the same site indices

  Nₘₚₛ = length(ψ⃗)

  @assert all(ψᵢ -> length(ψ⃗[1]) == length(ψᵢ), ψ⃗)

  N = length(ψ⃗[1])

  ψ⃗ = copy.(ψ⃗)

  X⃗ = site_combiners(ψ⃗[1])

  for ψᵢ in ψ⃗
    @preserve_ortho ψᵢ ψᵢ .*= X⃗
  end

  ψ⃗ = convert.(MPS, ψ⃗)

  s = siteinds(ψ⃗[1])

  ψ⃗ = orthogonalize.(ψ⃗, N)

  ψ = MPS(N)

  ρ⃗ₙ = [prime(ψᵢ[N], s[N]) * dag(ψᵢ[N]) for ψᵢ in ψ⃗]
  ρₙ = sum(ρ⃗ₙ)

  # Maximum theoretical link dimensions
  add_maxlinkdims = _add_maxlinkdims(ψ⃗...)

  C⃗ₙ = last.(ψ⃗)
  for n in reverse(2:N)
    Dₙ, Vₙ, spec = eigen(
      ρₙ;
      ishermitian=true,
      tags=tags(linkind(ψ⃗[1], n - 1)),
      cutoff=cutoff,
      maxdim=add_maxlinkdims[n - 1],
      kwargs...,
    )
    lₙ₋₁ = commonind(Dₙ, Vₙ)

    # Update the total state
    ψ[n] = Vₙ

    # Compute the new density matrix
    C⃗ₙ₋₁ = [ψ⃗[i][n - 1] * C⃗ₙ[i] * dag(Vₙ) for i in 1:Nₘₚₛ]
    C⃗ₙ₋₁′ = [prime(Cₙ₋₁, (s[n - 1], lₙ₋₁)) for Cₙ₋₁ in C⃗ₙ₋₁]
    ρ⃗ₙ₋₁ = C⃗ₙ₋₁′ .* dag.(C⃗ₙ₋₁)
    ρₙ₋₁ = sum(ρ⃗ₙ₋₁)

    C⃗ₙ = C⃗ₙ₋₁
    ρₙ = ρₙ₋₁
  end

  ψ[1] = sum(C⃗ₙ)
  ψ .*= dag.(X⃗)

  set_ortho_lims!(ψ, 1:1)

  return convert(MPST, ψ)
end

+(ψ::AbstractMPS) = ψ

add(ψ⃗::AbstractMPS...; kwargs...) = +(ψ⃗...; kwargs...)

-(ψ₁::AbstractMPS, ψ₂::AbstractMPS; kwargs...) = +(ψ₁, -ψ₂; kwargs...)

add(A::T, B::T; kwargs...) where {T<:AbstractMPS} = +(A, B; kwargs...)

"""
    sum(A::Vector{MPS}; kwargs...)

    sum(A::Vector{MPO}; kwargs...)

Add multiple MPS/MPO with each other, with some optional
truncation.
"""
function sum(ψ⃗::Vector{T}; kwargs...) where {T<:AbstractMPS}
  length(ψ⃗) == 0 && return T()
  length(ψ⃗) == 1 && return A[1]
  return +(ψ⃗...; kwargs...)
end

"""
    orthogonalize!(M::MPS, j::Int; kwargs...)
    orthogonalize(M::MPS, j::Int; kwargs...)

    orthogonalize!(M::MPO, j::Int; kwargs...)
    orthogonalize(M::MPO, j::Int; kwargs...)

Move the orthogonality center of the MPS
to site `j`. No observable property of the
MPS will be changed, and no truncation of the
bond indices is performed. Afterward, tensors
`1,2,...,j-1` will be left-orthogonal and tensors
`j+1,j+2,...,N` will be right-orthogonal.

Either modify in-place with `orthogonalize!` or
out-of-place with `orthogonalize`.
"""
function orthogonalize!(M::AbstractMPS, j::Int; kwargs...)
  while leftlim(M) < (j - 1)
    (leftlim(M) < 0) && setleftlim!(M, 0)
    b = leftlim(M) + 1
    linds = uniqueinds(M[b], M[b + 1])
    lb = linkind(M, b)
    if !isnothing(lb)
      ltags = tags(lb)
    else
      ltags = TagSet("Link,l=$b")
    end
    L, R = factorize(M[b], linds; tags=ltags, kwargs...)
    M[b] = L
    M[b + 1] *= R
    setleftlim!(M, b)
    if rightlim(M) < leftlim(M) + 2
      setrightlim!(M, leftlim(M) + 2)
    end
  end

  N = length(M)

  while rightlim(M) > (j + 1)
    (rightlim(M) > (N + 1)) && setrightlim!(M, N + 1)
    b = rightlim(M) - 2
    rinds = uniqueinds(M[b + 1], M[b])
    lb = linkind(M, b)
    if !isnothing(lb)
      ltags = tags(lb)
    else
      ltags = TagSet("Link,l=$b")
    end
    L, R = factorize(M[b + 1], rinds; tags=ltags, kwargs...)
    M[b + 1] = L
    M[b] *= R

    setrightlim!(M, b + 1)
    if leftlim(M) > rightlim(M) - 2
      setleftlim!(M, rightlim(M) - 2)
    end
  end
  return M
end

function orthogonalize(ψ0::AbstractMPS, args...; kwargs...)
  ψ = copy(ψ0)
  orthogonalize!(ψ, args...; kwargs...)
  return ψ
end

"""
    truncate!(M::MPS; kwargs...)
    truncate!(M::MPO; kwargs...)

Perform a truncation of all bonds of an MPS/MPO,
using the truncation parameters (cutoff,maxdim, etc.)
provided as keyword arguments.
"""
function truncate!(M::AbstractMPS; kwargs...)
  N = length(M)

  # Left-orthogonalize all tensors to make
  # truncations controlled
  orthogonalize!(M, N)

  # Perform truncations in a right-to-left sweep
  for j in reverse(2:N)
    rinds = uniqueinds(M[j], M[j - 1])
    ltags = tags(commonind(M[j], M[j - 1]))
    U, S, V = svd(M[j], rinds; lefttags=ltags, kwargs...)
    M[j] = U
    M[j - 1] *= (S * V)
    setrightlim!(M, j)
  end
  return M
end

function truncate(ψ0::AbstractMPS; kwargs...)
  ψ = copy(ψ0)
  truncate!(ψ; kwargs...)
  return ψ
end

# Make `*` and alias for `contract` of two `AbstractMPS`
*(A::AbstractMPS, B::AbstractMPS; kwargs...) = contract(A, B; kwargs...)

"""
    α::Number * ψ::MPS/MPO

Scales the MPS or MPO by the provided number.

Currently, this works by scaling one of the sites within the orthogonality limits.
"""
function (α::Number * ψ::AbstractMPS)
  limsψ = ortho_lims(ψ)
  n = first(limsψ)
  αψ = copy(ψ)
  αψ[n] = α * ψ[n]
  return αψ
end

(ψ::AbstractMPS * α::Number) = α * ψ

-(ψ::AbstractMPS) = -1 * ψ

"""
    setindex!(::Union{MPS, MPO}, ::Union{MPS, MPO},
              r::UnitRange{Int64})

Sets a contiguous range of MPS/MPO tensors
"""
function setindex!(ψ::MPST, ϕ::MPST, r::UnitRange{Int64}) where {MPST<:AbstractMPS}
  @assert length(r) == length(ϕ)
  # TODO: accept r::Union{AbstractRange{Int}, Vector{Int}}
  # if r isa AbstractRange
  #   @assert step(r) = 1
  # else
  #   all(==(1), diff(r))
  # end
  llim = leftlim(ψ)
  rlim = rightlim(ψ)
  for (j, n) in enumerate(r)
    ψ[n] = ϕ[j]
  end
  if llim + 1 ≥ r[1]
    setleftlim!(ψ, leftlim(ϕ) + r[1] - 1)
  end
  if rlim - 1 ≤ r[end]
    setrightlim!(ψ, rightlim(ϕ) + r[1] - 1)
  end
  return ψ
end

_isodd_fermionic_parity(s::Index, ::Integer) = false

function _isodd_fermionic_parity(s::QNIndex, n::Integer)
  qn_n = qn(space(s)[n])
  fermionic_qn_pos = findfirst(q -> isfermionic(q), qn_n)
  isnothing(fermionic_qn_pos) && return false
  return isodd(val(qn_n[fermionic_qn_pos]))
end

function _fermionic_swap(s1::Index, s2::Index)
  T = ITensor(QN(), s1', s2', dag(s1), dag(s2))
  for b in nzblocks(T)
    dval = 1.0
    # Must be a diagonal block
    ((b[1] ≠ b[3]) || (b[2] ≠ b[4])) && continue
    n1, n2 = b[1], b[2]
    if _isodd_fermionic_parity(s1, n1) && _isodd_fermionic_parity(s2, n2)
      dval = -1.0
    end
    Tb = blockview(tensor(T), b)
    mat_dim = prod(dims(Tb)[1:2])
    Tbr = reshape(Tb, mat_dim, mat_dim)
    for i in diagind(Tbr)
      NDTensors.setdiagindex!(Tbr, dval, i)
    end
  end
  return T
end

# TODO: add a version that determines the sites
# from common site indices of ψ and A
"""
    setindex!(ψ::Union{MPS, MPO}, A::ITensor, r::UnitRange{Int};
              orthocenter::Int = last(r), perm = nothing, kwargs...)
    replacesites!([...])
    replacesites([...])

Replace the sites in the range `r` with tensors made
from decomposing `A` into an MPS or MPO.

The MPS or MPO must be orthogonalized such that
```
firstsite ≤ ITensors.orthocenter(ψ) ≤ lastsite
```

Choose the new orthogonality center with `orthocenter`, which
should be within `r`.

Optionally, permute the order of the sites with `perm`.
"""
function setindex!(
  ψ::MPST,
  A::ITensor,
  r::UnitRange{Int};
  orthocenter::Integer=last(r),
  perm=nothing,
  kwargs...,
) where {MPST<:AbstractMPS}
  # Replace the sites of ITensor ψ
  # with the tensor A, splitting up A
  # into MPS tensors
  firstsite = first(r)
  lastsite = last(r)
  @assert firstsite ≤ ITensors.orthocenter(ψ) ≤ lastsite
  @assert firstsite ≤ leftlim(ψ) + 1
  @assert rightlim(ψ) - 1 ≤ lastsite

  # TODO: allow orthocenter outside of this
  # range, and orthogonalize/truncate as needed
  @assert firstsite ≤ orthocenter ≤ lastsite

  # Check that A has the proper common
  # indices with ψ
  lind = linkind(ψ, firstsite - 1)
  rind = linkind(ψ, lastsite)

  sites = [siteinds(ψ, j) for j in firstsite:lastsite]

  #s = collect(Iterators.flatten(sites))
  indsA = filter(x -> !isnothing(x), [lind, Iterators.flatten(sites)..., rind])
  @assert hassameinds(A, indsA)

  # For MPO case, restrict to 0 prime level
  #sites = filter(hasplev(0), sites)

  if !isnothing(perm)
    sites0 = sites
    sites = sites0[[perm...]]
    # Check if the site indices
    # are fermionic
    if !using_auto_fermion() && any(anyfermionic, sites)
      if length(sites) == 2 && ψ isa MPS
        if all(allfermionic, sites)
          s0 = Index.(sites0)

          # TODO: the Fermionic swap is could be diagonal,
          # if we combine the site indices
          #C = combiner(s0[1], s0[2])
          #c = combinedind(C)
          #AC = A * C
          #AC = noprime(AC * _fermionic_swap(c))
          #A = AC * dag(C)

          FSWAP = _fermionic_swap(s0[1], s0[2])
          A = noprime(A * FSWAP)
        end
      elseif ψ isa MPO
        @warn "In setindex!(MPO, ::ITensor, ::UnitRange), " *
          "fermionic signs are only not handled properly for non-trivial " *
          "permutations of sites. Please inform the developers of ITensors " *
          "if you require this feature (otherwise, fermionic signs can be " *
          "put in manually with fermionic swap gates)."
      else
        @warn "In setindex!(::Union{MPS, MPO}, ::ITensor, ::UnitRange), " *
          "fermionic signs are only handled properly for permutations involving 2 sites. " *
          "The original sites are $sites0, with a permutation $perm. " *
          "To have the fermion sign handled correctly, we recommend performing your permutation " *
          "pairwise."
      end
    end
  end

  ψA = MPST(A, sites; leftinds=lind, orthocenter=orthocenter - first(r) + 1, kwargs...)
  #@assert prod(ψA) ≈ A

  ψ[firstsite:lastsite] = ψA

  return ψ
end

function setindex!(
  ψ::MPST, A::ITensor, r::UnitRange{Int}, args::Pair{Symbol}...; kwargs...
) where {MPST<:AbstractMPS}
  return setindex!(ψ, A, r; args..., kwargs...)
end

replacesites!(ψ::AbstractMPS, args...; kwargs...) = setindex!(ψ, args...; kwargs...)

replacesites(ψ::AbstractMPS, args...; kwargs...) = setindex!(copy(ψ), args...; kwargs...)

_number_inds(s::Index) = 1
_number_inds(s::IndexSet) = length(s)
_number_inds(sites) = sum(_number_inds(s) for s in sites)

"""
    MPS(A::ITensor, sites; kwargs...)
    MPO(A::ITensor, sites; kwargs...)

Construct an MPS/MPO from an ITensor `A` by decomposing it site
by site according to the site indices `sites`.

# Keywords

- `leftinds = nothing`: optional left dangling indices. Indices that are not in `sites` and `leftinds` will be dangling off of the right side of the MPS/MPO.
- `orthocenter::Integer = length(sites)`: the desired final orthogonality center of the output MPS/MPO.
- `cutoff`: the desired truncation error at each link.
- `maxdim`: the maximum link dimension.
"""
function (::Type{MPST})(
  A::ITensor, sites; leftinds=nothing, orthocenter::Integer=length(sites), kwargs...
) where {MPST<:AbstractMPS}
  N = length(sites)
  for s in sites
    @assert hasinds(A, s)
  end
  @assert isnothing(leftinds) || hasinds(A, leftinds)

  @assert 1 ≤ orthocenter ≤ N

  ψ = Vector{ITensor}(undef, N)
  Ã = A
  l = leftinds
  # TODO: To minimize work, loop from
  # 1:orthocenter and reverse(orthocenter:N)
  # so the orthogonality center is set correctly.
  for n in 1:(N - 1)
    Lis = IndexSet(sites[n])
    if !isnothing(l)
      Lis = unioninds(Lis, l)
    end
    L, R = factorize(Ã, Lis; kwargs..., tags="Link,n=$n", ortho="left")
    l = commonind(L, R)
    ψ[n] = L
    Ã = R
  end
  ψ[N] = Ã
  M = MPST(ψ)
  setleftlim!(M, N - 1)
  setrightlim!(M, N + 1)
  orthogonalize!(M, orthocenter)
  return M
end

function (::Type{MPST})(A::AbstractArray, sites; kwargs...) where {MPST<:AbstractMPS}
  return MPST(itensor(A, sites...), sites; kwargs...)
end

"""
    swapbondsites(ψ::Union{MPS, MPO}, b::Integer; kwargs...)

Swap the sites `b` and `b+1`.
"""
function swapbondsites(ψ::AbstractMPS, b::Integer; kwargs...)
  ortho = get(kwargs, :ortho, "right")
  ψ = copy(ψ)
  if ortho == "left"
    orthocenter = b + 1
  elseif ortho == "right"
    orthocenter = b
  end
  if leftlim(ψ) < b - 1
    orthogonalize!(ψ, b)
  elseif rightlim(ψ) > b + 2
    orthogonalize!(ψ, b + 1)
  end
  ψ[b:(b + 1), orthocenter=orthocenter, perm=[2, 1], kwargs...] = ψ[b] * ψ[b + 1]
  return ψ
end

"""
    movesite(::Union{MPS, MPO}, n1n2::Pair{Int, Int})

Create a new MPS/MPO where the site at `n1` is moved to `n2`,
for a pair `n1n2 = n1 => n2`.

This is done with a series a pairwise swaps, and can introduce
a lot of entanglement into your state, so use with caution.
"""
function movesite(
  ψ::AbstractMPS, n1n2::Pair{Int,Int}; orthocenter::Integer=last(n1n2), kwargs...
)
  n1, n2 = n1n2
  n1 == n2 && return copy(ψ)
  ψ = orthogonalize(ψ, n2)
  r = n1:(n2 - 1)
  ortho = "left"
  if n1 > n2
    r = reverse(n2:(n1 - 1))
    ortho = "right"
  end
  for n in r
    ψ = swapbondsites(ψ, n; ortho=ortho, kwargs...)
  end
  ψ = orthogonalize(ψ, orthocenter)
  return ψ
end

# Helper function for permuting a vector for the
# movesites function.
function _movesite(ns::Vector{Int}, n1n2::Pair{Int,Int})
  n1, n2 = n1n2
  n1 == n2 && return copy(ns)
  r = n1:(n2 - 1)
  if n1 > n2
    r = reverse(n2:(n1 - 1))
  end
  for n in r
    ns = replace(ns, n => n + 1, n + 1 => n)
  end
  return ns
end

function _movesites(ψ::AbstractMPS, ns::Vector{Int}, ns′::Vector{Int}; kwargs...)
  ψ = copy(ψ)
  N = length(ns)
  @assert N == length(ns′)
  for i in 1:N
    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
    ns = _movesite(ns, ns[i] => ns′[i])
  end
  return ψ, ns
end

# TODO: make a permutesites(::MPS/MPO, perm)
# function that takes a permutation of the sites
# p(1:N) for N sites
function movesites(ψ::AbstractMPS, nsns′::Vector{Pair{Int,Int}}; kwargs...)
  ns = first.(nsns′)
  ns′ = last.(nsns′)
  ψ = copy(ψ)
  N = length(ns)
  @assert N == length(ns′)
  p = sortperm(ns′)
  ns = ns[p]
  ns′ = ns′[p]
  ns = collect(ns)
  while ns ≠ ns′
    ψ, ns = _movesites(ψ, ns, ns′; kwargs...)
  end
  return ψ
end

# TODO: call the Vector{Pair{Int, Int}} version
function movesites(ψ::AbstractMPS, ns, ns′; kwargs...)
  ψ = copy(ψ)
  N = length(ns)
  @assert N == length(ns′)
  p = sortperm(ns′)
  ns = ns[p]
  ns′ = ns′[p]
  ns = collect(ns)
  for i in 1:N
    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
    ns = _movesite(ns, ns[i] => ns′[i])
  end
  return ψ
end

"""
    apply(o::ITensor, ψ::Union{MPS, MPO}, [ns::Vector{Int}]; kwargs...)
    product([...])

Get the product of the operator `o` with the MPS/MPO `ψ`,
where the operator is applied to the sites `ns`. If `ns`
are not specified, the sites are determined by the common indices
between `o` and the site indices of `ψ`.

If `ns` are non-contiguous, the sites of the MPS are
moved to be contiguous. By default, the sites are moved
back to their original locations. You can leave them where
they are by setting the keyword argument `move_sites_back`
to false.

# Keywords

- `cutoff::Real`: singular value truncation cutoff.
- `maxdim::Int`: maximum MPS/MPO dimension.
- `apply_dag::Bool = false`: apply the gate and the dagger of the gate (only relevant for MPO evolution).
- `move_sites_back::Bool = true`: after the ITensors are applied to the MPS or MPO, move the sites of the MPS or MPO back to their original locations.
"""
function product(
  o::ITensor,
  ψ::AbstractMPS,
  ns=findsites(ψ, o);
  move_sites_back::Bool=true,
  apply_dag::Bool=false,
  kwargs...,
)
  N = length(ns)
  ns = sort(ns)

  # TODO: make this smarter by minimizing
  # distance to orthogonalization.
  # For example, if ITensors.orthocenter(ψ) > ns[end],
  # set to ns[end].
  ψ = orthogonalize(ψ, ns[1])
  diff_ns = diff(ns)
  ns′ = ns
  if any(!=(1), diff_ns)
    ns′ = [ns[1] + n - 1 for n in 1:N]
    ψ = movesites(ψ, ns .=> ns′; kwargs...)
  end
  ϕ = ψ[ns′[1]]
  for n in 2:N
    ϕ *= ψ[ns′[n]]
  end
  ϕ = product(o, ϕ; apply_dag=apply_dag)
  ψ[ns′[1]:ns′[end], kwargs...] = ϕ
  if move_sites_back
    # Move the sites back to their original positions
    ψ = movesites(ψ, ns′ .=> ns; kwargs...)
  end
  return ψ
end

"""
    apply(As::Vector{<:ITensor}, M::Union{MPS, MPO}; kwargs...)
    product([...])

Apply the ITensors `As` to the MPS or MPO `M`, treating them as gates or matrices from pairs of prime or unprimed indices.

# Keywords

- `cutoff::Real`: singular value truncation cutoff.
- `maxdim::Int`: maximum MPS/MPO dimension.
- `apply_dag::Bool = false`: apply the gate and the dagger of the gate (only relevant for MPO evolution).
- `move_sites_back::Bool = true`: after the ITensor is applied to the MPS or MPO, move the sites of the MPS or MPO back to their original locations.

# Examples

Apply one-site gates to an MPS:
```julia
N = 3

ITensors.op(::OpName"σx", ::SiteType"S=1/2", s::Index) =
  2*op("Sx", s)

ITensors.op(::OpName"σz", ::SiteType"S=1/2", s::Index) =
  2*op("Sz", s)

# Make the operator list.
os = [("σx", n) for n in 1:N]
append!(os, [("σz", n) for n in 1:N])

@show os

s = siteinds("S=1/2", N)
gates = ops(os, s)

# Starting state |↑↑↑⟩
ψ0 = MPS(s, "↑")

# Apply the gates.
ψ = apply(gates, ψ0; cutoff = 1e-15)

# Test against exact (full) wavefunction
prodψ = apply(gates, prod(ψ0))
@show prod(ψ) ≈ prodψ

# The result is:
# σz₃ σz₂ σz₁ σx₃ σx₂ σx₁ |↑↑↑⟩ = -|↓↓↓⟩
@show inner(ψ, MPS(s, "↓")) == -1
```
Apply nonlocal two-site gates and one-site gates to an MPS:
```julia
# 2-site gate
function ITensors.op(::OpName"CX", ::SiteType"S=1/2", s1::Index, s2::Index)
  mat = [1 0 0 0
         0 1 0 0
         0 0 0 1
         0 0 1 0]
  return itensor(mat, s2', s1', s2, s1)
end

os = [("CX", 1, 3), ("σz", 3)]

@show os

# Start with the state |↓↑↑⟩
ψ0 = MPS(s, n -> n == 1 ? "↓" : "↑")

# The result is:
# σz₃ CX₁₃ |↓↑↑⟩ = -|↓↑↓⟩
ψ = apply(ops(os, s), ψ0; cutoff = 1e-15)
@show inner(ψ, MPS(s, n -> n == 1 || n == 3 ? "↓" : "↑")) == -1
```
Perform TEBD-like time evolution:
```julia
# Define the nearest neighbor term `S⋅S` for the Heisenberg model
function ITensors.op(::OpName"expS⋅S", ::SiteType"S=1/2",
                     s1::Index, s2::Index; τ::Number)
  O = 0.5 * op("S+", s1) * op("S-", s2) +
      0.5 * op("S-", s1) * op("S+", s2) +
            op("Sz", s1) * op("Sz", s2)
  return exp(τ * O)
end

τ = -0.1im
os = [("expS⋅S", (1, 2), (τ = τ,)),
      ("expS⋅S", (2, 3), (τ = τ,))]
ψ0 = MPS(s, n -> n == 1 ? "↓" : "↑")
expτH = ops(os, s)
ψτ = apply(expτH, ψ0)
```
"""
function product(
  As::Vector{<:ITensor}, ψ::AbstractMPS; move_sites_back::Bool=true, kwargs...
)
  Aψ = ψ
  for A in As
    Aψ = product(A, Aψ; move_sites_back=false, kwargs...)
  end
  if move_sites_back
    s = siteinds(Aψ)
    ns = 1:length(ψ)
    ñs = [findsite(ψ, i) for i in s]
    Aψ = movesites(Aψ, ns .=> ñs; kwargs...)
  end
  return Aψ
end

#
# QN functions
#

"""
    hasqns(M::MPS)

    hasqns(M::MPO)

Return true if the MPS or MPO has
tensors which carry quantum numbers.
"""
hasqns(M::AbstractMPS) = any(hasqns, data(M))

# Trait type version of hasqns
# Note this is not inferrable, so hasqns would be preferred
symmetrystyle(M::AbstractMPS) = symmetrystyle(data(M))

"""
    flux(M::MPS)

    flux(M::MPO)

    totalqn(M::MPS)

    totalqn(M::MPO)

For an MPS or MPO which conserves quantum
numbers, compute the total QN flux. For
a tensor network such as an MPS or MPO,
the flux is the sum of fluxes of each of
the tensors in the network. The name
`totalqn` is an alias for `flux`.
"""
function flux(M::AbstractMPS)
  hasqns(M) || return nothing
  q = QN()
  for j in (M.llim + 1):(M.rlim - 1)
    q += flux(M[j])
  end
  return q
end

totalqn(M::AbstractMPS) = flux(M)

function checkflux(M::AbstractMPS)
  for m in M
    checkflux(m)
  end
  return nothing
end

"""
    splitblocks[!](::typeof(linkinds), M::AbstractMPS; tol = 0)

Split the QN blocks of the links of the MPS or MPO into dimension 1 blocks. Then, only keep the blocks with `norm(b) > tol`.

This can make the ITensors of the MPS/MPO more sparse, and is particularly helpful as a preprocessing step on a local Hamiltonian MPO for DMRG.
"""
function splitblocks!(::typeof(linkinds), M::AbstractMPS; tol=0)
  for i in eachindex(M)[1:(end - 1)]
    l = linkind(M, i)
    if !isnothing(l)
      @preserve_ortho M begin
        M[i] = splitblocks(M[i], l)
        M[i + 1] = splitblocks(M[i + 1], l)
      end
    end
  end
  return M
end

function splitblocks(::typeof(linkinds), M::AbstractMPS; tol=0)
  return splitblocks!(linkinds, copy(M); tol=0)
end

#
# Broadcasting
#

BroadcastStyle(MPST::Type{<:AbstractMPS}) = Style{MPST}()
function BroadcastStyle(::Style{MPST}, ::DefaultArrayStyle{N}) where {N,MPST<:AbstractMPS}
  return Style{MPST}()
end

broadcastable(ψ::AbstractMPS) = ψ
function copyto!(ψ::AbstractMPS, b::Broadcasted)
  copyto!(data(ψ), b)
  # In general, we assume the broadcast operation
  # will mess up the orthogonality
  # TODO: special case for `prime`, `settags`, etc.
  reset_ortho_lims!(ψ)
  return ψ
end

function similar(bc::Broadcasted{Style{MPST}}, ElType::Type) where {MPST<:AbstractMPS}
  return similar(Array{ElType}, axes(bc))
end

function similar(bc::Broadcasted{Style{MPST}}, ::Type{ITensor}) where {MPST<:AbstractMPS}
  # In general, we assume the broadcast operation
  # will mess up the orthogonality so we use
  # a generic constructor where we don't specify
  # the orthogonality limits.
  return MPST(similar(Array{ITensor}, axes(bc)))
end

#
# Printing functions
#

function Base.show(io::IO, M::AbstractMPS)
  print(io, "$(typeof(M))")
  (length(M) > 0) && print(io, "\n")
  for i in eachindex(M)
    if !isassigned(M, i)
      println(io, "#undef")
    else
      A = M[i]
      if order(A) != 0
        println(io, "[$i] $(inds(A))")
      else
        println(io, "[$i] ITensor()")
      end
    end
  end
end

#
# Old code for adding MPS/MPO
#

#function plussers(::Type{T},
#                  left_ind::Index,
#                  right_ind::Index,
#                  sum_ind::Index) where {T<:Array}
#  total_dim    = dim(left_ind) + dim(right_ind)
#  total_dim    = max(total_dim, 1)
#  # TODO: I am not sure if we should be using delta
#  # tensors for this purpose? I think we should consider
#  # not allowing them to be made with different index sizes
#  #left_tensor  = δ(left_ind, sum_ind)
#  left_tensor  = diagITensor(1.0,left_ind, sum_ind)
#  right_tensor = ITensor(right_ind, sum_ind)
#  for i in 1:dim(right_ind)
#    right_tensor[right_ind(i), sum_ind(dim(left_ind) + i)] = 1
#  end
#  return left_tensor, right_tensor
#end
#
#function Base.:+(A::T, B::T; kwargs...) where {T <: AbstractMPS}
#  A = copy(A)
#  B = copy(B)
#  N = length(A)
#  length(B) != N && throw(DimensionMismatch("lengths of MPOs A ($N) and B ($(length(B))) do not match"))
#  orthogonalize!(A, 1; kwargs...)
#  orthogonalize!(B, 1; kwargs...)
#  C = similar(A)
#  rand_plev = 13124
#  lAs = [linkind(A, i) for i in 1:N-1]
#  prime!(A, rand_plev, "Link")
#
#  first  = Vector{ITensor}(undef,N-1)
#  second = Vector{ITensor}(undef,N-1)
#  for i in 1:N-1
#    lA = linkind(A, i)
#    lB = linkind(B, i)
#    r  = Index(dim(lA) + dim(lB), tags(lA))
#    f, s = plussers(typeof(data(A[1])), lA, lB, r)
#    first[i]  = f
#    second[i] = s
#  end
#  C[1] = A[1] * first[1] + B[1] * second[1]
#  for i in 2:N-1
#      C[i] = dag(first[i-1]) * A[i] * first[i] + dag(second[i-1]) * B[i] * second[i]
#  end
#  C[N] = dag(first[N-1]) * A[N] + dag(second[N-1]) * B[N]
#  prime!(C, -rand_plev, "Link")
#  truncate!(C; kwargs...)
#  return C
#end


text:
abstract type Forecast <: TimeSeriesData end

# Subtypes of Forecast must implement
# - get_count
# - get_horizon
# - get_initial_times
# - get_initial_timestamp
# - get_name
# - get_scaling_factor_multiplier
# - get_window
# - iterate_windows

Base.length(ts::Forecast) = get_count(ts)

abstract type AbstractDeterministic <: Forecast end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic.
function eltype_data_common(forecast::Forecast)
    return eltype(first(values(get_data(forecast))))
end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic.
function get_count_common(forecast)
    return length(get_data(forecast))
end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic. Allows for optimized execution.
function get_horizon_common(forecast)
    return length(first(values(get_data(forecast))))
end

"""
Return the initial times in the forecast.
"""
function get_initial_times(f::Forecast)
    return get_initial_times(get_initial_timestamp(f), get_count(f), get_interval(f))
end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic. Allows for optimized execution.
function get_initial_times_common(forecast::Forecast)
    return keys(get_data(forecast))
end

"""
Return the total period covered by the forecast.
"""
function get_total_period(f::Forecast)
    return get_total_period(
        get_initial_timestamp(f),
        get_count(f),
        get_interval(f),
        get_horizon(f),
        get_resolution(f),
    )
end

"""
Return the forecast window corresponsing to interval index.
"""
function get_window(forecast::Forecast, index::Int; len = nothing)
    return get_window(forecast, index_to_initial_time(forecast, index); len = len)
end

function iterate_windows_common(forecast)
    return (get_window(forecast, it) for it in keys(get_data(forecast)))
end

"""
Return the Dates.DateTime corresponding to an interval index.
"""
function index_to_initial_time(forecast::Forecast, index::Int)
    return get_initial_timestamp(forecast) + get_interval(forecast) * index
end

"""
Return a TimeSeries.TimeArray for one forecast window.
"""
function make_time_array(
    forecast::Forecast,
    start_time::Dates.DateTime;
    len::Union{Nothing, Int} = nothing,
)
    return get_window(forecast, start_time; len = len)
end

function make_timestamps(forecast::Forecast, initial_time::Dates.DateTime, len = nothing)
    if len === nothing
        len = get_horizon(forecast)
    end

    return range(initial_time; length = len, step = get_resolution(forecast))
end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic. Allows for optimized execution.
function get_initial_timestamp_common(forecast)
    return first(keys(get_data(forecast)))
end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic. Allows for optimized execution.
function get_interval_common(forecast)
    its = get_initial_times(forecast)
    if length(its) == 1
        return Dates.Second(0)
    end
    first_it, state = iterate(its)
    second_it, state = iterate(its, state)
    return second_it - first_it
end

# This method requires that the forecast type implement a `get_data` method like
# Deterministic.
function get_window_common(
    forecast,
    initial_time::Dates.DateTime;
    len::Union{Nothing, Int} = nothing,
)
    horizon = get_horizon(forecast)
    if len === nothing
        len = horizon
    end

    data = get_data(forecast)[initial_time]
    if ndims(data) == 2
        # This is necessary because the Deterministic and Probabilistic are 3D Arrays
        # We need to do this to make the data a 2D TimeArray. In a get_window the data is always count = 1
        @assert_op size(data)[1] <= len
        data = @view data[1:len, :]
    else
        data = @view data[1:len]
    end

    return TimeSeries.TimeArray(make_timestamps(forecast, initial_time, len), data)
end

ext: jl


text:
"""
    process_bouts(df::AbstractDataFrame)

Take as input a processed pokes datafame `df` of 1 single session and returns a dataframe were each row are information about a poking bout.
A poking bout ends either when a reward occurs or the animal changes side
"""
function process_bouts(df::AbstractDataFrame)
    vals_per_bout = [:Stim, :StimFreq, :Wall, :Protocol, :Block, :StreakInBlock, :Side, :ReverseStreak]
    vals_per_day = [:Box,:StimDay];
    bout_table = combine(groupby(df, :Bout)) do dd
        dt = DataFrame(
        Omissions_plus_one = size(dd,1),
        Leave = dd[end,:Leave],
        Start_Reward = dd[1,:Reward],
        Side = dd[1,:Side],
        Stim = dd[1,:Stim],
        Trial_duration = (dd[end,:PokeOut]-dd[1,:PokeIn]),
        PreInterpoke = size(dd,1) > 1 ? maximum(skipmissing(dd[!,:PreInterpoke])) : missing,
        PostInterpoke = size(dd,1) > 1 ? maximum(skipmissing(dd[!,:PostInterpoke])) : missing,
        Wall = dd[1,:Wall],
        Correct_start = dd[1,:Correct],
        Block = dd[1,:Block],
        StreakInBlock = dd[1,:StreakInBlock],
        Streak = dd[1,:Streak],
        ReverseStreak = dd[1,:ReverseStreak],
        CorrectStart = dd[1,:Correct],
        CorrectLeave = !dd[end,:Correct]
        )

        for v in vals_per_bout
            if v in propertynames(dd)
                dt[!,v] .= dd[1, v]
            end
        end

        return dt
    end

    for v in vals_per_day
        if v in propertynames(df)
            bout_table[!,v] .= df[1, v]
        end
    end

    return bout_table
end


text:
using PyCall

const kwant = PyNULL()
const plotter = PyNULL()
const operator = PyNULL()
const physics = PyNULL()

function __init__()
    copy!(kwant, pyimport("kwant"))
    copy!(plotter, pyimport("kwant.plotter"))
    copy!(operator, pyimport("kwant.operator"))
    copy!(physics, pyimport("kwant.physics"))
end


text:
const Z3Expr = Z3.Expr

Tactic(s::String) = Z3.Tactic(main_ctx(), s)

function mk_solver()
    t = Tactic("simplify") & Tactic("propagate-values") & Tactic("solve-eqs") &
        Tactic("purify-arith") & Tactic("elim-term-ite") & Tactic("simplify") &
        Tactic("tseitin-cnf") & Tactic("nlsat")
    Z3.mk_solver(t)
end

function Int(x::Z3.Expr)
    @assert is_int(x)
    get_numeral_int(x)
end
Rational{Int}(x::Z3.Expr) = Int(numerator(x)) // Int(denominator(x))

function parse_model(m::Model)
    nlmodel = Dict{Symbol, Rational{Int}}()
    for (k,v) in consts(m)
        sym = Symbol(string(k))
        if is_int(v)
            push!(nlmodel, sym=>Int(v))
        elseif is_algebraic(v)
            approx = String(get_decimal_string(v, 10))[1:end-1]
            @warn "Algebraic numbers not yet supported, got $(v), returning approximation $(approx)"
            push!(nlmodel, sym=>parse(BigFloat, approx))
        elseif is_real(v)
            push!(nlmodel, sym=>try_int(Rational{Int}(v)))
        end
    end
    @debug "" nlmodel
    nlmodel
end

# ------------------------------------------------------------------------------

Z3.CxxWrap.@cxxdereference function Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, ex::XExpr)
    t = postwalk(ex) do x
        if x isa Float64
            # if isinteger(x)
            #     convert(Int, x)
            # else
                y = rationalize(x)
                real_val(ctx, numerator(y), denominator(y))
            # end
        elseif x isa Rational
            real_val(ctx, numerator(y), denominator(y))
        elseif x isa Int
            real_val(ctx, x)
        elseif issymbol(x)
            get(vs, x, x)
        else
            x
        end
    end
    eval(t)
end

Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Constraint{EQ})  = Z3Expr(ctx, vs, c.poly) == 0
Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Constraint{NEQ}) = Z3Expr(ctx, vs, c.poly) != 0
Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Constraint{LT})  = Z3Expr(ctx, vs, c.poly) <  0
Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Constraint{LEQ}) = Z3Expr(ctx, vs, c.poly) <= 0
Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Constraint{GT})  = Z3Expr(ctx, vs, c.poly) >  0
Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Constraint{GEQ}) = Z3Expr(ctx, vs, c.poly) >= 0
Z3Expr(ctx::Context, vs::Dict{Symbol,Z3Expr}, c::Clause) =
    length(c) > 1 ? or([Z3Expr(ctx, vs, x) for x in c]) : Z3Expr(ctx, vs, first(c))


text:
using StrFs, BenchmarkTools

primitive type String6 <: AbstractString 48 end
primitive type Bits48 48 end
primitive type Bits64 64 end
primitive type Bits32 32 end
primitive type Bits96 96 end

String(s::String6) = begin
    String(reinterpret(UInt8, [s]))
end

import Base:iterate
iterate(s::String6) = iterate(String(s))
iterate(s::String6, i) = iterate(String(s), i)

gender = [strf"abc", strf"female", strf"abcdefgh"]

#@time ass = ShortString7.(a);

#iterate(s::StrF, i::Integer) = iterate(String(s), i)
#iterate(s::StrF) = iterate(String(s))

using CSV
df = DataFrame(
    a=a,
    as = as
    );

@time CSV.write("c:/data/a.csv", df);

# Base.summarysize(a)
# Base.summarysize(as)

using BenchmarkTools, SortingLab, SortingAlgorithms

@time sort!(as, alg = RadixSort)

@elapsed a_sorted = SortingLab.fsort(a)
@elapsed as_sorted = sort(as, alg=RadixSort)
@elapsed ass_sorted = ShortStrings.fsort(ass)

using Blosc

@time ass_compressed = Blosc.compress(ass)

@time Base.summarysize(ass_compressed)

Ptr{UInt}(Ref(gender[1]))

unsafe_pointer_to_objref(gender[1])

Ptr{StrF{8}}

ss = gender[1]

reinterpret(UInt, UInt8.(ss.bytes))

using Blosc

compress(gender)

uint_mapping(_, s::StrF{S}) where S = begin
    Int(ceil(S/8))*8
end

uint_mapping(1, gender[1])

using SortingAlgorithms
sort(gender, alg = RadixSort)

reinterpret(String6, gender)

ptrs = Vector{UInt8}(undef, 8) |> pointer
x = unsafe_load(Ptr{Bits96}(ptrs))

reinterpret(UInt8, [x])

x = "abc"
using BenchmarkTools
@btime Base.zext_int(UInt64, unsafe_load(Ptr{Bits48}(pointer(x))))

UInt8.(Vector{Char}(x)

reinterpret(String6, UInt8[2, 3, 4, 6, 8, 9])

iterate(s::String6) = s |

Vector{String6}(undef, 222)

unsafe_load(Ptr{String6})


text:
abstract type Policy end

function policy(p::Policy, traj::ContactTraj, t::T) where T
    @warn "policy not defined"
    return nothing
end

"""
    no policy
"""
struct NoPolicy{T} <: Policy
    u::Vector{T}
end

function no_policy(model::ContactModel)
    NoPolicy(zeros(model.dim.u))
end

function policy(p::NoPolicy, x, traj::ContactTraj, t)
    return p.u
end

"""
    open-loop policy
"""
mutable struct OpenLoop{U} <: Policy
    u::Vector{U} # nominal controls
    idx::Int
    cnt::Int
    N_sample::Int
end

function open_loop_policy(u; N_sample = 1)
    OpenLoop(u, 0, N_sample, N_sample)
end

function policy(p::OpenLoop, x, traj, t)
    # reset
    if t == 1
        p.idx = 0
        p.cnt = p.N_sample
    end

    if p.cnt == p.N_sample
        p.idx += 1
        p.cnt = 0
    end

    p.cnt += 1

    return p.u[p.idx] ./ p.N_sample
end

"""
    control saturation
"""

control_saturation(u, uL, uU) = min.(max.(uL, u), uU)


text:
using Plots
pyplot()
using SolverBenchmark
using DataFrames

# We consider as solved the really solved (first_order) but also
# unbounded parameters (infeasible) or function values (unbounded)
# as well as exceptions (NaN in the evaluation of the model).

solved(df) = (df.status .== :first_order) .| (df.status .== :infeasible) .| (df.status .== :unbounded).| (df.status .== :exception)

costnames = ["time",
             "objective evals",
             "derivative evals",
             "second derivative evals",
             "obj + der + 2nd der"]

costs = [df -> .!solved(df) .* Inf .+ df.elapsed_time,
         df -> .!solved(df) .* Inf .+ df.neval_obj,
         df -> .!solved(df) .* Inf .+ df.neval_grad,
         df -> .!solved(df) .* Inf .+ df.neval_hess,
         df -> .!solved(df) .* Inf .+ df.neval_obj .+ df.neval_grad .+ df.neval_hprod]

#include("my_profile_tools.jl")

# profile on the sum of evaluations obj, grad and hess
#p = my_performance_profile(stats,costs[1], title = costnames[1])

# array of profiles
#p = my_profile_solvers(pruned_stats, costs, costnames)

#Plots.pdf(p, "test")


text:
for N in [Float64, Rational{Int}]
    # ====================================
    # Concrete convex hull for point sets
    # ====================================

    # corner cases in dimension 1
    @test convex_hull([Vector{N}(undef, 0)]) == [Vector{N}(undef, 0)]
    @test convex_hull([[N(0)]]) == [[N(0)]]
    @test ispermutation(convex_hull([N[2], N[1]]), [N[2], N[1]])
    @test convex_hull([[N(2)], [N(2)]]) == [[N(2)]]

    # corner cases in dimension 2
    @test convex_hull([[N(0), N(0)]]) == [[N(0), N(0)]]
    @test convex_hull([[N(1), N(0)], [N(0), N(1)]]) == [[N(1), N(0)], [N(0), N(1)]]
    @test convex_hull([[N(1), N(0)], [N(1), N(0)]]) == [[N(1), N(0)]]

    # test corner cases with one and two vectors (see #876)
    p1 = [1., 2.]
    p2 = [1., 3.]
    @test convex_hull([p1]) == [p1]
    @test ispermutation(convex_hull([p1, p2]), [p1, p2])

    # corner cases in higher dimension
    @test convex_hull([[N(0), N(0), N(0)]]) == [[N(0), N(0), N(0)]]
    # there is no sorting in higher dim
    @test convex_hull([[N(1), N(0), N(0)], [N(0), N(1), N(0)]]) == [[N(1), N(0), N(0)], [N(0), N(1), N(0)]]

    # dimension 1
    points_1D = [[N(2)], [N(1)], [N(1/2)], [N(-5)]]
    @test convex_hull(points_1D) == [[N(-5)], [N(2)]]
    convex_hull!(points_1D) # check in-place version
    @test points_1D == [[N(-5)], [N(2)]]

    # dimension 2
    points_2D = [[N(1), N(0)], [N(1), N(1)], [N(0), N(1)], [N(-1), N(-1)], [N(1/2), N(1/2)]]
    @test convex_hull(points_2D) == [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]
    convex_hull!(points_2D) # check in-place version
    @test points_2D == [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]

    # three vertex case in 2 dimensions
    ccw_points = [N[1, 1], N[-1, 1], N[-1, 0]]
    ccw_p = convex_hull!(ccw_points)
    ccw_expr = [N[1, 1], N[-1, 1], N[-1, 0]]
    @test is_cyclic_permutation(ccw_p, ccw_expr) # points sorted in a counter-clockwise fashion
    cw_points = [N[-1, 1], N[1, 1], N[-1, 0]]
    cw_p = convex_hull!(cw_points)
    cw_expr = [N[1, 1], N[-1, 1], N[-1, 0]]
    @test is_cyclic_permutation(cw_p, cw_expr) # points sorted in clockwise fashion
    @test ispermutation(convex_hull!([N[1, 1], N[2, 2], N[3, 3]]), [N[1, 1], N[3, 3]]) # points aligned
    @test ispermutation(convex_hull!([N[0, 1], N[0, 2], N[0, 3]]), [N[0, 1], N[0, 3]]) # three points on a vertical line
    @test convex_hull!([N[0, 1], N[0, 1], N[0, 1]]) == [N[0, 1]] # three equal points

    # four vertex case in 2 dimentions
    A = N[1, 0]
    B = N[1, 1]
    C = N[-1, 1]
    D = N[-1, 0]
    expr = [A, B, C, D]
    points = [A, B, C, D]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ABCD
    points = [A, D, C, B]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ADCB
    points = [A, B, D, C]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ABDC
    points = [A, D, B, C]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ADBC
    points = [D, A, C, B]
    @test is_cyclic_permutation(convex_hull!(points), expr) # DACB
    points = [A, C, D, B]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ACDB
    A = N[0, 1]
    B = N[-1, -1]
    C = N[1, -1]
    D = N[0, 0]
    expr = [A, B, C]
    points = [A, B, C, D]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ABC
    points = [A, C, B, D]
    @test is_cyclic_permutation(convex_hull!(points), expr) # CBA
    points = [D, B, C, A]
    @test is_cyclic_permutation(convex_hull!(points), expr) # BCD
    points = [D, C, B, A]
    @test is_cyclic_permutation(convex_hull!(points), expr) # DCB
    points = [B, D, C, A]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ACD
    points = [C, D, B, A]
    @test is_cyclic_permutation(convex_hull!(points), expr) # CAD
    points = [B, C, D, A]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ABD
    points = [C, B, D, A]
    @test is_cyclic_permutation(convex_hull!(points), expr) # ADB
    @test ispermutation(convex_hull!([N[1, 1], N[2, 2], N[3, 3], N[4, 4]]), [N[1, 1], N[4, 4]]) # points aligned
    @test  is_cyclic_permutation(convex_hull([[1., 1], [0., 0], [0., 1], [0., 2]]), [[1, 1.], [0., 2], [0., 0]]) # test that issue #1661 is fixed

    # five-vertices case in 2D
    points = [N[1//10, 3//10], N[1//5, 1//10], N[3//10, 7//25], N[2//5, 3//5],
              N[9//10, 1//5]]
    sorted = [N[1//10, 3//10], N[1//5, 1//10], N[9//10, 1//5], N[2//5, 3//5]]
    for algorithm in ["monotone_chain", "monotone_chain_sorted"]
        points_copy = copy(points)
        @test is_cyclic_permutation(
            convex_hull!(points_copy, algorithm=algorithm), sorted)
    end
    @test_throws ErrorException convex_hull!(points, algorithm="")

    # single-vertex case in 2D (five vertices to test the general algorithm)
    p = N[0, 0]
    @test convex_hull([p, p, p, p, p]) == [p]

    # higher dimension
    if test_suite_polyhedra && N != Float32 # no backend supporting Float32
        points_3D = [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)],
                     [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)],
                     [N(0), N(0), N(0)], [N(1), N(2), N(3)]]
        sorted = [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)],
            [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]]
        @test ispermutation(convex_hull(points_3D), sorted)
        convex_hull!(points_3D) # check in-place version
        @test ispermutation(points_3D, sorted)
    end

    # ============================
    # Binary concrete convex hull
    # ============================

    V1 = VPolygon([[N(1), N(0)], [N(1), N(1)], [N(0), N(1)]])
    V2 = VPolygon([[N(-1), N(-1)], [N(1/2), N(1/2)]])
    ch = convex_hull(V1, V2)
    @test ispermutation(vertices_list(ch), [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]])

    if test_suite_polyhedra && N != Float32 # no backend supporting Float32
        V1 = VPolytope([[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)]])
        V2 = VPolytope([[N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(0), N(0), N(0)],
                       [N(1), N(2), N(3)]])
        ch = convex_hull(V1, V2)
        @test ch isa VPolytope
        @test ispermutation(vertices_list(ch), [[N(1), N(0), N(4)], [N(1), N(1), N(5)],
                                    [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)],
                                    [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]])
    end

    # general LazySets
    for (n, T) in [(1, Interval), (2, VPolygon), (3, VPolytope)]
        if n > 2 && !test_suite_polyhedra  # dimensions > 2 require Polyhedra
            continue
        end
        X = BallInf(zeros(N, n), N(1))
        Y = Ball1(ones(N, n), N(3//2))
        ch = convex_hull(X, Y)
        @test ch isa T{N}
    end
    S = Singleton(N[1])
    ch = convex_hull(S, S)
    @test ch isa Singleton{N} && element(ch) == element(S)
    U = UnionSet(Interval(N(1), N(2)), Interval(N(3), N(4)))
    @test convex_hull(U) == Interval(N(1), N(4))

    # UnionSetArray
    V1 = VPolytope([N[0, 0], N[1, 0], N[0, 1]])
    V2 = VPolytope([N[1, 1], N[1, 0], N[0, 1]])
    U = UnionSetArray([V1, V2]);
    @test ispermutation(convex_hull(U), [N[0, 0], N[1, 0], N[0, 1], N[1, 1]])
end


text:
module LibCloud

include("common.jl")
include("lazyhelp.jl")
include("storage.jl")
include("dns.jl")
include("loadbalancer.jl")
include("compute.jl")
include("containers.jl")
include("backup.jl")

end # module


text:
# L-infinity norm

"""
  NormLinf(λ::Real=1.0)

Returns the function `g(x) = λ⋅max(abs(x))`, for a nonnegative parameter `λ ⩾ 0`.
"""

NormLinf{R <: Real}(lambda::R=1.0) = Postcompose(Conjugate(IndBallL1(one(R))), lambda)

function (f::Conjugate{IndBallL1{R}}){R <: Real, S <: RealOrComplex}(x::AbstractArray{S})
  # follows from the definition of conjugate function and properties of norms
  # although it is not really needed since with the above constructor one makes f.f.r = 1.0
  # but just in case one constructs the conjugate of a different L1 ball...
  return vecnorm(x, Inf)*(f.f.r)
end

fun_name{R <: Real}(f::Postcompose{Conjugate{IndBallL1{R}}, R}) = "weighted L-infinity norm"
fun_expr{R <: Real}(f::Postcompose{Conjugate{IndBallL1{R}}, R}) = "x ↦ λ||x||_∞ = λ⋅max(abs(x))"
fun_params{R <: Real}(f::Postcompose{Conjugate{IndBallL1{R}}, R}) = "λ = $(f.a*(f.f.f.r))"


text:
## the learning algorithm: PLS2 - multiple targets
function trainer(model::PLS2Model{T},
                X::AbstractArray{T}, Y::Matrix{T}) where T<:AbstractFloat
    W,b,P,Q  = model.W,model.b,model.P,model.Q
    model.ntargetcols = size(Y,2)
    nfactors = model.nfactors

    for i = 1:nfactors
        b[:,i]    = Y[:,1] #u: arbitrary col. Thus, I set to the first.
        Rold = b[:,i]
        local R::Vector{T}
        while true

            W[:,i]  = X'b[:,i]
            W[:,i] /= norm(W[:,i])#sqrt.(W[:,i]'*W[:,i])
            R       = X*W[:,i]

            Rold      = R
            Q[:,i]    = Y'R
            Q[:,i]   /= norm(Q[:,i])
            b[:,i]    = Y*Q[:,i]
            if  all(abs.(R - Rold) .<= 1.0e-3)
                break
            end

        end
        Rn     = R'/(R'R) # change to use function...
        P[:,i] = Rn*X

        X      = X - R * P[:,i]'
        c      = Rn'b[:,i:i]'./(R'R)
        Yp     = c*R*Q[:,i]'
        Y      = Y - Yp
        if  all(abs.(Y - Yp) .<= 1.0e-3)
            print("PLS2 converged. No need learning with more than $(i) factors")
            model.nfactors = i
            return model
        end

    end

    return model

end

function predictor(model::PLS2Model{T},
                   X::AbstractArray{T}) where T<:AbstractFloat

    W,Q,P    = model.W,model.Q,model.P
    nfactors   = model.nfactors
    Y          = zeros(T,size(X,1),model.ntargetcols)
    #println("nfactors: ",nfactors)
    for i = 1:nfactors
        R      = X*W[:,i]
        X      = X - R * P[:,i]'
        Y      = Y + R * Q[:,i]'
    end

    return Y
end


text:
export compute_rho,  compute_e
export compute_rho!, compute_e!

import Statistics:mean
import FFTW:fft, ifft, fft!, ifft!

"""

    compute_rho( mesh, f)

    Compute charge density

    ρ(x,t) = ∫ f(x,v,t) delta2

    return ρ - ρ̄ if neutralized=true

"""
function compute_rho(mesh2::UniformMesh, f, neutralized=true)

   local delta2  = mesh2.step
   ρ = delta2 * sum(f, dims=2)
   if (neutralized)
       ρ .- mean(ρ)
   else
       ρ
   end

end

"""

    compute_rho!(rho, mesh, f)

    Inplace computation of charge density

"""
function compute_rho!(rho::Vector{Complex{Float64}},
		      mesh2::UniformMesh,
		      f::Array{Complex{Float64},2})

   local delta2 = mesh2.step
   rho .= delta2 * vec(sum(f, dims=2))
   rho .= rho .- mean(rho)

end

"""

    compute_e( mesh, ρ)

    Compute 1d electric field using that -ik * e = ρ

"""
function compute_e(mesh1::UniformMesh, ρ)

   n1 = mesh1.length
   k =  2π / (mesh1.stop - mesh1.start)
   modes  = zeros(Float64, n1)
   modes .= k * vcat(0:n1÷2-1,-n1÷2:-1)
   modes[1] = 1.0
   ρ̂ = fft(ρ)./modes
   vec(real(ifft(-1im*ρ̂)))

end

"""

    compute_e!( e, mesh, ρ)

    Inplace computation of electric field

"""
function compute_e!(e::Vector{Complex{Float64}},
		    mesh1::UniformMesh,
		    ρ::Vector{Complex{Float64}})

   n1 = mesh1.length
   k =  2π / (mesh1.stop - mesh1.start)
   modes  = zeros(Float64, n1)
   modes .= k * vcat(0:n1÷2-1,-n1÷2:-1)
   modes[1] = 1.0
   fft!(ρ)
   e .= -1im * ρ ./ modes
   ifft!(e)

end

export RectMesh2D

"""

    mesh = RectMesh2D( mesh1, mesh2)

    Data type that represents a 2D Rectangular Mesh

"""
struct RectMesh2D

    x1min :: Float64
    x1max :: Float64
    n1   :: Int64
    delta1   :: Float64
    ymin :: Float64
    ymax :: Float64
    n2   :: Int64
    delta2   :: Float64

    function RectMesh2D( mesh1::UniformMesh, mesh2::UniformMesh )

        x1min = mesh1.start
        x1max = mesh1.stop
        n1   = mesh1.length
        delta1   = mesh1.step

        ymin = mesh2.start
        ymax = mesh2.stop
        n2   = mesh2.length
        delta2   = mesh2.step

        new( x1min, x1max, n1, delta1, ymin, ymax, n2, delta2)

    end

end

export meshgrid
"""

    X, Y = meshgrid( x, y )

    Utility function to mimic the numpy function

"""
function meshgrid( x, y )

   n1, n2 = size(x)[1], size(y)[1]
   repeat(x,1,n2), repeat(y',n1,1)

end

export poisson!
"""

   poisson!(ρ, mesh, ex, ey)

   Solve the equation Δ Φ = - ρ

   ex = ∂ Φ / ∂ x
   ey = ∂ Φ / ∂ y

   Be careful the ρ array is destroyed

"""
function poisson!( ρ::Array{Complex{Float64},2},
		   mesh::RectMesh2D,
		   ex::Array{Complex{Float64},2},
		   ey::Array{Complex{Float64},2} )

    kx0 =  2π / (mesh.x1max - mesh.x1min)
    ky0 =  2π / (mesh.ymax - mesh.ymin)

    fft!(ρ,[1,2])

    kx = kx0 * vcat(0:mesh.n1÷2-1,-mesh.n1÷2:-1)
    kx[1] = 1.0
    ky = ky0 * vcat(0:mesh.n2÷2-1,-mesh.n2÷2:-1)
    kx[1] = 1.0

    for i = 1:mesh.n1
       kx2 = kx[i]*kx[i]
       for j =  1:mesh.n2÷2+1
          k2 = kx2 +ky[j]*ky[j]
          ex[i,j] = -1im * kx[i]/k2 * ρ[i,j]
	  ey[i,j] = -1im * ky[j]/k2 * ρ[i,j]
       end
       for j = mesh.n2÷2+2:mesh.n2
          k2 = kx2 +ky[j]*ky[j]
          ex[i,j] = -1im * kx[i]/k2 * ρ[i,j]
          ey[i,j] = -1im * ky[j]/k2 * ρ[i,j]
       end
    end

    ifft!(ex,[1,2])
    ifft!(ey,[1,2])

end


text:
function BondNSxPs(m,b0,b1,b2,tau,b3=0.0,tau2=1.0)
#BondNSxPs    Extended Nelson and Siegel (1987) spot rate, forward rate, and discount function
#
#
#
#
#  Usage:    (s,f,d) = BondNSxPs(m,b0,b1,b2,tau[,b3[,tau2]])
#
#  Input:    m         NxK matrix  times to maturity in years (don't use less than 1e-8 or so)
#            b0        scalar, beta0 parameter in NS
#            b1        scalar, beta1 parameter in NS
#            b2        scalar, beta2 parameter in NS
#            tau       scalar, tau   parameter in NS
#            b3        (optional) scalar, beta2 parameter in extended NS [0]
#            tau2      (optional) scalar, tau2  parameter in extended NS [1]
#
#  Output:   s         NxK matrix,  spot rate (continously compounded)
#            f         NxK matrix,  forward rate (continously compounded)
#            d         NxK matrix,  discount function
#
#
#
#  Note:     For more details, see
#            (a) Svensson (1995), "Estimating Forward Interest Rates with
#                the Extended Nelson & Siegel Method," Quarterly Review,
#                Sveriges Riksbank, 1995:3, 13-26.
#            (b) Soderlind and Svensson (1997), "New Techniques to Extract
#                Market Expectations from Financial Instruments,"
#                Journal of Monetary Economics 40, 383-429.
#
#
#
#  Paul.Soderlind@unisg.ch, April 2002, to Julia 2015
#----------------------------------------------------------------------------*/

                                                 #forward rate
  f =  b0 .+ b1*exp.(-m/tau) + b2*(m/tau).*exp.(-m/tau) + b3*(m/tau2).*exp.(-m/tau2)

  s =  b0 .+ b1* (1 .- exp.(-m/tau)) ./(m/tau) +
             b2*((1 .- exp.(-m/tau)) ./(m/tau)  - exp.(-m/tau)) +
             b3*((1 .- exp.(-m/tau2))./(m/tau2) - exp.(-m/tau2))        #spot rate

  d = exp.(-s.*m)             #discount function

  return s,f,d

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
function BondNSxEstPs(par0,Q,tm,c,s0,ytmLoss=0,weight=1.0)
#BondNSxEstPs  Estimates parameters in (extended) Nelson-Siegel yield curve model by
#              non-linear least squares (minimizing squared differences between
#              actual and fitted bond prices or yield to maturities).
#
#
#
#  Usage:    NSb = BondNSxEstPs(par0,Q,tm,c,s0[,ytmLoss[,weight]])
#
#
#  Input:    par0       4x1 or 6x1 vector, initial guess of paremeters
#                         if 4x1: standard NS with par0 = [b0,b1,b2,tau]
#                         if 6x1: extended NS with par0 = [b0,b1,b2,tau,b3,tau2]
#            Q          nx1 vector, data on bond prices, eg. 1.01
#            tm         nx1 vector, data on time to maturity (in years), eg. 2.54
#            c          nx1 vector, data on bond coupons, eg. 0.06
#            s0         scalar, restricted value of short rate:
#                       if [] no restriction, else b0 + b1 = s0 is imposed
#            ytmLoss    scalar, scalar, 0: mimimize squared price errors;
#                                       1: minimize squared ytm errors (default 0)
#            weight     scalar or nx1 vector, weights in the loss function
#
#  Output:   NSb        4x1 or 6x1 vector, estimated parameters
#
#
#
#  Note:     For more details, see
#            (a) Svensson (1995), "Estimating Forward Interest Rates with
#                the Extended Neolson & Siegel Method," Quarterly Review,
#                Sveriges Riksbank, 1995:3, 13-26.
#            (b) Soderlind and Svensson (1997), "New Techniques to Extract
#                Market Expectations from Financial Instruments,"
#                Journal of Monetary Economics 40, 383-429.
#
#
#  Uses:  Optim
#
#
#
#  Paul.Soderlind@unisg.ch, April 2002, to Julia Nov 2015
#------------------------------------------------------------------------------

  Qtc = [Q tm c]

  if in(length(par0),[4,6]) && !isempty(s0)
      par0 = deleteat!(par0,2)   #b1 = s0 - b0 is imposed by BondNSxLossPs
  end

  Sol = optimize(b->BondNSxLossPs(b,Qtc,s0,ytmLoss,weight),par0)
  NSb = Optim.minimizer(Sol)
  if !Optim.converged(Sol)
    println("no convergence")
    return
  end

  if length(NSb) == 3
    NSb[[1,3]] = abs.(NSb[[1,3]])
    NSb        = [NSb[1]; s0-NSb[1]; NSb[2:3]]
  elseif length(NSb) == 4
    NSb[[1,4]] = abs.(NSb[[1,4]])
  elseif length(NSb) == 5
    NSb[[1,3,5]] = abs.(NSb[[1,3,5]])
    NSb = [NSb[1]; s0-NSb[1]; NSb[2:5]]
  elseif length(NSb) == 6
    NSb[[1,4,6]] = abs.(NSb[[1,4,6]])
  end

  return NSb

end
#----------------------------------------------------------------------------

#------------------------------------------------------------------------------
function BondNSxLossPs(b,Qtc,s0,ytmLoss=0,weight=1.0)
#BondNSxLossPs    Defines loss function for bond prices in extended Nelson-Siegel model.
#                 Used for estimation of the parameters in the model by minimizing
#                 Loss (squared price deviations, possibly with weights).
#
#
#  Usage:    Loss = BondNSxLossPs(b,Qtc,s0,ytmLoss,weight)    or
#                 = BondNSxLossPs(b,Qtc,s0,ytmLoss)           or
#                 = BondNSxLossPs(b,Qtc,s0)
#
#  Input:    b         3x1, 4x1, 5x1, or 6x1 vector with parameters in Nelson-Siegel model
#                        if 3x1: [b0,b2,tau], NS with restriction that b1 = g_s0 - b0
#                        if 4x1: [b0,b1,b2,tau], NS without restrictions
#                        if 5x1: [b0,b2,tau,b3,tau2], extended NS with restriction that b1 = g_s0 - b0
#                        if 6x1: [b0,b1,b2,tau,b3,tau2], extended NS without restrictions
#            Qtc       nx3 matrix, data on [bond prices,time to maturity,coupons]. Use ytm instead
#                      of bond prices if ytmLoss==1
#            s0        scalar, restricted value of short rate
#            ytmLoss   scalar, 0: mimimize squared price errors;
#                      1: minimize squared ytm errors (default 0)
#            weight    scalar or nx1 vector, weights in the loss function
#
#  Output:   Loss    scalar, sum of squared differences between implied and actual bond
#                    prices (or yields)
#
#  Calls on: BondNSxPs, BondYieldToMatPs (if ytmLoss==1)
#
#------------------------------------------------------------------------------

  (Q,tm,c) = [Qtc[:,i] for i=1:3]       #data on bond prices, time to maturity, coupons
  n        = length(c)           #number of bonds

  if length(b) == 3             #standard Nelson-Siegel with restriction b1 = s0-b0
    (b0,b1,b2,tau,b3,tau2) = (abs(b[1]),s0-abs(b[1]),b[2],abs(b[3]),0.0,1.0)
  elseif length(b) == 4         #standard Nelson-Siegel
   (b0,b1,b2,tau,b3,tau2)  = (abs(b[1]),b[2],b[3],abs(b[4]),0.0,1.0)
  elseif length(b) == 5         #extended Nelson-Siegel with restriction b1 = s0 - b0
    (b0,b1,b2,tau,b3,tau2) = (abs(b[1]),s0-abs(b[1]),b[2],abs(b[3]),b[4],abs(b[5]))
  elseif length(b) == 6         #extended Nelson-Siegel
   (b0,b1,b2,tau,b3,tau2)  = (abs(b[1]),b[2],b[3],abs(b[4]),b[5],abs(b[6]))
  end

  QNS = fill(NaN,n)
  for i = 1:n                            #loop over bonds
    ti     = filter(z->z>0,mod(tm[i],1):tm[i])      #time to coupon payments, >0 only
    d      = BondNSxPs(ti,b0,b1,b2,tau,b3,tau2)[3]  #NSx: spot, forward, discount fn
    QNS[i] = sum(d.*c[i]) + d[end]                  #fitted bond price
    if ytmLoss == 1
      QNS[i] = BondYieldToMatPs(QNS[i],c[i],ti,1,1,0.05,1e-7)  #fitted ytm
    end
  end
  Loss = 1.0 + 100*sum( weight.*(QNS - Q).^2 )  #weighted sum of squared deviations of fitted from actual

  return Loss

end
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
# messing with the tree ordering

# using Revise

using Caesar, IncrementalInference
using DistributedFactorGraphs

using Combinatorics
using Suppressor

# import IncrementalInference: resetBuildTreeFromOrder!

filepath = "/media/dehann/temp2/caesar/2019-12-12T00:46:15.438/fg_before_x35"

fg = initfg()
loadDFG(filepath, Main, fg)

drawGraph(fg)

vo = getEliminationOrder(fg)

vo2 = reverse(vo)

tree = resetBuildTreeFromOrder!(fg, vo2)

drawTree(tree, show=true)

### what we want for brute force cost on tree ordering...

# - generate the next ordering
vo = nextInPermutation

# our own
tree = buildTreeFromOrdering!(vo)
getTreeCostBruteForce(tree; revision=1, kargs...)

# comparible to
getIncidenceMatrixCost_MD(vo)
getIncidenceMatrixCost_AMD(incidence wrt vo)
getIncidenceMatrixCost_COLAMD
getIncidenceMatrixCost_CCOLAMD

# using Revise

using Caesar, IncrementalInference
using DistributedFactorGraphs

using Combinatorics
using Suppressor

# import IncrementalInference: resetBuildTreeFromOrder!

filepath = "/media/dehann/temp2/caesar/2019-12-12T00:46:15.438/fg_before_x35"

fg = initfg()
loadDFG(filepath, Main, fg)

function main(fg::AbstractDFG; upto::Int=6, num::Int=100)

  fgs = deepcopy(fg)

  delVars = filter(x->isSolvable(getVariable(fgs, x))==0, ls(fgs))
  todel = setdiff(lsf(fgs, solvable=0), lsf(fgs, solvable=1))
  delFcts = intersect(lsf(fgs), todel)
  allMags = filter(x->:MAGNETOMETER in getTags(getFactor(fgs, x)), lsfPriors(fgs) )
  union!(delFcts, filter(x->length(ls(fgs, x))==0, allMags) )

  union!(delVars, (ls(fgs, r"x\d") |> sortDFG)[upto:end])
  union!(delFcts, map(x->ls(fgs, x), delVars)...)

  map(x->deleteFactor!(fgs, x), delFcts)
  map(x->deleteVariable!(fgs, x), delVars)

  vars = ls(fgs)

  vos = permutations(vars) |> collect
  svos = deepcopy(vos)

  svos = round.(Int, length(vos)*rand(num) )
  TREES = Dict{Int, BayesTree}()
  @suppress begin
    ind = 0
    for vo in vos[svos]
      ind += 1
      TREES[ind] = resetBuildTreeFromOrder!(fgs, vo)
    end
  end

  return TREES
end

@time trees = main(fg)

#
# vos = rand(1000)
# svos = round.(Int, length(vos)*rand(100) )
# vos[svos]

#


text:
module TestGP
using GaussianProcesses, ScikitLearnBase
using Test, Random

Random.seed!(1)

@testset "GP" begin
    d, n = 3, 10

    X = 2π * rand(d, n)
    y = [sum(sin, view(X, :, i)) / d for i in 1:n]
    mZero = MeanZero()
    kern = SE(0.0, 0.0)

    gp = GP(X, y, mZero, kern)

    # Verify that predictive mean at input observations
    # are the same as the output observations
    @testset "Predictive mean" begin
        y_pred, sig = predict_y(gp, X)
        @test maximum(abs, gp.y - y_pred) ≈ 0.0 atol=0.1
    end

    # ScikitLearn interface test
    @testset "ScikitLearn interface" begin
        gp_sk = ScikitLearnBase.fit!(GPE(), X', y)
        y_pred = ScikitLearnBase.predict(gp_sk, X')
        @test maximum(abs, gp_sk.y - y_pred) ≈ 0.0 atol=0.1
    end

    # Modify kernel and update
    @testset "Update" begin
        gp.kernel.ℓ2 = 4.0
        X_pred = 2π * rand(d, n)

        GaussianProcesses.update_target!(gp)
        y_pred, sig = predict_y(gp, X_pred)
    end

    #Check that the rand function works
    @testset "Random GP sampling" begin
        X_test = 2π * rand(d, n)
        samples = rand(gp, X_test)
    end

end
end


text:
#!/usr/bin/env julia

# Script to parse and compare the libdevice PDF manual against our list of intrinsics

function parse_intrinsics(cb)
    fn = joinpath(@__DIR__, "..", "src", "device", "libdevice.jl")
    open(fn) do f
        for ln in eachline(f)
            m = match(r"@wrap ([\w.]+\(.+?\)::\w+)", ln)
            if m != nothing
                cb(replace(m.captures[1], r"\w+::", "::"))
            end
        end
    end
end

function parse_libdevice(fn, cb)
    open(fn) do f
        next_proto = false
        number = 0

        for ln in eachline(f)
            if (m = match(r"^\d\.(\d+)\..", ln); m != nothing)
                number = parse(Int, m.captures[1])
            elseif occursin(r"^Prototype:", ln)
                next_proto = true
            elseif next_proto
                cb(chomp(ln), number)
                next_proto = false
            end
        end
    end
end

function main(args)
    if length(args) != 1
        println("Usage: $(basename(@__FILE__)) LIBDEVICE_PDF")
        exit(1)
    end
    pdf = args[1]
    isfile(pdf) || error("input PDF does not exist")

    wrapped = Set{String}()
    parse_intrinsics(intr -> push!(wrapped, intr))

    intrinsics = Set{String}()
    numbering = Dict{String,Number}()
    txt = tempname()
    run(`pdftotext $pdf $txt`)
    parse_libdevice(txt, (proto, number) -> begin
        m = match(r"^(\w+) (@[\w.]+)\((.*?)\)", proto)
        if m != nothing
            rettype = m.captures[1]
            fn = m.captures[2]
            arglist = m.captures[3]

            argpairs = split(arglist, ", ")
            argtypes, args = zip(map(argpair -> split(argpair, " "), argpairs)...)

            wrap_fn = strip(fn, '@')
            wrap_argtypes = map(argtyp -> endswith(argtyp, '*') ? "Ptr{$(argtyp[1:end-1])}"
                                                                : argtyp, argtypes)
            wrap_args = map(arg -> strip(arg, '%'), args)
            wrap_arglist = join(["$arg::$argtyp" for (arg, argtyp) in zip(wrap_args, wrap_argtypes)], ", ")

            intr = "$wrap_fn($wrap_arglist)::$rettype"
            push!(intrinsics, intr)
            numbering[intr] = number
        end
    end)
    rm(txt)

    missing = setdiff(intrinsics, wrapped)
    superfluous = setdiff(wrapped, intrinsics)

    println("Missing intrinsics:")
    for intr in sort(collect(missing), lt=(a,b)->numbering[a]<numbering[b])
        println(" $(numbering[intr]). $intr")
    end

    println()

    println("Superfluous intrinsics:")
    for intr in sort(collect(superfluous))
        println(" - $intr")
    end
end

main(ARGS)


text:
using Test, Parsers, MaGe

@testset "Hits" begin
    h = Hit(1.1, 2.2, 3.3, 4.4, 5.5, 6, 7, 8)

    @test location(h) == (1.1, 2.2, 3.3)
    @test energy(h) == 4.4
    @test time(h) == 5.5
    @test particleid(h) == 6
    @test trackid(h) == 7
    @test trackparentid(h) == 8

    @test_throws InexactError Hit(1, 2, 3, 4, 5, 1.1, 2.2, 3.3)
end

@testset "Events" begin
    v = [Hit(i*i, sqrt(i), 3.5, 4, 0.5, 1, 2, 3) for i in 1:5]
    e = Event(1, 5, 2, v)

    @test eltype(typeof(e)) == Hit

    @test energy(e) == 4*5

    @test hits(e) == v
    @test eventnum(e) == 1
    @test length(e) == hitcount(e) == 5
    @test primarycount(e) == 2

    # Wrong number of hits for hitcount
    @test_throws ArgumentError Event(1, 3, 2, v)
end

@testset "Parsing .root.hits files" begin
    path = tempname()

    open(path, "w") do io
        truncate(io, 0)
        write(io, "123 456 789\n")
    end

    open(path) do io
        @test MaGe.parse_meta(io) == (123, 456, 789)
    end

    open(path, "w") do io
        truncate(io, 0)
        write(io, "1.2 -2.3 -3.4 4.5 5.6 67891011 121314 1516 physiDet
                  0.11 0.1234321 -333 -0.44 0.001 1 2 3 physiDet")
    end

    open(path) do io
        @test MaGe.parse_hit(io) == Hit(1.2, -2.3, -3.4, 4.5, 5.6, 67891011, 121314, 1516)
        @test MaGe.parse_hit(io) == Hit(0.11, 0.1234321, -333, -0.44, 0.001, 1, 2, 3)
    end

    open(path, "w") do io
        truncate(io, 0)
        write(io, "624 2 3
                  1.60738 -2.07026 -201.594 0.1638 0 22 187 4 physiDet
                  1.91771 -2.52883 -201.842 0.24458 0 22 187 4 physiDet")
    end

    testreader = MaGe.loadstreaming(path)

    e = read(testreader)
    @test eventnum(e) == 624
    @test hitcount(e) == 2
    @test primarycount(e) == 3

    @test hits(e)[1] == Hit(1.60738, -2.07026, -201.594, 0.1638, 0, 22, 187, 4)
    @test hits(e)[2] == Hit(1.91771, -2.52883, -201.842, 0.24458, 0, 22, 187, 4)

    for (h1, h2) in zip(e, hits(e))
        @test h1 === h2
    end

    @test isnothing(iterate(testreader))

    testpath = realpath(joinpath(dirname(pathof(MaGe)), "..", "test", "test.root.hits"))

    @test MaGe.is_root_hit_file(testpath)

    events = MaGe.load(testpath)

    lastevent = events[end]

    @test eventnum(lastevent) == 999851

    @test hitcount(lastevent) == 319

    @test primarycount(lastevent) == 4

    @test hits(lastevent)[1] == Hit(-1.2325, -2.44103, -196.814, 0.07764, 0.0, 22, 9, 6)

    @test hits(lastevent)[end] == Hit(-1.18915, -2.28214, -198.958, 8.4419, 0.0, 11, 165, 16)

    @test MaGe.loadstreaming(testpath) do stream
        for e in stream
            if eventnum(e) == 999851
                return hitcount(e) == 319
            end
        end
        return false
    end

    @test MaGe.loadstreaming(testpath) do stream
        while !eof(stream)
            e = read(stream)
            if eventnum(e) == 999851
                return hitcount(e) == 319
            end
        end
        return false
    end

    @test all(extrema(energy, MaGe.loadstreaming(testpath)) .≈ (17.05186, 3068.0469217))

    badpath = realpath(joinpath(dirname(pathof(MaGe)), "..", "test", "bad.root.hits"))
    badreader = MaGe.loadstreaming(badpath)
    @test_throws Parsers.Error read(badreader)
    @test_throws Parsers.Error read(badreader)
end


text:
function free_energy_nopool(rbm::ConvRBM, v::AbstractArray; β::Real = true)
    E = energy(visible(rbm), v)
    I = inputs_v_to_h(rbm, v)
    F = free_energy(hidden(rbm), I; β)
    @assert size(E) == (vsizes(rbm, v).input_size..., vsizes(rbm, v).batch_size...)
    @assert size(F) == (hsizes(rbm, I).output_size..., vsizes(rbm, v).batch_size...)
    # reduce over kernel dims
    E_conv = reshape_maybe(sum(E; dims=1:kernel_ndims(rbm)), vsizes(rbm, v).batch_size)
    F_conv = reshape_maybe(sum(F; dims=1:kernel_ndims(rbm)), vsizes(rbm, v).batch_size)
    return E_conv + F_conv
end

function free_energy_pooled(rbm::ConvRBM, v::AbstractArray; β::Real = 1)
    E = energy(visible(rbm), v)
    I = inputs_v_to_h(rbm, v)
    F = RBMs.free_energies(hidden(rbm), I; β)
    @assert size(E) == (vsizes(rbm, v).input_size..., batch_size(rbm, v)...)
    @assert size(F) == size(I) == (size(hidden(rbm))..., output_size(rbm, v)..., batch_size(rbm, v)...)
    G = -logsumexp(-β * F; dims=output_dims(rbm)) / β # pooling step <=> only one h active per group
    @assert size(G)[1:ndims(hidden(rbm))] == size(hidden(rbm))
    @assert all(size(G)[output_dims(rbm)] .== 1)
    # reduce over remaining dims
    E_reduced = reshape_maybe(sum(E; dims=1:kernel_ndims(rbm)), batch_size(rbm, v))
    F_reduced = reshape_maybe(sum(G; dims=1:ndims(hidden(rbm))), batch_size(rbm, v))
    return E_reduced + F_reduced
end

function sample_h_from_v_nopool(rbm::ConvRBM, v::AbstractArray; β::Real = true)
    inputs = RBMs.inputs_v_to_h(rbm, v)
    return RBMs.transfer_sample(hidden(rbm), inputs; β)
end

function sample_h_from_v_pooled(rbm::ConvRBM, v::AbstractArray; β::Real = 1)
    inputs = RBMs.inputs_v_to_h(rbm, v)
    F = RBMs.free_energies(hidden(rbm), inputs; β)
    @assert size(F) == size(inputs) == (size(hidden(rbm))..., output_size(rbm, v)..., batch_size(rbm, v)...)
    selected_size = (size(hidden(rbm))..., batch_size(rbm, v)...)
    k = reshape(cartesian_sample_from_logits_gumbel(-β * F; dims=output_dims(rbm)), selected_size)

    selected_inputs = similar(inputs, selected_size)
    for n in CartesianIndices(batch_size(rbm, v)), μ in CartesianIndices(size(hidden(rbm)))
        @inbounds selected_inputs[μ,n] = inputs[k[μ,n]]
    end

    selected_h = RBMs.transfer_sample(hidden(rbm), selected_inputs; β)
    @assert size(selected_h) == size(selected_inputs) == selected_size

    h = fill!(similar(selected_h, size(F)), 0)
    for n in CartesianIndices(batch_size(rbm, v)), μ in CartesianIndices(size(hidden(rbm)))
        h[k[μ,n]] = selected_h[μ,n]
    end
    return h
end