text:
###############################################################################
#
#   fmpz_mod_rel_series.jl : Power series over flint fmpz integers mod n
#
###############################################################################

export fmpz_mod_rel_series, FmpzModRelSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fmpz_mod_rel_series)
   val = pol_length(a) + valuation(a) - 1
   val < 0 && throw(DomainError())
   z = fmpz_mod_rel_series(modulus(a), Array{fmpz}(0), 0, val, val)
   z.parent = parent(a)
   return z
end

elem_type(::Type{FmpzModRelSeriesRing}) = fmpz_mod_rel_series

parent_type(::Type{fmpz_mod_rel_series}) = FmpzModRelSeriesRing

base_ring(R::FmpzModRelSeriesRing) = R.base_ring

var(a::FmpzModRelSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FmpzModRelSeriesRing) = R.prec_max

function normalise(a::fmpz_mod_rel_series, len::Int)
   if len > 0
      c = fmpz()
      ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_mod_rel_series}, Int), c, a, len - 1)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz_mod_rel_series}, Int), c, a, len - 1)
      end
   end
   return len
end

function pol_length(x::fmpz_mod_rel_series)
   return ccall((:fmpz_mod_poly_length, :libflint), Int, (Ref{fmpz_mod_rel_series},), x)
end

precision(x::fmpz_mod_rel_series) = x.prec

function polcoeff(x::fmpz_mod_rel_series, n::Int)
   R = base_ring(x)
   if n < 0
      return R(0)
   end
   z = fmpz()
   ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz_mod_rel_series}, Int), z, x, n)
   return R(z)
end

zero(R::FmpzModRelSeriesRing) = R(0)

one(R::FmpzModRelSeriesRing) = R(1)

function gen(R::FmpzModRelSeriesRing)
   z = fmpz_mod_rel_series(modulus(R), [fmpz(1)], 1, max_precision(R) + 1, 1)
   z.parent = R
   return z
end

modulus(R::FmpzModRelSeriesRing) = modulus(base_ring(R))

function deepcopy_internal(a::fmpz_mod_rel_series, dict::ObjectIdDict)
   z = fmpz_mod_rel_series(a)
   z.prec = a.prec
   z.val = a.val
   z.parent = parent(a)
   return z
end

function renormalize!(z::fmpz_mod_rel_series)
   i = 0
   zlen = pol_length(z)
   zval = valuation(z)
   zprec = precision(z)
   while i < zlen && iszero(polcoeff(z, i))
      i += 1
   end
   z.prec = zprec
   if i == zlen
      z.val = zprec
   else
      z.val = zval + i
      ccall((:fmpz_mod_poly_shift_right, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int), z, z, i)
   end
   return nothing
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FmpzModRelSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fmpz_mod_rel_series}) = show_minus_one(fmpz)

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fmpz_mod_rel_series)
   z = parent(x)()
   ccall((:fmpz_mod_poly_neg, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}),
               z, x)
   z.prec = x.prec
   z.val = x.val
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fmpz_mod_rel_series, b::fmpz_mod_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpz_mod_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpz_mod_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function -(a::fmpz_mod_rel_series, b::fmpz_mod_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   lenz = max(lena, lenb)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpz_mod_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpz_mod_poly_neg, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}), z, z)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpz_mod_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:fmpz_mod_poly_sub_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpz_mod_poly_sub_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function *(a::fmpz_mod_rel_series, b::fmpz_mod_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z = parent(a)()
   z.val = a.val + b.val
   z.prec = prec + z.val
   if lena == 0 || lenb == 0
      return z
   end
   lenz = min(lena + lenb - 1, prec)
   ccall((:fmpz_mod_poly_mullow, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, a, b, lenz)
   renormalize!(z)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Generic.Res{fmpz}, y::fmpz_mod_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fmpz_mod_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz}),
               z, y, x.data)
   renormalize!(z)
   return z
end

*(x::fmpz_mod_rel_series, y::Generic.Res{fmpz}) = y * x

function *(x::fmpz, y::fmpz_mod_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   r = mod(x, modulus(y))
   ccall((:fmpz_mod_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz}),
               z, y, r)
   renormalize!(z)
   return z
end

*(x::fmpz_mod_rel_series, y::fmpz) = y * x

*(x::Integer, y::fmpz_mod_rel_series) = fmpz(x)*y

*(x::fmpz_mod_rel_series, y::Integer) = y * x

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpz_mod_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   z = fmpz_mod_rel_series(x)
   z.prec = x.prec + len
   z.val = x.val + len
   z.parent = parent(x)
   return z
end

function shift_right(x::fmpz_mod_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   xval = valuation(x)
   z = parent(x)()
   if len >= xlen + xval
      z.prec = max(0, x.prec - len)
      z.val = max(0, x.prec - len)
   else
      z.prec = max(0, x.prec - len)
      z.val = max(0, xval - len)
      zlen = min(xlen + xval - len, xlen)
      ccall((:fmpz_mod_poly_shift_right, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, x, xlen - zlen)
      renormalize!(z)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fmpz_mod_rel_series, prec::Int)
   prec < 0 && throw(DomainError())
   xlen = pol_length(x)
   xprec = precision(x)
   xval = valuation(x)
   if xprec + xval <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   if prec <= xval
      z = parent(x)()
      z.val = prec
      z.prec = prec
   else
      z.val = xval
      ccall((:fmpz_mod_poly_set_trunc, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, x, min(prec - xval, xlen))
   end
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpz_mod_rel_series, b::Int)
   b < 0 && throw(DomainError())
   if isgen(a)
      z = parent(a)()
      z = setcoeff!(z, 0, fmpz(1))
      z.prec = a.prec + b - 1
      z.val = b
   elseif pol_length(a) == 0
      z = parent(a)()
      z.prec = b*valuation(a)
      z.val = b*valuation(a)
   elseif pol_length(a) == 1
      z = parent(a)([polcoeff(a, 0)^b], 1,
                           (b - 1)*valuation(a) + precision(a), b*valuation(a))
      renormalize!(z)
      return z
   elseif b == 0
      return one(parent(a))
   else
      z = parent(a)()
      z.prec = a.prec + (b - 1)*valuation(a)
      z.val = b*valuation(a)
      ccall((:fmpz_mod_poly_pow_trunc, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int, Int),
               z, a, b, z.prec - z.val)
   end
   renormalize!(z)
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fmpz_mod_rel_series, y::fmpz_mod_rel_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)
   if prec <= x.val && prec <= y.val
      return true
   end
   if x.val != y.val
      return false
   end
   xlen = normalise(x, min(pol_length(x), prec - x.val))
   ylen = normalise(y, min(pol_length(y), prec - y.val))
   if xlen != ylen
      return false
   end
   return Bool(ccall((:fmpz_mod_poly_equal_trunc, :libflint), Cint,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               x, y, xlen))
end

function isequal(x::fmpz_mod_rel_series, y::fmpz_mod_rel_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || x.val != y.val || pol_length(x) != pol_length(y)
      return false
   end
   return Bool(ccall((:fmpz_mod_poly_equal, :libflint), Cint,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               x, y, pol_length(x)))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpz_mod_rel_series, y::Generic.Res{fmpz})
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         z = fmpz()
         ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
                       (Ref{fmpz}, Ref{fmpz_mod_rel_series}, Int), z, x, 0)
         return ccall((:fmpz_equal, :libflint), Bool,
               (Ref{fmpz}, Ref{fmpz}), z, y.data)
      else
         return false
      end
   else
      return iszero(y.data)
   end
end

==(x::Generic.Res{fmpz}, y::fmpz_mod_rel_series) = y == x

function ==(x::fmpz_mod_rel_series, y::fmpz)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         z = fmpz()
         ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
                       (Ref{fmpz}, Ref{fmpz_mod_rel_series}, Int), z, x, 0)
         r = mod(y, modulus(x))
         return ccall((:fmpz_equal, :libflint), Bool,
               (Ref{fmpz}, Ref{fmpz}, Int), z, r, 0)
      else
         return false
      end
   else
      r = mod(y, modulus(x))
      return iszero(r)
   end
end

==(x::fmpz, y::fmpz_mod_rel_series) = y == x

==(x::fmpz_mod_rel_series, y::Integer) = x == fmpz(y)

==(x::Integer, y::fmpz_mod_rel_series) = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpz_mod_rel_series, y::fmpz_mod_rel_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   yval = valuation(y)
   xval = valuation(x)
   if yval != 0
      if xval >= yval
         x = shift_right(x, yval)
         y = shift_right(y, yval)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec - x.val, y.prec - y.val)
   z = parent(x)()
   z.val = xval - yval
   z.prec = prec + z.val
   if prec != 0
      ccall((:fmpz_mod_poly_div_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, x, y, prec)
   end
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpz_mod_rel_series, y::Generic.Res{fmpz})
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.val = x.val
   ccall((:fmpz_mod_poly_scalar_div_fmpz, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz}),
               z, x, y.data)
   return z
end

function divexact(x::fmpz_mod_rel_series, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   r = mod(y, modulus(x))
   ccall((:fmpz_mod_poly_scalar_div_fmpz, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz}),
               z, x, r)
   return z
end

divexact(x::fmpz_mod_rel_series, y::Integer) = divexact(x, fmpz(y))

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fmpz_mod_rel_series)
   iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ainv.val = 0
   ccall((:fmpz_mod_poly_inv_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               ainv, a, a.prec)
   return ainv
end

###############################################################################
#
#   Special functions
#
###############################################################################

function exp(a::fmpz_mod_rel_series)
   if iszero(a)
      z = one(parent(a))
      z.prec = precision(a)
      z.val = valuation(a)
      return z
   end
   z = parent(a)()
   R = base_ring(a)
   vala = valuation(a)
   preca = precision(a)
   d = Array{fmpz}(preca)
   c = vala == 0 ? polcoeff(a, 0) : R()
   d[1] = exp(c).data
   len = pol_length(a) + vala
   z0 = fmpz()
   for k = 1 : preca - 1
      s = fmpz()
      for j = 1 : min(k + 1, len) - 1
         c = j >= vala ? polcoeff(a, j - vala).data : z0
         s += j * c * d[k - j + 1]
      end
      !isunit(base_ring(a)(k)) && error("Unable to divide in exp")
      d[k + 1] = divexact(base_ring(a)(s), k).data
   end
   z = parent(a)(d, preca, preca, 0)
   ccall((:_fmpz_mod_poly_set_length, :libflint), Void,
         (Ref{fmpz_mod_rel_series}, Int), z, normalise(z, preca))
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(x::fmpz_mod_rel_series)
  ccall((:fmpz_mod_poly_zero, :libflint), Void,
                   (Ref{fmpz_mod_rel_series},), x)
  x.prec = parent(x).prec_max
  return x
end

function fit!(x::fmpz_mod_rel_series, n::Int)
  ccall((:fmpz_mod_poly_fit_length, :libflint), Void,
                   (Ref{fmpz_mod_rel_series}, Int), x, n)
  return nothing
end

function setcoeff!(z::fmpz_mod_rel_series, n::Int, x::fmpz)
   ccall((:fmpz_mod_poly_set_coeff_fmpz, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Int, Ref{fmpz}),
               z, n, x)
   return z
end

function setcoeff!(z::fmpz_mod_rel_series, n::Int, x::Generic.Res{fmpz})
   ccall((:fmpz_mod_poly_set_coeff_fmpz, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Int, Ref{fmpz}),
               z, n, x.data)
   return z
end

function mul!(z::fmpz_mod_rel_series, a::fmpz_mod_rel_series, b::fmpz_mod_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z.val = a.val + b.val
   z.prec = prec + c.val
   lenz = min(lena + lenb - 1, prec)
   if lena <= 0 || lenb <= 0
      lenz = 0
   end
   ccall((:fmpz_mod_poly_mullow, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               z, a, b, lenz)
   renormalize!(z)
   return z
end

function addeq!(a::fmpz_mod_rel_series, b::fmpz_mod_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   modulus = modulus(a)
   if a.val < b.val
      z = fmpz_mod_rel_series(modulus)
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpz_mod_poly_set_trunc, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               a, a, z, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpz_mod_poly_truncate, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Int),
            a, max(0, lenz - a.val + b.val))
      ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
            (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
            a, a, a.val - b.val)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               a, a, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               a, a, b, lenz)
   end
   a.prec = prec
   a.val = val
   renormalize!(a)
   return a
end

function add!(c::fmpz_mod_rel_series, a::fmpz_mod_rel_series, b::fmpz_mod_rel_series)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenc = max(lena, lenb)
   c.prec = prec
   ccall((:fmpz_mod_poly_add_series, :libflint), Void,
                (Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Ref{fmpz_mod_rel_series}, Int),
               c, a, b, lenc)
   return c
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fmpz_mod_rel_series}, ::Type{T}) where {T <: Integer} = fmpz_mod_rel_series

promote_rule(::Type{fmpz_mod_rel_series}, ::Type{fmpz}) = fmpz_mod_rel_series

promote_rule(::Type{fmpz_mod_rel_series}, ::Type{Generic.Res{fmpz}}) = fmpz_mod_rel_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FmpzModRelSeriesRing)()
   z = fmpz_mod_rel_series(modulus(a))
   z.prec = a.prec_max
   z.val = a.prec_max
   z.parent = a
   return z
end

function (a::FmpzModRelSeriesRing)(b::Integer)
   if b == 0
      z = fmpz_mod_rel_series(modulus(a))
      z.prec = a.prec_max
      z.val = a.prec_max
   else
      z = fmpz_mod_rel_series(modulus(a), [fmpz(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpzModRelSeriesRing)(b::fmpz)
   if iszero(b)
      z = fmpz_mod_rel_series(modulus(a))
      z.prec = a.prec_max
      z.val = a.prec_max
   else
      z = fmpz_mod_rel_series(modulus(a), [b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpzModRelSeriesRing)(b::Generic.Res{fmpz})
   if iszero(b)
      z = fmpz_mod_rel_series(modulus(a))
      z.prec = a.prec_max
      z.val = a.prec_max
   else
      z = fmpz_mod_rel_series(modulus(a), [b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpzModRelSeriesRing)(b::fmpz_mod_rel_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FmpzModRelSeriesRing)(b::Array{fmpz, 1}, len::Int, prec::Int, val::Int)
   z = fmpz_mod_rel_series(modulus(a), b, len, prec, val)
   z.parent = a
   return z
end

function (a::FmpzModRelSeriesRing)(b::Array{Generic.Res{fmpz}, 1}, len::Int, prec::Int, val::Int)
   z = fmpz_mod_rel_series(modulus(a), b, len, prec, val)
   z.parent = a
   return z
end


text:
using UnicodePlots, ReferenceTests, Test
using ReferenceTests: BeforeAfterFull
using Dates: Date, Day
using SparseArrays
using StableRNGs
using StatsBase
using Random

import Random: seed!

tests = [
    "tst_issues.jl",
    "tst_common.jl",
    "tst_graphics.jl",
    "tst_canvas.jl",
    "tst_plot.jl",
    "tst_barplot.jl",
    "tst_histogram.jl",
    "tst_scatterplot.jl",
    "tst_lineplot.jl",
    "tst_spy.jl",
    "tst_boxplot.jl",
    "tst_heatmap.jl"
]

const RNG = StableRNG(1337)
const TERM_SIZE = (24, 80)

# see JuliaTesting/ReferenceTests.jl/pull/91
test_ref(reference, actual) = @test_reference(reference, actual, render = BeforeAfterFull(), format = "TXT")

# helpers
macro show_col(p, kv...)
    :(@io2str($(Expr(:call, :show, Expr(:call, :IOContext, :(::IO), :color=>true, kv...), esc(p)))))
end
macro show_nocol(p, kv...)
    :(@io2str($(Expr(:call, :show, Expr(:call, :IOContext, :(::IO), :color=>false, kv...), esc(p)))))
end
macro print_col(p, kv...)
    :(@io2str($(Expr(:call, :print, Expr(:call, :IOContext, :(::IO), :color=>true, kv...), esc(p)))))
end
macro print_nocol(p, kv...)
    :(@io2str($(Expr(:call, :print, Expr(:call, :IOContext, :(::IO), :color=>false, kv...), esc(p)))))
end

# sprand or sprandn is not stable across versions (e.g. 1.0 vs 1.6)
function stable_sprand(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat)
  I = Int[]
  J = Int[]
  for li in randsubseq(r, 1:(m*n), density)
    j, i = divrem(li - 1, m)
    push!(I, i + 1)
    push!(J, j + 1)
  end
  V = rand(r, length(I))
  return sparse(I, J, V)
end

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

ext: jl


text:
"""
    WeightedMean(T = Float64)

Simple weighted mean, tracked as type `T`.

# Example:
    o = fit!(WeightedMean(), rand(100), rand(100))
    sum(o)
    mean(o)
"""
mutable struct WeightedMean{T} <: WeightedOnlineStat{T}
    μ::T
    W::T
    n::Int
    function WeightedMean{T}(μ = T(0), W = T(0), n = 0) where T
        new{T}(T(μ), T(W), Int(n))
    end
end

WeightedMean(μ::T, W::T, n::Int) where T = WeightedMean{T}(μ, W, n)
WeightedMean(::Type{T}) where T = WeightedMean(T(0), T(0), 0)
WeightedMean() = WeightedMean(Float64)

function OnlineStatsBase._fit!(o::WeightedMean{T}, x, w) where T
    xx = convert(T, x)
    ww = convert(T, w)

    o.n += 1
    o.W = smooth(o.W, ww, T(1) / o.n)
    o.μ = smooth(o.μ, xx, ww / (o.W * o.n))

    o
end

function OnlineStatsBase._merge!(o::WeightedMean{T}, o2::WeightedMean) where T
    o2_W = convert(T, o2.W)
    o2_μ = convert(T, o2.μ)

    o.n += o2.n
    o.W = smooth(o.W, o2_W, o2.n / o.n)
    o.μ = smooth(o.μ, o2_μ, (o2_W * o2.n) / (o.W * o.n))

    o
end

OnlineStatsBase.value(o::WeightedMean) = o.μ
Statistics.mean(o::WeightedMean) = value(o)
Base.sum(o::WeightedMean) = mean(o) * weightsum(o)
Base.copy(o::WeightedMean) = WeightedMean(o.μ, o.W, o.n)


text:
# We expects the data to be in npy format, return a dict of {group: items}, each file is a different group
function load_data(path::String,groupcount::Number; prefix::String="", swapDimension::Bool = true)
    groups_dict = Dict()
    for i = 1:groupcount
        arr = npzread(path * prefix * string(i) * ".npy")
        for (index, value) in enumerate(arr)
            if isnan(value)
                arr[index] = 0.0
            end
        end
        groups_dict[i] = swapDimension ? transpose(arr) : arr
    end
    return groups_dict
end

# Preprocessing on the samples, global_preprocessing and local_preprocessing are functions. only same dimesions input output are supported atm
function preprocessing!(samples_dict,local_dim::Number, global_preprocessing, local_preprocessing)
    gp = x->Base.invokelatest(global_preprocessing,x)
    lp = x->Base.invokelatest(local_preprocessing,x)
    if global_preprocessing == nothing && local_preprocessing == nothing
        return
    end
    for (k,v) in samples_dict
        if global_preprocessing != nothing
            samples_dict[k][1 : local_dim-1,:] = mapslices(gp,v[1 : local_dim-1,:], dims= [2])
        end
        if local_preprocessing != nothing
            samples_dict[k][local_dim : end,:] = mapslices(lp,v[local_dim : end,:], dims= [2])
        end
    end
end

function dcolwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    n = length(r)
    for j = 1:n
        v = zero(promote_type(eltype(a), eltype(b)))
        for i = 1:size(a, 1)
            @inbounds v += a[i, j]*b[i, j]
        end
        r[j] = v
    end
end

function distributer_factory(arr::AbstractArray)
    function distributer!(a,b)
        #for arg in args
        try
            show(a)
            show(b)
            arr[b[1]] = b[2]
        catch y
            show(a)
            show(b)
            println("Exception: ", y)
        end
        return -1
        #end

    end
    return distributer!
end

# Note that we expect the log_likelihood_array to be in rows (samples) x columns (clusters) , this is due to making it more efficent that way.
# function sample_log_cat_array!(labels::AbstractArray{Int64,2}, log_likelihood_array::AbstractArray{Float64,2})
#     # println("lsample log cat" * string(log_likelihood_array))
#     max_log_prob_arr = maximum(log_likelihood_array, dims = 2)
#     log_likelihood_array .-= max_log_prob_arr
#     map!(exp,log_likelihood_array,log_likelihood_array)
#     # println("lsample log cat2" * string(log_likelihood_array))
#     #sum_prob_arr = sum(log_likelihood_array, dims =[2])
#     sum_prob_arr = (cumsum(log_likelihood_array, dims =2))
#     randarr = rand(length(labels)) .* sum_prob_arr[:,(size(sum_prob_arr,2))]
#     sum_prob_arr .-= randarr
#     sum_prob_arr[sum_prob_arr .< 0] .= maxintfloat()
#     #replace!(x -> x < 0 ? maxintfloat() : x, sum_prob_arr)
#     labels .= mapslices(argmin, sum_prob_arr, dims= [2])
# end

function sample_log_cat_array!(labels::AbstractArray{Int64,2}, log_likelihood_array::AbstractArray{Float64,2})
    # println("lsample log cat" * string(log_likelihood_array))
    log_likelihood_array[isnan.(log_likelihood_array)] .= -Inf #Numerical errors arent fun
    max_log_prob_arr = maximum(log_likelihood_array, dims = 2)
    log_likelihood_array .-= max_log_prob_arr
    map!(exp,log_likelihood_array,log_likelihood_array)
    # println("lsample log cat2" * string(log_likelihood_array))
    sum_prob_arr = sum(log_likelihood_array, dims =[2])
    log_likelihood_array ./=  sum_prob_arr
    for i=1:length(labels)
        labels[i,1] = sample(1:size(log_likelihood_array,2), ProbabilityWeights(log_likelihood_array[i,:]))
    end
end

function sample_log_cat(logcat_array::AbstractArray{Float64, 1})
    max_logprob::Float64 = maximum(logcat_array)
    for i=1:length(logcat_array)
        logcat_array[i] = exp(logcat_array[i]-max_logprob)
    end
    sum_logprob::Float64 = sum(logcat_array)
    i::Int64 = 1
    c::Float64 = logcat_array[1]
    u::Float64 = rand()*sum(logcat_array)
    while c < u && i < length(logcat_array)
        c += logcat_array[i += 1]
    end
    return i
end

function create_sufficient_statistics(dist::distribution_hyper_params, pts::Array{Any,1})
    return create_sufficient_statistics(dist,dist, Array{Float64}(undef, 0, 0))
end

function get_labels_histogram(labels)
    hist_dict = Dict()
    for v in labels
        if haskey(hist_dict,v) == false
            hist_dict[v] = 0
        end
        hist_dict[v] += 1
    end
    return sort(collect(hist_dict), by=x->x[1])
end

function create_global_labels(group::local_group)
    clusters_dict = Dict()
    for (i,v) in enumerate(group.local_clusters)
        clusters_dict[i] = v.globalCluster
    end
    return [clusters_dict[i] for i in group.labels]
end

function print_global_sub_cluster(group::local_group)
    println([v.globalCluster for v in group.local_clusters])
    println([v.globalCluster_subcluster for v in group.local_clusters])
end

function print_groups_global_clusters(model::hdp_shared_features)
    for (k,g) in model.groups_dict
        println("Group: " * string(k) * "Global Clusters: " * string([x.globalCluster for x in g.local_clusters]))
    end
end

function axes_swapper(groups_pts_dict::Dict, axes_swap_vector)
    for (k,v) in groups_pts_dict
        groups_pts_dict[k] = v[axes_swap_vector,:]
    end
    return groups_pts_dict
end

function create_params_jld(jld_path,
     random_seed,
     data_path,
     data_prefix,
     groups_count,
     global_preprocessing,
     local_preprocessing,
     iterations,
     hard_clustering,
     total_dim,
     local_dim,
     split_stop,
     argmax_sample_stop,
     α,
     γ,
     global_weight,
     local_weight,
     initial_global_clusters,
     initial_local_clusters,
     global_hyper_params,
     local_hyper_params)
     @save jld_path random_seed data_path data_prefix groups_count global_preprocessing local_preprocessing iterations hard_clustering total_dim local_dim split_stop argmax_sample_stop α γ global_weight local_weight initial_global_clusters initial_local_clusters global_hyper_params local_hyper_params
end

function print_params_to_files(file_path,
     random_seed,
     iterations,
     hard_clustering,
     split_stop,
     argmax_sample_stop,
     α,
     γ,
     global_weight,
     local_weight,
     initial_global_clusters,
     initial_local_clusters,
     global_hyper_params,
     local_hyper_params,
     global_multiplier = 0,
     local_multiplier = 0)
    io = open(file_path, "w+")
    println(io, "random_seed = " * string(random_seed))
    println(io, "iterations = " * string(iterations))
    println(io, "hard_clustering = " * string(hard_clustering))
    println(io, "split_stop = " * string(split_stop))
    println(io, "argmax_sample_stop = " * string(argmax_sample_stop))
    println(io, "α = " * string(α))
    println(io, "γ = " * string(γ))
    println(io, "global_weight = " * string(global_weight))
    println(io, "local_weight = " * string(local_weight))
    println(io, "initial_global_clusters = " * string(initial_global_clusters))
    println(io, "initial_local_clusters = " * string(initial_local_clusters))
    println(io, "global_hyper_params = " * string(global_hyper_params))
    println(io, "local_hyper_params = " * string(local_hyper_params))
    println(io, "global_multiplier = " * string(global_multiplier))
    println(io, "local_multiplier = " * string(local_multiplier))
    close(io)
 end

 function get_node_leaders_dict()
     leader_dict = Dict()
     cur_leader = 2
     leader_dict[cur_leader] = []
     for i in workers()
         if i in procs(cur_leader)
             push!(leader_dict[cur_leader], i)
         else
             cur_leader = i
             leader_dict[cur_leader] = [i]
         end
     end
     return leader_dict
 end

function assign_group_leaders(groups_count, leader_dict)
    group_assignments = zeros(groups_count)
    group_leaders = collect(keys(leader_dict))
    for i=1:length(group_assignments)
        group_assignments[i] = group_leaders[i%length(group_leaders) +1 ]
    end
    return group_assignments
end

function log_multivariate_gamma(x::Number, D::Number)
    res::Float64 = D*(D-1)/4*log(pi)
    for d = 1:D
        res += logabsgamma(x+(1-d)/2)[1]
    end
    return res
end


text:
@testset "env_table" begin
    c = raw"""
        \begin{table}{caption="a table", class="table-striped"}
        | A   | B   |
        | --- | --- |
        | 0   | 1   |
        \end{table}""" |> F.fd2html
    @test isapproxstr(c, """

a table

            A
            B

            0
            1

""")
end

@testset "env_column" begin
    s = raw"""
        AB
        \begin{column}{}
        AAA
        \end{column}
        BA
        """ |> F.fd2html
    @test isapproxstr(s, """

AB

 AAA

BA

        """)
end

@testset "env_section" begin
    F.def_LOCAL_VARS!()
    F.def_GLOBAL_VARS!()
    F.LOCAL_VARS["sections"] = Pair{String,String}[] => (Pair{String,String},)
    F.GLOBAL_VARS["section_counter"] = 1 => (Int,)
    counter = F.globvar("section_counter")
    c = F.fd2html(raw"""
        \begin{section}{title="aaa"}
        Hello
        \end{section}
        """, internal=true)
    @test isapproxstr(c, """

aaa

                Hello

        """)
end


text:
# TODO: see julia iterators

"""
Implementation of iterator used in `LogMgrs` module.

Iterates within a `Page` backwards block by block until reaching the 0th block,
at which point I am appointed emperor of China.
"""
module LogIterators
  using HelloDB.Blocks: Block
  using HelloDB.Pages: Page, getint
  using HelloDB: DB
  using HelloDB.Data: INT_SIZE

  ##############################################################################
  # Exports
  ##############################################################################

  export LogIterator
  export hasnext, next

  ##############################################################################
  # Implementation
  ##############################################################################

  type LogIterator
    block::Block
    page::Page
    currentrec::Int

    function LogIterator(block::Block)
      page = Page()
      read(page, block)

      new(block, page, getint(page, DB.logmgr.lastpos))
    end
  end

  function hasnext(itr::LogIterator)
    itr.currentrec > 0 || itr.block.number > 0
  end

  function next(itr::LogIterator)
    if itr.currentrec == 0
      _moveto_next_block(itr)
    end
    itr.currentrec = getint(itr.page, currentrec)
    BasicLogRecord(itr.page, itr.currentrec + INT_SIZE)
  end

  ##############################################################################
  # Private
  ##############################################################################

  function _moveto_next_block(itr)
    block = Block(itr.block.filename, itr.block.number - 1)
    read(itr.page, block)
    itr.currentrec = getint(page, DB.LogMgr.lastpos)
  end
end


text:
using Documenter

push!(LOAD_PATH, "../src/")
using ResizableArrays

DEPLOYDOCS = (get(ENV, "CI", nothing) == "true")

makedocs(
    sitename = "ResizableArrays.jl Package",
    format = Documenter.HTML(
        prettyurls = DEPLOYDOCS,
    ),
    authors = "Éric Thiébaut and contributors",
    pages = ["index.md", "install.md", "usage.md", "library.md"]
)

if DEPLOYDOCS
    deploydocs(
        repo = "github.com/emmt/ResizableArrays.jl.git",
    )
end


text:
using AES, Test

k = [0x2b, 0x7e, 0x15, 0x16, 0x28, 0xae, 0xd2, 0xa6, 0xab, 0xf7, 0x15, 0x88, 0x09, 0xcf, 0x4f, 0x3c]
key = AES128Key(k)

plaintext = "The quick brown fox jumps over the lazy dog."
plaintext8 = "The quick brown fox jumps over the 😴 dog."
iv = [0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f]

cipher = AESCipher(;key_length=128, mode=AES.CBC, key=key)
ct = encrypt(plaintext, cipher; iv=iv)
pt = decrypt(ct, cipher)
ct8 = encrypt(plaintext8, cipher)
pt8 = decrypt(ct8, cipher)

ciphertext = [0xbd, 0x13, 0x20, 0x4f, 0x67, 0xd8, 0x16, 0x7f,
			  0x20, 0x21, 0x1c, 0x99, 0xb0, 0xa7, 0xcc, 0x05,
			  0x06, 0xd5, 0xc7, 0x03, 0xea, 0xfb, 0x01, 0xa7,
			  0xd0, 0x47, 0x3b, 0x5c, 0xc9, 0x99, 0xaa, 0xa2,
			  0x1e, 0x01, 0x7b, 0x76, 0x6e, 0x52, 0x98, 0x5f,
			  0x5e, 0xf7, 0xf0, 0xbd, 0xc5, 0x72, 0xb7, 0x2e]

@testset "ASCII" begin
	@test pt == transcode(UInt8,plaintext)
	@test ct.data == ciphertext
end

@testset "UTF-8" begin
	@test pt8 == transcode(UInt8,plaintext8)
end


text:
abstract type AbstractMPS end

"""
    length(::MPS/MPO)

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

"""
    size::MPS/MPO)

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

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

The Vector storage of an MPS/MPO.

This is 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::Int)
  m.llim = new_ll
end

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

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

function orthocenter(m::T) where {T<:AbstractMPS}
  !isortho(m) && error("$T has no well-defined orthogonality center")
  return leftlim(m)+1
end

Base.getindex(M::AbstractMPS,
              n::Integer) = getindex(data(M),n)

function Base.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 Base.setindex!(M::MPST, v::MPST,
                        ::Colon) where {MPST <: AbstractMPS}
  setleftlim!(M, leftlim(v))
  setrightlim!(M, rightlim(v))
  data(M)[:] = data(v)
  return M
end

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

Base.copy(m::AbstractMPS) = typeof(m)(copy(data(m)),
                                      leftlim(m),
                                      rightlim(m))

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

Base.deepcopy(m::AbstractMPS) = typeof(m)(deepcopy(data(m)),
                                          leftlim(m),
                                          rightlim(m))

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

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

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

"""
    unique_siteind(A::MPO, B::MPS, j::Int)
    unique_siteind(A::MPO, B::MPO, j::Int)

Get the site index of MPO `A` that is unique to `A` (not shared with MPS/MPO `B`).
"""
function unique_siteind(A::AbstractMPS, B::AbstractMPS, j::Int)
  N = length(A)
  j == 1 && return uniqueind(A[j], A[j+1], B[j])
  j == N && return uniqueind(A[j], A[j-1], B[j])
  return uniqueind(A[j], A[j-1], A[j+1], B[j])
end

"""
    unique_siteinds(A::MPO, B::MPS)
    unique_siteinds(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 unique_siteinds(A::AbstractMPS, B::AbstractMPS)
  return [unique_siteind(A, B, j) for j in eachindex(A)]
end

"""
    common_siteind(A::MPO, B::MPS, j::Int)
    common_siteind(A::MPO, B::MPO, j::Int)

Get the site index of MPO `A` that is shared with MPS/MPO `B`.
"""
function common_siteind(A::AbstractMPS, B::AbstractMPS, j::Int)
  return commonind(A[j], B[j])
end

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

Get the site indices of MPO `A` that are shared with MPS/MPO `B`, as a `Vector{<:Index}`.
"""
function common_siteinds(A::AbstractMPS, B::AbstractMPS)
  return [common_siteind(A, B, j) for j in eachindex(A)]
end

Base.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))

"""
    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(ψ::ITensors.AbstractMPS, is) =
  findall(hascommoninds(is), ψ)

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

# TODO: depracate in favor of findsite.
"""
    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`.
"""
function findfirstsiteind(ψ::AbstractMPS,
                          s::Index)
  return findfirst(hasind(s), ψ)
end

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

    findfirstsiteind(M::MPO, is)

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

"""
    firstsiteind(M::Union{MPS,MPO}, j::Int; 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 firstsiteind(M::AbstractMPS, j::Int;
                      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::Int; 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::Int; 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

# TODO: change kwarg from `set_limits` to `preserve_ortho`
function Base.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`
Base.map(f::Function, M::AbstractMPS; set_limits::Bool = true) =
  map!(f, copy(M); set_limits = set_limits)

for fname in (:dag,
              :prime,
              :setprime,
              :noprime,
              :addtags,
              :removetags,
              :replacetags,
              :settags)
  fname_bang = 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.
    """
    $fname(M::AbstractMPS, args...;
           set_limits::Bool = false, kwargs...) =
      map(m -> $fname(m, args...; kwargs...), M;
          set_limits = set_limits)

    """
        $($fname_bang)(M::MPS, args...; kwargs...)

        $($fname_bang)(M::MPO, args...; kwargs...)

    Apply $($fname) to all ITensors of an MPS/MPO in-place.
    """
    $fname_bang(M::AbstractMPS, args...;
                set_limits::Bool = false, kwargs...) =
      map!(m -> $fname(m, args...; kwargs...), M;
           set_limits = set_limits)
  end
end

function map_linkinds!(f::Function, M::AbstractMPS)
  for i in eachindex(M)[1:end-1]
    l = linkind(M, i)
    if !isnothing(l)
      l̃ = f(l)
      M[i, set_limits = false] = replaceind(M[i], l, l̃)
      M[i+1, set_limits = false] = replaceind(M[i+1], l, l̃)
    end
  end
  return M
end

map_linkinds(f::Function, M::AbstractMPS) = map_linkinds!(f, copy(M))

function map_common_siteinds!(f::Function, M1::AbstractMPS,
                                           M2::AbstractMPS)
  length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
  for i in eachindex(M1)
    s = common_siteind(M1, M2, i)
    if !isnothing(s)
      s̃ = f(s)
      M1[i, set_limits = false] = replaceind(M1[i], s, s̃)
      M2[i, set_limits = false] = replaceind(M2[i], s, s̃)
    end
  end
  return M1, M2
end

function map_common_siteinds(f::Function, M1::AbstractMPS,
                                          M2::AbstractMPS)
  return map_common_siteinds!(f, copy(M1), copy(M2))
end

function map_unique_siteinds!(f::Function, M1::AbstractMPS,
                                           M2::AbstractMPS)
  length(M1) != length(M2) && error("MPOs/MPSs must be the same length")
  for i in eachindex(M1)
    s = unique_siteind(M1, M2, i)
    if !isnothing(s)
      s̃ = f(s)
      M1[i, set_limits = false] = replaceind(M1[i], s, s̃)
    end
  end
  return M1
end

function map_unique_siteinds(f::Function, M1::AbstractMPS,
                                        M2::AbstractMPS)
  return map_unique_siteinds!(f, copy(M1), M2)
end

for fname in (:sim,
              :prime,
              :setprime,
              :noprime,
              :addtags,
              :removetags,
              :replacetags,
              :settags)
  fname_linkinds = Symbol(fname, :_linkinds)
  fname_linkinds_inplace = Symbol(fname_linkinds, :!)
  fname_common_siteinds = Symbol(fname, :_common_siteinds)
  fname_common_siteinds_inplace = Symbol(fname_common_siteinds, :!)
  fname_unique_siteinds = Symbol(fname, :_unique_siteinds)
  fname_unique_siteinds_inplace = Symbol(fname_unique_siteinds, :!)

  @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.
    """
    $fname_linkinds(M::AbstractMPS, args...; kwargs...) =
      map_linkinds(i -> $fname(i, args...; kwargs...), M)

    """
        $($fname_linkinds)!(M::MPS, args...; kwargs...)

        $($fname_linkinds)!(M::MPO, args...; kwargs...)

    Apply $($fname) to all link indices of the ITensors of an MPS/MPO in-place.
    """
    $fname_linkinds_inplace(M::AbstractMPS, args...; kwargs...) =
      map_linkinds!(i -> $fname(i, args...; kwargs...), M)

    """
        $($fname_common_siteinds)(M1::MPO, M2::MPS, args...; kwargs...)

        $($fname_common_siteinds)(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_common_siteinds(M1::AbstractMPS,
                                    M2::AbstractMPS,
                                    args...;
                                    kwargs...)
      return map_common_siteinds(i -> $fname(i, args...;
                                             kwargs...), M1, M2)
    end

    """
        $($fname_common_siteinds)!(M1::MPO, M2::MPS, args...; kwargs...)

        $($fname_common_siteinds)!(M1::MPO, M2::MPO, args...; kwargs...)

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

    Modifies the input MPSs/MPOs in-place.
    """
    $fname_common_siteinds_inplace(M1::AbstractMPS, M2::AbstractMPS,
                                   args...; kwargs...) =
      map_common_siteinds!(i -> $fname(i, args...; kwargs...), M1, M2)

    """
        $($fname_unique_siteinds)(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_unique_siteinds(M1::AbstractMPS,
                                    M2::AbstractMPS,
                                    args...;
                                    kwargs...)
      return map_unique_siteinds(i -> $fname(i, args...;
                                             kwargs...), M1, M2)
    end

    """
        $($fname_unique_siteinds)!(M1::MPO, M2::MPS, args...; kwargs...)

    Apply $($fname) to the site indices of `M1` that are not shared with `M2`. Modifies the input MPSs/MPOs in-place.
    """
    function $fname_unique_siteinds_inplace(M1::AbstractMPS,
                                            M2::AbstractMPS,
                                            args...;
                                            kwargs...)
      return map_unique_siteinds!(i -> $fname(i, args...;
                                              kwargs...), M1, M2)
    end
  end
end

"""
    maxlinkdim(M::MPS)

    maxlinkdim(M::MPO)

Get the maximum link dimension of the MPS or MPO.
"""
function maxlinkdim(M::AbstractMPS)
  md = 0
  for b ∈ eachindex(M)[1:end-1]
    l = linkind(M, b)
    linkdim = isnothing(l) ? 0 : dim(l)
    md = max(md, linkdim)
  end
  md
end

"""
    linkind(M::MPS, j::Int)

    linkind(M::MPO, j::Int)

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::Int)
  N = length(M)
  (j ≥ length(M) || j < 1) && return nothing
  return commonind(M[j], M[j+1])
end

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

"""
    linkdim(M::MPS, j::Int)

    linkdim(M::MPO, j::Int)

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::Int)
  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)
  if make_inds_match
    replace_siteinds!(M1dag, siteinds(M2))
  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 LinearAlgebra.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...)

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

"""
    norm(A::MPS)

    norm(A::MPO)

Compute the norm of the MPS or MPO.

See also `lognorm`.
"""
function LinearAlgebra.norm(M::AbstractMPS)
  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` and `loginner`/`logdot`.
"""
function lognorm(M::AbstractMPS)
  return 0.5 * logdot(M, M)
end

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

"""
    add(A::MPS, B::MPS; kwargs...)
    +(A::MPS, B::MPS; kwargs...)

    add(A::MPO, B::MPO; kwargs...)
    +(A::MPO, B::MPO; kwargs...)

Add two MPS/MPO with each other, with some optional
truncation.
"""
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{2}}(undef,N-1)
  second = Vector{ITensor{2}}(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

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 Base.sum(A::Vector{T};
                  kwargs...) where {T <: AbstractMPS}
  length(A) == 0 && return T()
  length(A) == 1 && return A[1]
  length(A) == 2 && return +(A[1], A[2]; kwargs...)
  nsize = isodd(length(A)) ? (div(length(A) - 1, 2) + 1) : div(length(A), 2)
  newterms = Vector{T}(undef, nsize)
  np = 1
  for n in 1:2:length(A) - 1
    newterms[np] = +(A[n], A[n+1]; kwargs...)
    np += 1
  end
  if isodd(length(A))
    newterms[nsize] = A[end]
  end
  return sum(newterms; 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])
    ltags = tags(linkind(M, b))
    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])
    ltags = tags(linkind(M, b))
    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 NDTensors.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])
    U,S,V = svd(M[j],rinds;kwargs...)
    M[j] = U
    M[j-1] *= (S*V)
    setrightlim!(M,j)
  end

end

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

"""
    *(x::Number, M::MPS)

    *(x::Number, M::MPO)

Scale the MPS or MPO by the provided number.

Note: right now this just naively scales the
middle tensor. In the future, the plan would be
to scale the tensors between the left limit
and right limit by `x^(1/N)` where `N` is the distance
between the left limit and right limit.
"""
function Base.:*(x::Number, M::AbstractMPS)
  N = deepcopy(M)
  c = div(length(N), 2)
  N[c] .*= x
  return N
end

Base.:-(M::AbstractMPS) = Base.:*(-1,M)

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

Sets a contiguous range of MPS/MPO tensors
"""
function Base.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

# 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 Base.setindex!(ψ::MPST,
                        A::ITensor,
                        r::UnitRange{Int};
                        orthocenter::Int = 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)
    sites = sites[[perm...]]
  end

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

  ψ[firstsite:lastsite] = ψA

  return ψ
end

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

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; <keyword arguments>)

    MPO(A::ITensor, sites; <keyword arguments>)

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

# Arguments
- `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::Int = 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::Int = 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..., 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

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

Swap the sites `b` and `b+1`.
"""
function swapbondsites(ψ::AbstractMPS, b::Int; 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::Int = 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

# 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)
  for i in 1:N
    ψ = movesite(ψ, ns[i] => ns′[i]; kwargs...)
    ns = _movesite(ns, ns[i] => ns′[i])
  end
  return ψ
end

# TODO: make a permutesites(::MPS/MPO, perm)
# function that that a permutation of the sites
# p(1:N) for N sites
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

"""
    hasqns(M::MPS)

    hasqns(M::MPO)

Return true if the MPS or MPO has
tensors which carry quantum numbers.
"""
hasqns(M::AbstractMPS) = hasqns(M[1])

"""
    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=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

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

@deprecate orthoCenter(args...;
                       kwargs...) orthocenter(args...; kwargs...)

import .NDTensors.store

@deprecate store(m::AbstractMPS) data(m)

@deprecate replacesites!(args...;
                         kwargs...) ITensors.replace_siteinds!(args...; kwargs...)

@deprecate applyMPO(args...; kwargs...) contract(args...; kwargs...)

@deprecate applympo(args...; kwargs...) contract(args...; kwargs...)

@deprecate errorMPOprod(args...;
                        kwargs...) error_contract(args...;
                                                  kwargs...)

@deprecate error_mpoprod(args...;
                         kwargs...) error_contract(args...;
                                                   kwargs...)

@deprecate error_mul(args...;
                     kwargs...) error_contract(args...;
                                               kwargs...)

@deprecate multMPO(args...; kwargs...) contract(args...; kwargs...)

import Base.sum

@deprecate sum(A::AbstractMPS,
               B::AbstractMPS; kwargs...) add(A, B; kwargs...)

@deprecate multmpo(args...;
                   kwargs...) contract(args...; kwargs...)

@deprecate set_leftlim!(args...;
                        kwargs...) ITensors.setleftlim!(args...;
                                                        kwargs...)

@deprecate set_rightlim!(args...;
                         kwargs...) ITensors.setrightlim!(args...;
                                                          kwargs...)

@deprecate tensors(args...;
                   kwargs...) ITensors.data(args...; kwargs...)

@deprecate primelinks!(args...;
                       kwargs...) ITensors.prime_linkinds!(args...;
                                                          kwargs...)

@deprecate simlinks!(args...;
                     kwargs...) ITensors.sim_linkinds!(args...;
                                                      kwargs...)

@deprecate mul(A::AbstractMPS,
               B::AbstractMPS;
               kwargs...) contract(A, B; kwargs...)


text:
using MPIMeasurements
using Unitful

# Create Isel BaseScanner object
iR = IselRobot("/dev/ttyS0")

initRefZYX(iR)
################# Use case 1 Basic Move #######################################

# Move absolue LowLevel no safetyCheck!
movePark(iR)
moveCenter(iR)
getPos(iR)
moveAbs(iR, 1.0Unitful.mm, 0.0Unitful.mm, 0.0Unitful.mm)

# Move absolute MidLevel
posXYZ = [1.0Unitful.mm,0.0Unitful.mm,0.0Unitful.mm]
moveAbs(iR, posXYZ)

################# Use case 2 Move and Measure #################################


text:
# ---
# title: Figure tweaking
# cover: assets/figure_tweaking.png
# description: Setting figure attributes
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

# To tweak one or more axes, simply use the `axis` keyword when plotting. For example

df = (x=rand(100), y=rand(100), z=rand(100), c=rand(["a", "b"], 100))
xyc = data(df) * mapping(:x, :y, layout=:c)
layers = linear() + mapping(color=:z)
plt = xyc * layers
fg = draw(plt, axis=(aspect=1,), figure=(resolution=(800, 400),))

#

fg = draw(plt, axis=(aspect=1, xticks=0:0.1:1, yticks=0:0.1:1, ylabel="custom label"))

# save cover image #src
mkpath("assets") #src
save("assets/figure_tweaking.png", fg) #src