text:
export iszero_row, howell_form, kernel_basis, isdiagonal, diagonal

import Nemo.matrix

import Base.vcat

import LinearAlgebra
LinearAlgebra.dot(a::RingElem, b::RingElem) = a*b

################################################################################
#
#  Dense matrix types
#
################################################################################

dense_matrix_type(::Type{T}) where {T} = Generic.MatSpaceElem{T}

# TODO (easy): Check if the following can be simplified to
# coefficient_type(::Type{<:Mat{T}}} where {T} = T
coefficient_type(::Type{fmpz_mat}) = fmpz

coefficient_type(::Type{fmpq_mat}) = fmpq

coefficient_type(::Type{nmod_mat}) = nmod

coefficient_type(::Type{fq_nmod_mat}) = fq_nmod

coefficient_type(::Type{fq_mat}) = fq

coefficient_type(::Type{arb_mat}) = arb

coefficient_type(::Type{acb_mat}) = acb

coefficient_type(::Type{gfp_mat}) = gfp_elem

coefficient_type(::Type{Generic.Mat{T}}) where {T} = T

################################################################################
#
#  Unsafe functions for generic matrices
#
################################################################################

#function zero!(a::MatElem)
#  for i in 1:nrows(a)
#    for j in 1:ncols(a)
#      a[i, j] = zero!(a[i, j])
#    end
#  end
#  return a
#end

function mul!(c::MatElem, a::MatElem, b::MatElem)
  ncols(a) != nrows(b) && error("Incompatible matrix dimensions")
  nrows(c) != nrows(a) && error("Incompatible matrix dimensions")
  ncols(c) != ncols(b) && error("Incompatible matrix dimensions")

  if c === a || c === b
    d = parent(a)()
    return mul!(d, a, b)
  end

  t = base_ring(a)()
  for i = 1:nrows(a)
    for j = 1:ncols(b)
      c[i, j] = zero!(c[i, j])
      for k = 1:ncols(a)
        c[i, j] = addmul_delayed_reduction!(c[i, j], a[i, k], b[k, j], t)
      end
      c[i, j] = reduce!(c[i, j])
    end
  end
  return c
end

function add!(c::MatElem, a::MatElem, b::MatElem)
  parent(a) != parent(b) && error("Parents don't match.")
  parent(c) != parent(b) && error("Parents don't match.")
  for i = 1:nrows(c)
    for j = 1:ncols(c)
      c[i, j] = add!(c[i, j], a[i, j], b[i, j])
    end
  end
  return c
end

################################################################################
#
#  Saturation
#
################################################################################

@doc Markdown.doc"""
    saturate(A::fmpz_mat) -> fmpz_mat

Computes the saturation of $A$, that is, a basis for $\mathbf{Q}\otimes M \meet
\mathbf{Z}^n$, where $M$ is the row span of $A$ and $n$ the number of rows of
$A$.

Equivalently, return $TA$ for an invertiable rational matrix $T$ such that $TA$
is integral and the elementary divisors of $TA$ are all trivial.
"""
function saturate(A::fmpz_mat)
  #row saturation: want
  #  TA in Z, T in Q and elem_div TA = [1]
  #
  #  AT = H (in HNF), then A = HT^-1 and H^-1A = T^-1
  # since T is uni-mod, H^-1 A is in Z with triv. elm. div

  H = hnf(A')
  H = H'
  Hi, d = pseudo_inv(sub(H, 1:nrows(H), 1:nrows(H)))
  S = Hi*A
  Sd = divexact(S, d)
#  @hassert :HNF 1  d*Sd == S
  return Sd
end

################################################################################
#
#  Zero matrix constructors
#
################################################################################

function zero_matrix(::Type{MatElem}, R::Ring, n::Int)
  return zero_matrix(R, n)
end

function zero_matrix(::Type{MatElem}, R::Ring, n::Int, m::Int)
  return zero_matrix(R, n, m)
end

function matrix(A::Array{fmpz, 2})
  m = matrix(FlintZZ, A)
  return m
end

function matrix(A::Array{T, 2}) where T <: RingElem
  r, c = size(A)
  (r < 0 || c < 0) && error("Array must be non-empty")
  m = matrix(parent(A[1, 1]), A)
  return m
end

function matrix(A::Array{T, 1}) where T <: RingElem
  return matrix(reshape(A,length(A),1))
end

function scalar_matrix(R::Ring, n::Int, a::RingElement)
  b = R(a)
  z = zero_matrix(R, n, n)
  for i in 1:n
    z[i, i] = b
  end
  return z
end

function Array(a::fmpz_mat; S::Type{T} = fmpz) where T
  A = Array{T}(undef, nrows(a), ncols(a))
  for i = 1:nrows(a)
    for j = 1:ncols(a)
      A[i,j] = T(a[i,j])
    end
  end
  return A
end

function iszero_row(M::fmpz_mat, i::Int)
  for j = 1:ncols(M)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function iszero_row(M::nmod_mat, i::Int)
  zero = UInt(0)
  for j in 1:ncols(M)
    t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ref{nmod_mat}, Int, Int), M, i - 1, j - 1)
    if t != zero
      return false
    end
  end
  return true
end

function iszero_row(M::MatElem{T}, i::Int) where T
  for j in 1:ncols(M)
    if !iszero(M[i,j])
      return false
    end
  end
  return true
end

function iszero_row(M::Array{T, 2}, i::Int) where T <: Integer
  for j = 1:Base.size(M, 2)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function iszero_row(M::Array{fmpz, 2}, i::Int)
  for j = 1:Base.size(M, 2)
    if M[i,j] != 0
      return false
    end
  end
  return true
end

function iszero_row(M::Array{T, 2}, i::Int) where T <: RingElem
  for j in 1:Base.size(M, 2)
    if !iszero(M[i,j])
      return false
    end
  end
  return true
end

function divexact!(a::fmpz_mat, b::fmpz_mat, d::fmpz)
  ccall((:fmpz_mat_scalar_divexact_fmpz, :libflint), Nothing,
               (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), a, a, d)
end

function mul!(a::fmpz_mat, b::fmpz_mat, c::fmpz)
  ccall((:fmpz_mat_scalar_mul_fmpz, :libflint), Nothing,
                  (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), a, b, c)
end

function _hnf(x::fmpz_mat, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = hnf(reverse_cols(x))
    reverse_cols!(h)
    reverse_rows!(h)
    return h::fmpz_mat
  end
  return hnf(x)::fmpz_mat
end

function _hnf!(x::fmpz_mat, shape::Symbol = :upperright)
  if shape == :lowerleft
    reverse_cols!(x)
    hnf!(x)
    reverse_cols!(x)
    reverse_rows!(x)
    return x::fmpz_mat
  end
  hnf!(x)
  return x::fmpz_mat
end

function hnf!(x::fmpz_mat)
  ccall((:fmpz_mat_hnf, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}), x, x)
  return x
end

function _hnf_modular_eldiv(x::fmpz_mat, m::fmpz, shape::Symbol = :upperright)
  if shape == :lowerleft
    h = hnf_modular_eldiv!(reverse_cols(x), m)
    reverse_cols!(h)
    reverse_rows!(h)
    return h
  elseif shape == :upperright
    return hnf_modular_eldiv(x, m)
  else
    error("Shape $shape not supported")
  end
end

function hnf_modular_eldiv!(x::fmpz_mat, d::fmpz, shape::Symbol = :upperright)
   (nrows(x) < ncols(x)) &&
                error("Matrix must have at least as many rows as columns")
   if shape == :upperright
     ccall((:fmpz_mat_hnf_modular_eldiv, :libflint), Nothing,
                  (Ref{fmpz_mat}, Ref{fmpz}), x, d)
     return x
   elseif shape == :lowerleft
     reverse_cols!(x)
     ccall((:fmpz_mat_hnf_modular_eldiv, :libflint), Nothing,
                 (Ref{fmpz_mat}, Ref{fmpz}), x, d)
     reverse_cols!(x)
     reverse_rows!(x)
     return x
   else
     error("shape $shape is not supported")
   end
end

function ishnf(x::fmpz_mat, shape::Symbol)
  if shape == :upperright
    return ishnf(x)
  elseif shape == :lowerleft
    r = nrows(x)
    i = 0
    j_old = ncols(x) + 1

    for outer i in nrows(x):-1:1

      if iszero_row(x, i)
        break
      end

      j = ncols(x)
      while iszero(x[i, j])
        j = j - 1
      end
      x[i, j] < 0 && return false
      j >= j_old && return false
      for k in i+1:r
        x[k, j] < 0 && return false
        x[k, j] >= x[i, j] && return false
      end
      j_old = j
      i = i - 1
    end

    for l in i:-1:1
      !iszero_row(x, l) && return false
    end
    return true
  end
end

################################################################################
#
#  Is LLL?
#
################################################################################

function islll_reduced(x::fmpz_mat, ctx::lll_ctx = lll_ctx(0.99, 0.51))
  b = ccall((:fmpz_lll_is_reduced_d, :libflint), Cint,
            (Ref{fmpz_mat}, Ref{lll_ctx}), x, ctx)
  return Bool(b)
end

################################################################################
#
################################################################################

function maximum(f::typeof(abs), a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, 0,0)
  for i=1:nrows(a)
    for j=1:ncols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, i-1, j-1)
      if ccall((:fmpz_cmpabs, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) < 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_abs, :libflint), Nothing, (Ref{fmpz}, Ptr{fmpz}), r, m)
  return r
end

function maximum(a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, 0,0)
  for i=1:nrows(a)
    for j=1:ncols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, i-1, j-1)
      if ccall((:fmpz_cmp, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) < 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_set, :libflint), Nothing, (Ref{fmpz}, Ptr{fmpz}), r, m)
  return r
end

function minimum(a::fmpz_mat)
  m = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, 0,0)
  for i=1:nrows(a)
    for j=1:ncols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), a, i-1, j-1)
      if ccall((:fmpz_cmp, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}), m, z) > 0
        m = z
      end
    end
  end
  r = fmpz()
  ccall((:fmpz_set, :libflint), Nothing, (Ref{fmpz}, Ptr{fmpz}), r, m)
  return r
end

################################################################################
#
#  Reduce mod hnf
#
################################################################################

# H must be in upper right HNF form
function reduce_mod_hnf_ur!(a::fmpz_mat, H::fmpz_mat)
  for c = 1:nrows(a)
    j = 1
    for i=1:min(nrows(H), ncols(H))
      while j <= ncols(H) && iszero(H[i, j])
        j += 1
      end
      if j > ncols(H)
        break
      end
      q = fdiv(a[c, j], H[i, j])
      if iszero(q)
        continue
      end
      for k=j:ncols(a)
        a[c, k] = a[c, k] - q * H[i, k]
      end
    end
  end
  return nothing
end

function reduce_mod_hnf_ll!(a::fmpz_mat, H::fmpz_mat)
  for c = 1:nrows(a)
    j = ncols(H)
    for i = nrows(H):-1:1
      while j > 0 && iszero(H[i, j])
        j -= 1
      end
      if iszero(j)
        break
      end
      q = fdiv(a[c, j], H[i, j])
      if iszero(q)
        continue
      end
      for k = 1:j
        a[c, k] = a[c, k] - q * H[i, k]
      end
    end
  end
  return nothing
end

################################################################################
#
#  Lift of matrices to overrings
#
################################################################################

@doc Markdown.doc"""
    lift(a::Generic.Mat{Generic.Res{fmpz}}) -> fmpz_mat

It returns a lift of the matrix to the integers.
"""
function lift(a::Generic.Mat{Generic.Res{fmpz}})
  z = zero_matrix(FlintZZ, nrows(a), ncols(a))
  for i in 1:nrows(a)
    for j in 1:ncols(a)
      z[i, j] = lift(a[i, j])
    end
  end
  return z
end

if Nemo.version() > v"0.15.1"
  function lift(a::Generic.Mat{Nemo.fmpz_mod})
    z = zero_matrix(FlintZZ, nrows(a), ncols(a))
    for i in 1:nrows(a)
      for j in 1:ncols(a)
        z[i, j] = lift(a[i, j])
      end
    end
    return z
  end
end

function lift_unsigned(a::nmod_mat)
  z = zero_matrix(FlintZZ, nrows(a), ncols(a))
  ccall((:fmpz_mat_set_nmod_mat_unsigned, :libflint), Nothing,
          (Ref{fmpz_mat}, Ref{nmod_mat}), z, a)
  return z
end

################################################################################
#
#  Kernel function
#
################################################################################
@doc Markdown.doc"""
    kernel_basis(a::MatElem{T}, side:: Symbol) -> Vector{Vector{T}} where {T <: AbstractAlgebra.FieldElem}

It returns a basis for the kernel of the matrix defined over a field. If side is $:right$ or not
specified, the right kernel is computed. If side is $:left$, the left kernel is computed.
"""
function kernel_basis(A::MatElem{T}, side::Symbol = :right) where T<: AbstractAlgebra.FieldElem
  if side == :right
    return right_kernel_basis(A)
  elseif side == :left
    return left_kernel_basis(A)
  else
    error("Unsupported argument: :$side for side: Must be :left or :right")
  end
end

@doc Markdown.doc"""
    right_kernel_basis(a::MatElem{T}) -> Vector{Vector{T}} where {T <: AbstractAlgebra.FieldElem}

It returns a basis for the right kernel of the matrix defined over a field
"""
function right_kernel_basis(a::MatElem{T}) where T <: AbstractAlgebra.FieldElem
  R = base_ring(a)
  n, z = nullspace(a)
  ar = typeof(Array{T}(undef, nrows(z)))[]
  for i in 1:n
    t = Array{T}(undef, nrows(z))
    for j in 1:nrows(z)
      t[j] = R(z[j, i])
    end
    push!(ar,t)
  end
  return ar
end

@doc Markdown.doc"""
    left_kernel_basis(a::MatElem{T}) -> Vector{Vector{T}}

It returns a basis for the left kernel of the matrix
"""
left_kernel_basis(a::MatElem{T}) where T <: AbstractAlgebra.FieldElem = right_kernel_basis(transpose(a))

@doc Markdown.doc"""
    kernel(a::MatElem{T}; side::Symbol = :right) -> Int, MatElem{T}

It returns a tuple $(n, M)$, where n is the rank of the kernel and $M$ is a basis for it. If side is $:right$ or not
specified, the right kernel is computed. If side is $:left$, the left kernel is computed.
"""
function kernel(A::MatElem; side::Symbol = :right)
  if side == :right
    return right_kernel(A)
  elseif side == :left
    return left_kernel(A)
  else
    error("Unsupported argument: :$side for side: Must be :left or :right")
  end
end

@doc Markdown.doc"""
    right_kernel(a::gfp_mat) ->  Int, gfp_mat

It returns a tuple (n, M) where M is a matrix whose columns generate
the kernel and n is the dimension of the kernel.
"""
function right_kernel(x::gfp_mat)
  z = zero_matrix(base_ring(x), ncols(x), max(nrows(x),ncols(x)))
  n = ccall((:nmod_mat_nullspace, :libflint), Int, (Ref{gfp_mat}, Ref{gfp_mat}), z, x)
  return n, z
end

function left_kernel(x::gfp_mat)
  n, M = right_kernel(transpose(x))
  return n, transpose(M)
end

@doc Markdown.doc"""
    left_kernel(a::fmpz_mat) -> Int, fmpz_mat

It returns a tuple (n, M) where M is a matrix whose rows generate
the kernel of M and n is the rank of the kernel.
"""
function left_kernel(x::fmpz_mat)
  H, U = hnf_with_transform(x)
  i = 1
  for outer i in 1:nrows(H)
    if iszero_row(H, i)
      break
    end
  end
  if iszero_row(H, i)
    return nrows(U)-i+1, view(U, i:nrows(U), 1:ncols(U))
  else
    return 0, zero_matrix(FlintZZ, 0, ncols(U))
  end
end

right_kernel(M::MatElem) = nullspace(M)

function left_kernel(M::MatElem)
  rk, M1 = nullspace(transpose(M))
  return rk, transpose(M1)
end

@doc Markdown.doc"""
    right_kernel(a::fmpz_mat) -> Int, fmpz_mat

It returns a tuple (n, M) where M is a matrix whose rows generate
the kernel of M and n is the rank of the kernel.
"""
function right_kernel(x::fmpz_mat)
  n, M = left_kernel(transpose(x))
  return n, transpose(M)
end

@doc Markdown.doc"""
    right_kernel(a::nmod_mat) -> Int, nmod_mat

It returns a tuple (n, M) where M is a matrix whose rows generate
the kernel of M and n is the rank of the kernel.
"""
function right_kernel(M::nmod_mat)
  R = base_ring(M)
  if isprime(modulus(R))
    k = zero_matrix(R, ncols(M), ncols(M))
    n = ccall((:nmod_mat_nullspace, :libflint), Int, (Ref{nmod_mat}, Ref{nmod_mat}), k, M)
    return n, k
  end

  H = hcat(M', identity_matrix(R, ncols(M)))
  if nrows(H) < ncols(H)
    H = vcat(H, zero_matrix(R, ncols(H) - nrows(H), ncols(H)))
  end
  howell_form!(H)
  nr = 1
  while nr <= nrows(H) && !iszero_row(H, nr)
    nr += 1
  end
  nr -= 1
  h = sub(H, 1:nr, 1:nrows(M))
  for i=1:nrows(h)
    if iszero_row(h, i)
      k = sub(H, i:nrows(h), nrows(M)+1:ncols(H))
      return nrows(k), k'
    end
  end
  return 0, zero_matrix(R,nrows(M),0)
end

function left_kernel(a::nmod_mat)
  n, M = right_kernel(transpose(a))
  return n, transpose(M)
end

if Nemo.version() > v"0.15.1"
  function right_kernel(M::fmpz_mod_mat)
    R = base_ring(M)
    N = hcat(M', identity_matrix(R, ncols(M)))
    if nrows(N) < ncols(N)
      N = vcat(N, zero_matrix(R, ncols(N) - nrows(N), ncols(N)))
    end
    howell_form!(N)
    H = N
    nr = 1
    while nr <= nrows(H) && !iszero_row(H, nr)
      nr += 1
    end
    nr -= 1
    h = sub(H, 1:nr, 1:nrows(M))
    for i=1:nrows(h)
      if iszero_row(h, i)
        k = sub(H, i:nrows(h), nrows(M)+1:ncols(H))
        return nrows(k), k'
      end
    end
    return 0, zero_matrix(R,nrows(M),0)
  end

  function left_kernel(a::fmpz_mod_mat)
    n, M = right_kernel(transpose(a))
    return n, transpose(M)
  end
else
  function right_kernel(M::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
    R = base_ring(M)
    N = hcat(M', identity_matrix(R, ncols(M)))
    if nrows(N) < ncols(N)
      N = vcat(N, zero_matrix(R, ncols(N) - nrows(N), ncols(N)))
    end
    H = howell_form(N)
    nr = 1
    while nr <= nrows(H) && !iszero_row(H, nr)
      nr += 1
    end
    nr -= 1
    h = sub(H, 1:nr, 1:nrows(M))
    for i=1:nrows(h)
      if iszero_row(h, i)
        k = sub(H, i:nrows(h), nrows(M)+1:ncols(H))
        return nrows(k), k'
      end
    end
    return 0, zero_matrix(R,nrows(M),0)
  end

  function left_kernel(a::Generic.Mat{Nemo.Generic.Res{Nemo.fmpz}})
    n, M = right_kernel(transpose(a))
    return n, transpose(M)
  end
end

################################################################################
#
#  Kernel over different rings
#
################################################################################

@doc Markdown.doc"""
    kernel(a::MatrixElem{T}, R::Ring; side::Symbol = :right) -> n, MatElem{elem_type(R)}

It returns a tuple $(n, M)$, where n is the rank of the kernel over $R$ and $M$ is a basis for it. If side is $:right$ or not
specified, the right kernel is computed. If side is $:left$, the left kernel is computed.
"""
function kernel(M::MatrixElem, R::Ring; side::Symbol = :right)
  MP = change_base_ring(R, M)
  return kernel(MP, side = side)
end

################################################################################
#
#  Diagonal (block) matrix creation
#
################################################################################

@doc Markdown.doc"""
    diagonal_matrix(x::Vector{T}) where T <: RingElem -> MatElem{T}

Returns a diagonal matrix whose diagonal entries are the element of $x$.
"""
function diagonal_matrix(x::Vector{T}) where T <: RingElem
  M = zero_matrix(parent(x[1]), length(x), length(x))
  for i = 1:length(x)
    M[i, i] = x[i]
  end
  return M
end

function diagonal_matrix(x::T...) where T <: RingElem
  return diagonal_matrix(collect(x))
end

@doc Markdown.doc"""
    diagonal_matrix(x::Vector{T}) where T <: MatElem -> MatElem

Returns a block diagonal matrix whose diagonal blocks are the matrices in $x$.
"""
function diagonal_matrix(x::Vector{T}) where T <: MatElem
  return cat(x..., dims = (1, 2))
end

function diagonal_matrix(x::T...) where T <: MatElem
  return cat(x..., dims = (1, 2))
end
################################################################################
#
#  Copy matrix into another matrix
#
################################################################################

# Copy B into A at position (i, j)
function _copy_matrix_into_matrix(A::fmpz_mat, i::Int, j::Int, B::fmpz_mat)
  for k in 0:nrows(B) - 1
    for l in 0:ncols(B) - 1
      d = ccall((:fmpz_mat_entry, :libflint),
                Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), B, k, l)
      t = ccall((:fmpz_mat_entry, :libflint),
                Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), A, i - 1 + k, j - 1 + l)
      ccall((:fmpz_set, :libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}), t, d)
    end
  end
end

@doc Markdown.doc"""
    isposdef(a::fmpz_mat) -> Bool

Tests if $a$ positive definite by testing if all principal minors
have positive determinant.
"""
function isposdef(a::fmpz_mat)
  for i=1:nrows(a)
    if det(sub(a, 1:i, 1:i)) <= 0
      return false
    end
  end
  return true
end

#scales the i-th column of a by 2^d[1,i]
function mult_by_2pow_diag!(a::Array{BigFloat, 2}, d::fmpz_mat)
  s = size(a)
  R = RealRing()::_RealRing
  tmp_mpz::BigInt = R.z1
  for i = 1:s[1]
    for j = 1:s[2]
      e = ccall((:mpfr_get_z_2exp, :libmpfr), Clong, (Ref{BigInt}, Ref{BigFloat}), tmp_mpz, a[i,j])
      ccall((:mpfr_set_z_2exp, :libmpfr), Nothing, (Ref{BigFloat}, Ref{BigInt}, Clong, Int32), a[i,j], tmp_mpz, e+Clong(Int(d[1,j])), __get_rounding_mode())
    end
  end
end

#converts BigFloat -> fmpz via round(a*2^l), in a clever(?) way
function round_scale(a::Array{BigFloat, 2}, l::Int)
  s = size(a)
  b = zero_matrix(FlintZZ, s[1], s[2])
  return round_scale!(b, a, l)
end

function round_scale!(b::fmpz_mat, a::Array{BigFloat, 2}, l::Int)
  s = size(a)
  R = RealRing()::_RealRing

  local tmp_mpz::BigInt, tmp_fmpz::fmpz
  tmp_mpz = R.z1
  tmp_fmpz = R.zz1
  tmp_mpfr = deepcopy(a[1,1])  #cannot use the R.?? tmp variable as it may/will
                               #have the wrong precision

  rd = __get_rounding_mode()
  for i = 1:s[1]
    for j = 1:s[2]
      e = a[i,j].exp
      a[i,j].exp += l
      ccall((:mpfr_round, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), tmp_mpfr, a[i,j], rd)
      a[i,j].exp = e
      f = ccall((:mpfr_get_z_2exp, :libmpfr), Clong, (Ref{BigInt}, Ref{BigFloat}),
        tmp_mpz, tmp_mpfr)
      ccall((:fmpz_set_mpz, :libflint), Nothing, (Ref{fmpz}, Ref{BigInt}), tmp_fmpz, tmp_mpz)
      if f > 0
        ccall((:fmpz_mul_2exp, :libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, UInt), tmp_fmpz, tmp_fmpz, f)
      else
        ccall((:fmpz_tdiv_q_2exp, :libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, UInt), tmp_fmpz, tmp_fmpz, -f);
      end
      setindex!(b, tmp_fmpz, i, j)
    end
  end
  return b
end

function round_scale!(b::fmpz_mat, a::arb_mat, l::Int)
  s = size(a)

  R = base_ring(a)
  r = R()
  for i = 1:s[1]
    for j = 1:s[2]
      v = ccall((:arb_mat_entry_ptr, :libarb), Ptr{arb},
                    (Ref{arb_mat}, Int, Int), a, i - 1, j - 1)
      ccall((:arb_mul_2exp_si, :libarb), Nothing, (Ref{arb}, Ptr{arb}, Int), r, v, l)
      b[i,j] = round(fmpz, r)
    end
  end
  return b
end

function shift!(g::fmpz_mat, l::Int)
  for i=1:nrows(g)
    for j=1:ncols(g)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), g, i-1, j-1)
      if l > 0
        ccall((:fmpz_mul_2exp, :libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, l)
      else
        ccall((:fmpz_tdiv_q_2exp, :libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Int), z, z, -l)
      end
    end
  end
  return g
end

################################################################################
#
#  Reduce the entries of a matrix modulo p
#
################################################################################

@doc Markdown.doc"""
    mod!(M::fmpz_mat, p::fmpz)
Reduces every entry modulo $p$ in-place, ie. applies the mod function to every entry.
Positive residue system.
"""
function mod!(M::fmpz_mat, p::fmpz)
  for i=1:nrows(M)
    for j=1:ncols(M)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), M, i - 1, j - 1)
      ccall((:fmpz_mod, :libflint), Nothing, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), z, z, p)
    end
  end
  nothing
end

@doc Markdown.doc"""
    mod(M::fmpz_mat, p::fmpz) -> fmpz_mat
Reduces every entry modulo $p$, ie. applies the mod function to every entry.
"""
function mod(M::fmpz_mat, p::fmpz)
  N = deepcopy(M)
  mod!(N, p)
  return N
end

@doc Markdown.doc"""
    mod!(M::fmpz_mat, p::fmpz)
Reduces every entry modulo $p$ in-place, into the symmetric residue system.
"""
function mod_sym!(M::fmpz_mat, B::fmpz)
  @assert !iszero(B)
  ccall((:fmpz_mat_scalar_smod, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), M, M, B)
end
mod_sym!(M::fmpz_mat, B::Integer) = mod_sym!(M, fmpz(B))

@doc Markdown.doc"""
    mod(M::fmpz_mat, p::fmpz) -> fmpz_mat
Reduces every entry modulo $p$ into the symmetric residue system.
"""
function mod_sym(M::fmpz_mat, B::fmpz)
  N = zero_matrix(FlintZZ, nrows(M), ncols(M))
  ccall((:fmpz_mat_scalar_smod, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), N, M, B)
  return N
end
mod_sym(M::fmpz_mat, B::Integer) = mod_sym(M, fmpz(B))
################################################################################
#
#  Concatenation of matrices
#
################################################################################

@doc Markdown.doc"""
    vcat(A::Array{Generic.Mat, 1}) -> Generic.Mat
    vcat(A::Array{fmpz_mat}, 1}) -> fmpz_mat
Forms a big matrix my vertically concatenating the matrices in $A$.
All component matrices need to have the same number of columns.
"""
function vcat(A::Array{T, 1})  where {S <: RingElem, T <: MatElem{S}}
  if any(x->ncols(x) != ncols(A[1]), A)
    error("Matrices must have same number of columns")
  end
  M = zero_matrix(base_ring(A[1]), sum(nrows, A), ncols(A[1]))
  s = 0
  for i=A
    for j=1:nrows(i)
      for k=1:ncols(i)
        M[s+j, k] = i[j,k]
      end
    end
    s += nrows(i)
  end
  return M
end

function vcat(A::Array{fmpz_mat, 1})
  if any(x->ncols(x) != ncols(A[1]), A)
    error("Matrices must have same number of columns")
  end
  M = zero_matrix(base_ring(A[1]), sum(nrows, A), ncols(A[1]))
  s = 0
  for i=A
    for j=1:nrows(i)
      for k=1:ncols(i)
        M[s+j, k] = i[j,k]
      end
    end
    s += nrows(i)
  end
  return M
end

function vcat(A::Array{nmod_mat, 1})
  if any(x->ncols(x) != ncols(A[1]), A)
    error("Matrices must have same number of columns")
  end
  M = zero_matrix(base_ring(A[1]), sum(nrows, A), ncols(A[1]))
  s = 0
  for i=A
    for j=1:nrows(i)
      for k=1:ncols(i)
        M[s+j, k] = i[j,k]
      end
    end
    s += nrows(i)
  end
  return M
end

function Base.vcat(A::MatElem...)
  r = nrows(A[1])
  c = ncols(A[1])
  R = base_ring(A[1])
  for i=2:length(A)
    @assert ncols(A[i]) == c
    @assert base_ring(A[i]) == R
    r += nrows(A[i])
  end
  X = zero_matrix(R, r, c)
  o = 1
  for i=1:length(A)
    for j=1:nrows(A[i])
      X[o, :] = A[i][j, :]
      o += 1
    end
  end
  return X
end

function Base.hcat(A::Array{T, 1}) where {S <: RingElem, T <: MatElem{S}}
  if any(x->nrows(x) != nrows(A[1]), A)
    error("Matrices must have same number of rows")
  end
  M = zero_matrix(base_ring(A[1]), nrows(A[1]), sum(ncols, A))
  s = 0
  for i = A
    for j=1:ncols(i)
      for k=1:nrows(i)
        M[k, s + j] = i[k,j]
      end
    end
    s += ncols(i)
  end
  return M
end

function Base.hcat(A::MatElem...)
  r = nrows(A[1])
  c = ncols(A[1])
  R = base_ring(A[1])
  for i=2:length(A)
    @assert nrows(A[i]) == r
    @assert base_ring(A[i]) == R
    c += ncols(A[i])
  end
  X = zero_matrix(R, r, c)
  o = 1
  for i=1:length(A)
    for j=1:ncols(A[i])
      X[:, o] = A[i][:, j]
      o += 1
    end
  end
  return X
end

function Base.cat(A::MatElem...;dims)
  @assert dims == (1,2) || isa(dims, Int)

  if isa(dims, Int)
    if dims == 1
      return hcat(A...)
    elseif dims == 2
      return vcat(A...)
    else
      error("dims must be 1, 2, or (1,2)")
    end
  end

  local X
  for i=1:length(A)
    if i==1
      X = hcat(A[1], zero_matrix(base_ring(A[1]), nrows(A[1]), sum(Int[ncols(A[j]) for j=2:length(A)])))
    else
      X = vcat(X, hcat(zero_matrix(base_ring(A[1]), nrows(A[i]), sum(ncols(A[j]) for j=1:i-1)), A[i], zero_matrix(base_ring(A[1]), nrows(A[i]), sum(Int[ncols(A[j]) for j=i+1:length(A)]))))
    end
  end
  return X
end

function Base.hvcat(rows::Tuple{Vararg{Int}}, A::MatElem...)
  B = hcat([A[i] for i=1:rows[1]]...)
  o = rows[1]
  for j=2:length(rows)
    C = hcat([A[i+o] for i=1:rows[j]]...)
    o += rows[j]
    B = vcat(B, C)
  end
  return B
end

################################################################################
#
#  Smith normal form with trafo
#
################################################################################

#=
g, e,f = gcdx(a, b)
U = [1 0 ; -divexact(b, g)*f 1]*[1 1; 0 1];
V = [e -divexact(b, g) ; f divexact(a, g)];

then U*[ a 0; 0 b] * V = [g 0 ; 0 l]
=#
@doc Markdown.doc"""
    snf_with_transform(A::fmpz_mat, l::Bool = true, r::Bool = true) -> fmpz_mat, fmpz_mat, fmpz_mat

Given some integer matrix A, compute the Smith normal form (elementary
divisor normal form) of A. If l and/ or r are true, then the corresponding
left and/ or right transformation matrices are computed as well.
"""
function snf_with_transform(A::fmpz_mat, l::Bool = true, r::Bool = true)
  if r
    R = identity_matrix(FlintZZ, ncols(A))
  end

  if l
    L = identity_matrix(FlintZZ, nrows(A))
  end
  # TODO: if only one trafo is required, start with the HNF that does not
  #       compute the trafo
  #       Rationale: most of the work is on the 1st HNF..
  #
  S = deepcopy(A)
  while !isdiagonal(S)
    if l
      S, T = hnf_with_transform(S)
      L = T*L
    else
      S = hnf(S)
    end

    if isdiagonal(S)
      break
    end
    if r
      S, T = hnf_with_transform(S')
      R = T*R
    else
      S = hnf(S')
    end
    S = S'
  end
  #this is probably not really optimal...
  for i=1:min(nrows(S), ncols(S))
    if S[i,i] == 1
      continue
    end
    for j=i+1:min(nrows(S), ncols(S))
      if S[j,j] == 0
        continue
      end
      if S[i,i] != 0 && S[j,j] % S[i,i] == 0
        continue
      end
      g, e,f = gcdx(S[i,i], S[j,j])
      a = divexact(S[i,i], g)
      S[i,i] = g
      b = divexact(S[j,j], g)
      S[j,j] *= a
      if l
        # U = [1 0; -b*f 1] * [ 1 1; 0 1] = [1 1; -b*f -b*f+1]
        # so row i and j of L will be transformed. We do it naively
        # those 2x2 transformations of 2 rows should be a c-primitive
        # or at least a Nemo/Hecke primitive
        for k=1:ncols(L)
          x = -b*f
#          L[i,k], L[j,k] = L[i,k]+L[j,k], x*L[i,k]+(x+1)*L[j,k]
          L[i,k], L[j,k] = L[i,k]+L[j,k], x*(L[i,k]+L[j,k])+L[j,k]
        end
      end
      if r
        # V = [e -b ; f a];
        # so col i and j of R will be transformed. We do it naively
        # careful: at this point, R is still transposed
        for k=1:nrows(R)
          R[i, k], R[j, k] = e*R[i,k]+f*R[j,k], -b*R[i,k]+a*R[j,k]
        end
      end
    end
  end

  if l
    if r
      return S, L, R'
    else
      # last is dummy
      return S, L, L
    end
  elseif r
    # second is dummy
    return S, R, R'
  else
    # last two are dummy
    return S, S, S
  end
end

################################################################################
#
#  IsUpper\Lower triangular
#
################################################################################

function isupper_triangular(M::MatElem)
  n = nrows(M)
  @assert n == ncols(M)
  for i = 2:n
    for j = 1:i-1
      if !iszero(M[i, j])
        return false
      end
    end
  end
  return true
end

function islower_triangular(M::MatElem)
  n = nrows(M)
  @assert n == ncols(M)
  for i = 1:n
    for j = i+1:n
      if !iszero(M[i, j])
        return false
      end
    end
  end
  return true
end

################################################################################
#
#  Is diagonal
#
################################################################################

@doc Markdown.doc"""
    isdiagonal(A::Mat)

Tests if A is diagonal.
"""
function isdiagonal(A::MatElem)
  for i = 1:ncols(A)
    for j = 1:nrows(A)
      if i != j && !iszero(A[j, i])
        return false
      end
    end
  end
  return true
end

################################################################################
#
#  Diagonal
#
################################################################################

@doc Markdown.doc"""
    diagonal(A::Mat{T}) -> Vector{T}

Returns the diagonal of `A` is an array.
"""
diagonal(A::Generic.Mat{T}) where {T} = T[A[i, i] for i in 1:nrows(A)]

################################################################################
#
#  Triangular solving
#
################################################################################

# Solves A x = b for A upper triangular m\times n matrix and b m\times 1.

@doc Markdown.doc"""
    solve_ut(A::MatElem{T}, b::MatElem{T}) -> MatElem{T})

Given an upper triangular $m \times m$ matrix $A$ and a matrix $b$ of size $m
\times k$, this function computes $x$ such that $Ax = b$. Might fail if
the pivots of $A$ are not invertible.
"""
function solve_ut(A::MatElem{T}, b::MatElem{T}) where T
  m = nrows(A)
  n = ncols(A)
  s = ncols(b)
  @assert m == nrows(b)
  @assert m <= n
  x = zero_matrix(base_ring(A), n, s)
  pivot_cols = Vector{Int}()
  r = 0
  last_pivot = n + 1
  for i = m:-1:1
    for j = 1:last_pivot - 1
      if iszero(A[i, j])
        continue
      end
      for z = 1:s
        x[j, z] = b[i, z]
        for k = 1:r
          if !iszero(A[i, pivot_cols[k]])
            x[j, z] -= A[i, pivot_cols[k]]*x[pivot_cols[k], z]
          end
        end
        q, re = divrem(x[j, z], A[i, j])
        @assert iszero(re)
        x[j, z] = q
      end
      last_pivot = j
      r += 1
      push!(pivot_cols, j)
      break
    end
  end
  return x
end

@doc Markdown.doc"""
    solve_lt(A::MatElem{T}, b::MatElem{T}) -> MatElem{T})

Given a lower triangular $m \times m$ matrix $A$ and a matrix $b$ of size
$m \times k$, this function computes $x$ such that $Ax = b$. Might fail if
that the pivots of $A$ are not invertible.
"""
function solve_lt(A::MatElem{T}, b::MatElem{T}) where T
  m = nrows(A)
  n = ncols(A)
  s = ncols(b)
  @assert m == nrows(b)
  @assert m <= n
  x = zero_matrix(base_ring(A), n, s)
  pivot_cols = Vector{Int}()
  r = 0
  last_pivot = 0
  for i = 1:m
    j = n
    while iszero(A[i, j])
      j -= 1
    end
    for z = 1:s
      x[j, z] = b[i, z]
      for k = 1:r
        if !iszero(A[i, pivot_cols[k]]) && !iszero(x[pivot_cols[k], z])
          x[j, z] -= A[i, pivot_cols[k]]*x[pivot_cols[k], z]
        end
      end
      q, re = divrem(x[j, z], A[i, j])
      @assert iszero(re)
      x[j, z] = q
    end
    last_pivot = j
    r += 1
    push!(pivot_cols, j)
  end
  return x
end

# =======================================
# Array interface for MatElem
# =======================================
#
# TODO: re-write for special types (fmpz_mat e.g.) to gain efficiency
#

length(A::MatElem) = nrows(A) * ncols(A)

Base.ndims(A::MatElem) = 2

function Base.size(A::MatElem, n::Int)
  if n == 1
    return nrows(A)
  elseif n == 2
    return ncols(A)
  elseif n < 1
    error("arraysize: dimension out of range")
  else
    return 1
  end
end

function Base.axes(A::MatElem)
  return (Base.OneTo(nrows(A)), Base.OneTo(ncols(A)))
end

function Base.axes(A::MatElem, n::Int)
  return Base.OneTo(size(A, n))
end

function Base.eachindex(A::MatElem)
  return Base.OneTo(length(A))
end

function Base.stride(A::MatElem, n::Int)
  if n <= 1
    return 1
  elseif n == 2
    return nrows(A)
  else
    return length(A)
  end
end

Base.eltype(A::MatElem{T}) where T <: RingElem = T

getindex(A::MatElem, n::Int) = A[1 + ((n-1) % nrows(A)), 1 + div((n-1), nrows(A))]

function setindex!(A::MatElem{T}, n::Int, s::T) where T <: RingElem
  A[1 + ((n-1) % nrows(A)), 1 + div((n-1), nrows(A))] = s
end

function Base.iterate(A::MatElem, state::Int = 0)
  s = size(A)
  if state < s[1]*s[2]
    state += 1
    return A[state], state
  end
  return nothing
end

Base.IteratorSize(M::MatElem) = Base.HasLength()
Base.IteratorEltype(M::MatElem) = Base.HasEltype()
Base.eltype(M::MatElem) = elem_type(base_ring(M))

if Nemo.version() <= v"0.15.1"
  function setindex!(A::MatElem{T}, b::MatElem{T}, ::Colon, i::Int) where T
    @assert ncols(b) == 1 && nrows(b) == nrows(A)
    for j=1:nrows(A)
      A[j,i] = b[j]
    end
    b
  end

  function setindex!(A::MatElem{T}, b::MatElem{T}, i::Int, ::Colon) where T
    @assert nrows(b) == 1 && ncols(b) == ncols(A)
    for j=1:ncols(A)
      A[i,j] = b[j]
    end
    b
  end
end

function setindex!(A::MatElem, b::Array{<: Any, 1}, ::Colon, i::Int)
  @assert length(b) == nrows(A)
  for j=1:nrows(A)
    A[j,i] = b[j]
  end
  b
end

function setindex!(A::MatElem, b::Array{ <: Any, 1}, i::Int, ::Colon)
  @assert length(b) == ncols(A)
  for j=1:ncols(A)
    A[i,j] = b[j]
  end
  b
end

function setindex!(A::MatElem, b, ::Colon, i::Int)
  for j=1:nrows(A)
    A[j,i] = b
  end
  b
end

function setindex!(A::MatElem, b, i::Int, ::Colon)
  for j=1:ncols(A)
    A[i,j] = b
  end
  b
end

@doc Markdown.doc"""
    reduce_mod!(A::MatElem{T}, B::MatElem{T}) where T <: FieldElem

For a reduced row echelon matrix $B$, reduce $A$ modulo $B$, ie. all the pivot
columns will be zero afterwards.
"""
function reduce_mod!(A::MatElem{T}, B::MatElem{T}) where T <: FieldElem
  if isrref(B)
    scale = false
  else
    scale = true
  end

  for h=1:nrows(A)
    j = 1
    for i=1:nrows(B)
      while iszero(B[i, j])
        j += 1
      end
      if scale
        A[h, :] -= A[h, j] * (inv(B[i, j]) * B[i, :])
      else
        A[h, :] -= A[h, j] * B[i, :]
      end
    end
  end
  return A
end

@doc Markdown.doc"""
    reduce_mod(A::MatElem{T}, B::MatElem{T}) where T <: FieldElem -> MatElem

For a reduced row echelon matrix $B$, reduce $A$ modulo $B$, ie. all the pivot
columns will be zero afterwards.
"""
function reduce_mod(A::MatElem{T}, B::MatElem{T}) where T <: FieldElem
  C = deepcopy(A)
  reduce_mod!(C, B)
  return C
end

@doc Markdown.doc"""
    find_pivot(A::MatElem{<:RingElem}) -> Array{Int, 1}

Find the pivot-columns of the reduced row echelon matrix $A$
"""
function find_pivot(A::MatElem{<:RingElem})
  @assert isrref(A)
  p = Int[]
  j = 0
  for i=1:nrows(A)
    j += 1
    if j > ncols(A)
      return p
    end
    while iszero(A[i,j])
      j += 1
      if j > ncols(A)
        return p
      end
    end
    push!(p, j)
  end
  return p
end

@doc Markdown.doc"""
    can_solve(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: FieldElem -> Bool, MatElem

Tries to solve $Ax = B$ for $x$ if `side = :right` and $xA = B$ if `side =
:left`.
"""
function can_solve(A::MatElem{T}, B::MatElem{T};
                                  side = :right) where T <: FieldElem
  @assert base_ring(A) == base_ring(B)

  if side === :right
    @assert nrows(A) == nrows(B)
    return _can_solve(A, B)
  elseif side === :left
    @assert ncols(A) == ncols(B)
    b, C = _can_solve(transpose(A), transpose(B))
    if b
      return true, transpose(C)
    else
      return false, C
    end
  else
    error("Unsupported argument :$side for side: Must be :left or :right")
  end
end

function _can_solve(A::MatElem{T}, B::MatElem{T}) where T <: FieldElem
  R = base_ring(A)
  mu = [A B]
  rk, mu = rref(mu)
  p = find_pivot(mu)
  if any(i -> i > ncols(A), p)
    return false, B
  end
  sol = zero_matrix(R, ncols(A), ncols(B))
  for i = 1:length(p)
    for j = 1:ncols(B)
      sol[p[i], j] = mu[i, ncols(A) + j]
    end
  end
  return true, sol
end

@doc Markdown.doc"""
    can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: FieldElem -> Bool, MatElem, MatElem

Tries to solve $Ax = B$ for $x$ if `side = :right` or $Ax = B$ if `side = :left`.
It returns the solution and the right respectively left kernel of $A$.
"""
function can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: FieldElem
  @assert base_ring(A) == base_ring(B)
  if side === :right
    @assert nrows(A) == nrows(B)
    return _can_solve_with_kernel(A, B)
  elseif side === :left
    b, C, K = _can_solve_with_kernel(transpose(A), transpose(B))
    @assert ncols(A) == ncols(B)
    if b
      return b, transpose(C), transpose(K)
    else
      return b, C, K
    end
  else
    error("Unsupported argument :$side for side: Must be :left or :right")
  end
end

function _can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}) where T <: FieldElem
  R = base_ring(A)
  mu = [A B]
  rk, mu = rref(mu)
  p = find_pivot(mu)
  if any(i->i>ncols(A), p)
    return false, B, B
  end
  sol = zero_matrix(R, ncols(A), ncols(B))
  for i = 1:length(p)
    for j = 1:ncols(B)
      sol[p[i], j] = mu[i, ncols(A) + j]
    end
  end
  n = zero_matrix(R, ncols(A), ncols(A) - length(p))
  np = sort(setdiff(1:ncols(A), p))
  i = 0
  push!(p, ncols(A)+1)
  for j = 1:length(np)
    if np[j] >= p[i+1]
      i += 1
    end
    if i > 0
      n[p[i], j] = -mu[i, np[j]]
    end
    n[np[j], j] = 1
  end
  return true, sol, n
end

@doc Markdown.doc"""
    can_solve(A::MatElem{T}, B::MatElem{T}, side = :right) where T <: RingElem -> Bool, MatElem

Tries to solve $Ax = B$ for $x$ if `side = :right` or $Ax = B$ if `side = :left`
over a euclidean ring.
"""
function can_solve(A::MatElem{T}, B::MatElem{T};
                                  side = :right) where T <: RingElem
  @assert base_ring(A) == base_ring(B)

  if side === :right
    @assert nrows(A) == nrows(B)
    return _can_solve(A, B)
  elseif side === :left
    @assert ncols(A) == ncols(B)
    b, C = _can_solve(transpose(A), transpose(B))
    if b
      return true, transpose(C)
    else
      return false, C
    end
  else
    error("Unsupported argument :$side for side: Must be :left or :right")
  end
end

function _can_solve(a::MatElem{S}, b::MatElem{S}, side = :left) where S <: RingElem
  H, T = hnf_with_transform(transpose(a))
  b = deepcopy(b)
  z = similar(a, ncols(b), ncols(a))
  l = min(nrows(a), ncols(a))
  for i = 1:ncols(b)
    for j = 1:l
      k = 1
      while k <= ncols(H) && iszero(H[j, k])
        k += 1
      end
      if k > ncols(H)
        continue
      end
      q, r = divrem(b[k, i], H[j, k])
      if !iszero(r)
        return false, b
      end
      for h = k:ncols(H)
        b[h, i] -= q*H[j, h]
      end
      z[i, j] = q
    end
  end
  if !iszero(b)
    return false, b
  end
  return true, transpose(z*T)
end

@doc Markdown.doc"""
    can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}) where T <: RingElem -> Bool, MatElem, MatElem

Tries to solve $Ax = B$ for $x$ if `side = :right` or $xA = B$ if `side = :left`.
It returns the solution and the right respectively left kernel of $A$.
"""
function can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: RingElem
  @assert base_ring(A) == base_ring(B)
  if side === :right
    @assert nrows(A) == nrows(B)
    return _can_solve_with_kernel(A, B)
  elseif side === :left
    b, C, K = _can_solve_with_kernel(transpose(A), transpose(B))
    @assert ncols(A) == ncols(B)
    if b
      return b, transpose(C), transpose(K)
    else
      return b, C, K
    end
  else
    error("Unsupported argument :$side for side: Must be :left or :right")
  end
end

function _can_solve_with_kernel(a::MatElem{S}, b::MatElem{S}) where S <: RingElem
  H, T = hnf_with_transform(transpose(a))
  z = similar(a, ncols(b), ncols(a))
  l = min(nrows(a), ncols(a))
  b = deepcopy(b)
  for i=1:ncols(b)
    for j=1:l
      k = 1
      while k <= ncols(H) && iszero(H[j, k])
        k += 1
      end
      if k > ncols(H)
        continue
      end
      q, r = divrem(b[k, i], H[j, k])
      if !iszero(r)
        return false, b, b
      end
      for h=k:ncols(H)
        b[h, i] -= q*H[j, h]
      end
      z[i, j] = q
    end
  end
  if !iszero(b)
    return false, b, b
  end

  for i = nrows(H):-1:1
    for j = 1:ncols(H)
      if !iszero(H[i,j])
        N = similar(a, ncols(a), nrows(H) - i)
        for k = 1:nrows(N)
          for l = 1:ncols(N)
            N[k,l] = T[nrows(T) - l + 1, k]
          end
        end
        return true, transpose(z*T), N
      end
    end
  end
  N =  similar(a, ncols(a), 0)

  return true, transpose(z*T), N
end

################################################################################
#
#  Elementary divisors
#
################################################################################

@doc Markdown.doc"""
    elementary_divisors(A::fmpz_mat) -> Vector{fmpz}

The elementary divisors of $A$, that is, the diagonal entries of the Smith
normal form of $A$.
"""
function elementary_divisors(A::fmpz_mat)
  s = snf(A)
  return fmpz[s[i,i] for i = 1:min(ncols(s), nrows(s))]
end

@doc Markdown.doc"""
    maximal_elementary_divisors(A::fmpz_mat) -> fmpz

The largest elementary divisor of $A$, that is, the last diagonal entry of the
Smith normal form of $A$.
"""
function maximal_elementary_divisor(A::fmpz_mat)
  return elementary_divisors(A)[end]
end

@doc Markdown.doc"""
    divisors(A::fmpz_mat, d::fmpz) -> fmpz

Returns the diagonal entries of a diagonal form of A. We assume that all the elementary
divisors are divisors of d.
"""
function divisors(M::fmpz_mat, d::fmpz)
  sp = fmpz[2, 3, 5, 7, 11, 13, 17, 19, 23]
  l = fmpz[]
  for p in sp
    c, d = ppio(d, p)
    if !isone(c)
      push!(l, p)
    end
  end
  d = ispower(d)[2]
  M1 = _hnf_modular_eldiv(M, d)
  while !isdiagonal(M1)
    M1 = M1'
    hnf_modular_eldiv!(M1, d)
  end

  for j = 1:nrows(M1)
    if !isone(M1[j,j])
      push!(l, M1[j, j])
    end
  end
  return l
end

@doc Markdown.doc"""
    largest_elementary_divisor(A::fmpz_mat) -> fmpz

The largest elementary divisor of $A$, that is, the last diagonal entry of the
Smith normal form of $A$.
"""
largest_elementary_divisor(A::fmpz_mat) = maximal_elementary_divisor(A)

################################################################################
#
#  Function to convert a matrix to array
#
################################################################################

function to_array(M::fmpq_mat)
  A = Vector{fmpq}(undef, ncols(M)*nrows(M))
  for i = 1:nrows(M)
    for j = 1:ncols(M)
      A[(i-1)*ncols(M) + j] = M[i, j]
    end
  end
  return A
end

################################################################################
#
#  Minpoly and Charpoly
#
################################################################################

function minpoly(M::MatElem)
  k = base_ring(M)
  kx, x = PolynomialRing(k, cached = false)
  return minpoly(kx, M)
end

function charpoly(M::MatElem)
  k = base_ring(M)
  kx, x = PolynomialRing(k, cached = false)
  return charpoly(kx, M)
end


text:
using Combinatorics: combinations

function build_excursion(arity, n_way, disallow, seed = nothing)
    test_list = [zeros(eltype(arity), length(arity), 1)]
    for way_idx in 1:n_way
        push!(test_list, all_combinations(arity .- 1, way_idx))
    end
    test_cases = hcat(test_list...)
    test_cases .+= 1

    if seed !== nothing
        seeded = unique(hcat(seed, test_cases), dims = 2)
    else
        seeded = test_cases
    end

    allowed = zeros(Int, length(seeded))
    allow_cnt = 0
    for case_idx in 1:size(seeded, 2)
        if !disallow(seeded[:, case_idx])
            allow_cnt += 1
            allowed[allow_cnt] = case_idx
        end
    end
    seeded[:, allowed[1:allow_cnt]]
end

function build_excursion_multi(arity, n_way, levels, disallow, seed = nothing)
    param_cnt = length(arity)
    levels[n_way] = [collect(1:length(arity))]
    # excursion_set = Set{Tuple{Vararg{Int64, N} where N}}()
    excursion_set = Set{Array{Int, 1}}()
    for level in keys(levels)
        for indices in levels[level]
            # combos = Set(collect(tuple(x...) for x in combinations(indices, level)))
            for subset_idx in 1:level
                combos = Set(collect(combinations(indices, subset_idx)))
                union!(excursion_set, combos)
            end
        end
    end
    essity = arity .- 1  # 1 is the origin of the excursion
    excursions = sort(collect(excursion_set), by = x -> (length(x), x))
    include_origin = 1
    test_cnt = sum(prod(essity[inds]) for inds in excursions) + include_origin

    test_cases = ones(Int, param_cnt, test_cnt)
    test_cases[:, 1] .= 1
    test_idx = include_origin
    for exc_indices in excursions
        walks = ones(Int, length(exc_indices))
        for sub_idx in 1:prod(essity[exc_indices])
            test_idx += 1
            test_cases[exc_indices, test_idx] = walks .+ 1
            next_multiplicative!(walks, essity[exc_indices])
        end
        @assert walks == ones(Int, length(exc_indices))
    end
    @assert test_idx == test_cnt

    if seed !== nothing
        seeded = unique(hcat(seed, test_cases), dims = 2)
    else
        seeded = test_cases
    end

    allowed = zeros(Int, length(seeded))
    allow_cnt = 0
    for case_idx in 1:size(seeded, 2)
        if !disallow(seeded[:, case_idx])
            allow_cnt += 1
            allowed[allow_cnt] = case_idx
        end
    end
    seeded[:, allowed[1:allow_cnt]]
end

ext: jl


text:
# This file was generated, do not modify it. # hide
PCA = @load PCA pkg=MultivariateStats

pca_mdl = PCA(pratio=1)
pca = machine(pca_mdl, X)
fit!(pca)
PCA
W = MLJ.transform(pca, X);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using Test, StatsBase, CUDA, FixedEffects, PooledArrays, CategoricalArrays

p1 = repeat(1:5, inner = 2)
p2 = repeat(1:5, outer = 2)
x = [ 0.5548445405298847 , 0.9444014472663531 , 0.0510866660400604 , 0.9415750229576445 , 0.697755708534771  , 0.9664962514198971 , 0.12752269572311858, 0.4633531422366297 , 0.03341608526498096, 0.1647934493047556]
fes = [FixedEffect(p1), FixedEffect(p2)]
r_ols =  [-0.2015993617092453,  0.2015993617092464, -0.2015993617092463,  0.2015993617092462, -0.2015993617092465,  0.2015993617092467, -0.2015993617092465,  0.2015993617092470, -0.2015993617092468,  0.20159936170924628]
(r, iter, conv) = solve_coefficients!(deepcopy(x), fes)

(r, iter, conv) = solve_residuals!(deepcopy(x), fes)
@test r ≈ r_ols

# PooledArrays
(r, iter, conv) = solve_residuals!(deepcopy(x), [FixedEffect(PooledArray(p1)), FixedEffect(PooledArray(p2))])
@test r ≈ r_ols

# CategorialArrays
(r, iter, conv) = solve_residuals!(deepcopy(x), [FixedEffect(categorical(p1)), FixedEffect(categorical(p2))])
@test r ≈ r_ols

(c, iter, conv) = solve_residuals!([x x], fes)

# test update_weights
weights = ones(10)
feM = FixedEffects.AbstractFixedEffectSolver{Float64}(fes, Weights(weights), Val{:cpu})
weights = Weights([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
FixedEffects.update_weights!(feM, weights)
solve_residuals!(deepcopy(x), feM)[1] ≈ solve_residuals!(deepcopy(x), fes, weights)[1]

method_s = [:cpu]
if FixedEffects.has_CUDA()
	push!(method_s, :gpu)
end
for method in method_s
	println("$method Float32")
	local (r, iter, conv) = solve_residuals!(deepcopy(x),fes, method=method, double_precision = false)
	@test Float32.(r) ≈ Float32.(r_ols)
end

fe = FixedEffect([1, 2])
@test_throws "FixedEffects must have the same length as y" ỹ = solve_residuals!(ones(100), [fe])


text:
################################################################################
#
#  Operation agnostic functionality
#
################################################################################

export closure, small_generating_set, find_identity

################################################################################
#
#  Find identity with respect to operation
#
################################################################################

# It is assumed that the elements have finite order with respect to op.
function find_identity(S, op, eq = ==)
  @assert length(S) > 0
  g = S[1]
  h = g
  while true
    hh = op(h, g)
    if eq(hh, g)
      return h
    end
    h = hh
  end
end

################################################################################
#
#  Small generating set
#
################################################################################

function small_generating_set(G, op)
  i = find_identity(G, op)
  return small_generating_set(G, op, i)
end

function _non_trivial_randelem(G, id)
  x = rand(G)::typeof(id)
  while x == id
    x = rand(G)::typeof(id)
  end
  return x
end

function small_generating_set(G::Vector, op, id)
  orderG = length(G)

  if length(G) == 1
    return G
  end

  firsttry = 10
  secondtry = 20
  thirdtry = 30

  # First try one element
  for i in 1:firsttry
    gen = _non_trivial_randelem(G, id)
    if length(closure([gen], op, id)) == orderG
      return [gen]
    end
  end

  for i in 1:secondtry
    gens = typeof(id)[_non_trivial_randelem(G, id), _non_trivial_randelem(G, id)]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end

  for i in 1:thirdtry
    gens = typeof(id)[_non_trivial_randelem(G, id), _non_trivial_randelem(G, id), _non_trivial_randelem(G, id)]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end

  # Now use that unconditionally log_2(|G|) elements generate G

  b = ceil(Int, log(2, orderG))
  @assert orderG <= 2^b

  j = 0
  while true
    if j > 2^20
      error("Something wrong with generator search")
    end
    j = j + 1
    gens = [_non_trivial_randelem(G, id) for i in 1:b]
    if length(closure(gens, op, id)) == orderG
      return unique(gens)
    end
  end
end

################################################################################
#
#  Computing closure under group operation
#
################################################################################

# It is assumed that S is nonempty and that the group generated by S under op
# is finite.
function closure(S, op; eq = ==)
  i = find_identity(S, op, eq)
  return closure(S, op, i; eq = eq)
end

function closure(S, op, id; eq = ==)
  if length(S) == 0
    return [id]
  elseif length(S) == 1
    return _closing_under_one_generator(S[1], op, id, eq = eq)
  else
    return _closing_under_generators_dimino(S, op, id, eq = eq)
  end
end

function _closing_under_generators_naive(S, op, id; eq = ==)
  list = push!(copy(S), id)
  stable = false
  while !stable
    stable = true
    for g in list
      for s in S
        m = op(g, s)
        if !(any(x -> eq(x, m), list))
          push!(list, m)
          stable = false
        end
      end
    end
  end
  return list
end

function _closing_under_one_generator(x, op, id; eq = ==)
  elements = [x]
  y = x
  while !(eq(y, id))
    y = op(y, x)
    push!(elements, y)
  end
  return elements
end

function _closing_under_generators_dimino(S, op, id; eq = ==)

  t = length(S)
  order = 1
  elements = [id]
  g = S[1]

  while !(eq(g, id))
    order = order +1
    push!(elements, g)
    g = op(g, S[1])
  end

  for i in 2:t
    if !(any(x -> eq(x, S[i]), elements))
      previous_order = order
      order = order + 1
      push!(elements, S[i])
      for j in 2:previous_order
        order = order + 1
        push!(elements, op(elements[j], S[i]))
      end

      rep_pos = previous_order + 1
      while rep_pos <= order
        for k in 1:i
          s = S[k]
          elt = op(elements[rep_pos], s)
          if !(any(x -> eq(x, elt), elements))
            order = order + 1
            push!(elements, elt)
            for j in 2:previous_order
              order = order + 1
              push!(elements, op(elements[j], elt))
            end
          end
        end
        rep_pos = rep_pos + previous_order
      end
    end
  end
  return elements
end

################################################################################
#
#  Multiplication table
#
################################################################################

# Construct multiplication table of G under op
function _multiplication_table(G, op)
  l = length(G)
  z = Matrix{Int}(undef, l, l)
  for i in 1:l
    for j in 1:l
      p = op(G[i], G[j])
      for k in 1:l
        if p == G[k]
          z[i, j] = k
          break
        end
      end
    end
  end
  return z
end


text:
# %%
using Pkg
Pkg.activate("./GraphFusedElasticNet.jl/")
Pkg.instantiate()
using GraphFusedElasticNet
using LightGraphs
using DelimitedFiles
using JSON
using DataStructures

# %%
fname ="./processed_data/spatiotemporal_graph.csv"
edgeinfo = readdlm(fname, ',', Int, skipstart=1)
edgeinfo[1:10,:]

# %% use our utils from GFEN library for trail extraction
graph = graph_from_edgelist(edgeinfo[:,1:2], from_zero=true)
trails = find_trails(graph, ntrails=64)

# %% identify temporal edges
# create temporal edge dict
tmpdict = DefaultDict{Tuple{Int, Int}, Bool}(false)
cols = [edgeinfo[:,i] for i in 1:3]
for (k, (i, j, t)) in enumerate(zip(cols...))
    x, y = min(i, j), max(i, j)
    tmpdict[x, y] = t
end
# now tag along ptr
ptr = trails.ptr
istemporal = zeros(Int, length(ptr))
brks_ = OrderedSet(trails.brks)
for (k, (i, j)) in enumerate(zip(ptr[1:end-1], ptr[2:end]))
    if k + 1 ∉ brks_
        x, y = min(i, j), max(i, j)
        istemporal[k] = tmpdict[x, y]
    end
end

# %% save fields as json
data = Dict("ptr" => trails.ptr,
            "brks" => trails.brks,
            "wts" => trails.wts,
            "istemp" => istemporal,
            "num_nodes" => trails.num_nodes)
open("./processed_data/trails.json","w") do f
    JSON.print(f, data)
end

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
# EXPERIMENTAL - Threaded version

import Base.Threads.@spawn

# Threading not possible in this case
corkendallthreads_v4(x::Union{RealVector,RealOrMissingVector}, y::Union{RealVector,RealOrMissingVector}) = corkendall(x, y)

function corkendallthreads_v4(X::Union{RealMatrix,RealOrMissingMatrix}, Y::Union{RealMatrix,RealOrMissingMatrix})
    nr = size(X, 2)
    nc = size(Y, 2)

    #speedup provided by threading is greater when the second argument has more columns, so switch.
    #TODO test that hypothesis more thoroughly!
    if nr > nc
        return(convert(RealMatrix, transpose(corkendallthreads_v4(Y, X))))
    end

    C = zeros(Float64, nr, nc)
    chunks = partitioncols(nc, false)

    tasks = Array{Task,1}(undef, length(chunks))

    for j = 1:length(chunks)
        tasks[j] = @spawn corkendall(X, Y[:,chunks[j]])
    end

    for (c,t) in zip(chunks,tasks)
        C[:,c] = fetch(t)
    end

    return C
end

function corkendallthreads_v4(X::Union{RealMatrix,RealOrMissingMatrix})
    nr = nc = size(X, 2)
    C = zeros(Float64, nr, nc)

    chunks = partitioncols(nc, true,max(Threads.nthreads(),nc/100))

    tasks = Array{Task,1}(undef, length(chunks))
    for j = 1:length(chunks)
        tasks[j] = @spawn ck_belowdiagonal(X, chunks[j])
    end

    for (c,t) in zip(chunks,tasks)
        C[:,c] = fetch(t)
    end

    for j = 1:nr
        C[j,j] = 1.0
        for i = (1 + j):nr
            C[j,i] = C[i,j]
        end
    end

    return C
end

function corkendallthreads_v4(x::Union{RealVector,RealOrMissingVector}, Y::Union{RealMatrix,RealOrMissingMatrix})
    nr = 1
    nc = size(Y, 2)

    C = zeros(Float64, nr, nc)

    chunks = partitioncols(nc, false,max(8,nc/20))

    tasks = Array{Task,1}(undef, length(chunks))
    for j = 1:length(chunks)
        tasks[j] = @spawn corkendall(x, view(Y,:,chunks[j]))
    end

    for (c,t) in zip(chunks,tasks)
        C[:,c] = fetch(t)
    end

    return C
end

function corkendallthreads_v4(X::Union{RealMatrix,RealOrMissingMatrix}, y::Union{RealVector,RealOrMissingVector})
    l = size(X, 2)
    #it is idiosyncratic that this method returns a vector, not a matrix.
    C = zeros(Float64, l)

    chunks = partitioncols(l, false)

    tasks = Array{Task,1}(undef, length(chunks))
    for j = 1:length(chunks)
        tasks[j] = @spawn corkendall(X[:,chunks[j]], y)
    end

    for (c,t) in zip(chunks,tasks)
        C[c] = fetch(t)
    end

    return C
end

"""
    partitioncols(nc::Int64, triangular::Bool)
Auxiliary function for task load balancing. Returns a vector of `UnitRange`s, which partition the columns of the
correlation matrix to be calculated, one "chunk" per task.
# Arguments
- `nc::Int64`: the number of columns in the correlation matrix.
- `triangular::Bool`: should be `true` when calculating below-the-diagonal elements of a correlation matrix. In this
case the partitions have (approximately) equal numbers of elements below the diagonal, so early partitions are narrower
than later partitions.
- `numtasks`: the number of tasks (or chunks) to partition to. The length of the return from the function
is `min(nc,numtasks)`.

# Examples
```
julia> length.(KendallTau.partitioncols(10,false,4))
4-element Vector{Int64}:
 3
 3
 2
 2

julia> length.(KendallTau.partitioncols(100,true,4))
4-element Vector{Int64}:
 14
 16
 21
 48
```
"""
function partitioncols(nc::Int64, triangular::Bool, numtasks=Threads.nthreads())

    chunks = Vector{UnitRange{Int64}}(undef, 0)
    if triangular
        lastcol = nc - 1
        numcorrels = nc * (nc - 1) ÷ 2
    else
        lastcol = nc
        numcorrels = nc * nc
    end

    correlsdone = tasksdone = 0
    starttask = true
    chunkfirstcol = 0
    target = 0.0
    for i = 1:lastcol
        if starttask
            chunkfirstcol = i
            target = correlsdone + (numcorrels - correlsdone) / (numtasks - tasksdone)
        end
        correlsdone += triangular ? (nc - i) : nc
        if correlsdone >= target || i == lastcol
            push!(chunks, chunkfirstcol:i)
            starttask = true
            tasksdone += 1
        else
            starttask = false
        end
    end

    return(chunks)
end

"""
    ck_belowdiagonal(X::Union{RealMatrix,RealOrMissingMatrix}, colnos::UnitRange{Int64})
For use from multi-threading code to avoid double calculation of elements. Returns corkendall(X)[:,colnos] but with NaNs
on and above the diagonal of corkendall(X).
"""
function ck_belowdiagonal(X::Union{RealMatrix,RealOrMissingMatrix}, colnos::UnitRange{Int64})
    nr = size(X, 2)
    nc = length(colnos)
    C = Matrix{Float64}(undef, nr, nc)
    for j = 1:nr
        permx = sortperm(X[:,j])
        sortedx = X[:,j][permx]
        for i = 1:nc
            if j > i + colnos[1] - 1
                C[j,i] = ck_sorted!(sortedx, X[:,colnos[i]], permx)
            else
                C[j,i] = NaN
            end
        end
    end
    return C
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Cairomm"
version = v"1.16.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.cairographics.org/releases/cairomm-$(version).tar.xz",
                  "6f6060d8e98dd4b8acfee2295fddbdd38cf487c07c26aad8d1a83bb9bff4a2c6")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cairomm*/
mkdir output && cd output/
meson --cross-file=${MESON_TARGET_TOOLCHAIN} ..
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct(["libcairomm-$(version.major).$(version.minor)", "libcairomm-$(version.major)"], :libcariomm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libsigcpp_jll", uuid="8855df87-3ca1-5a3d-a68b-4f0ddc198ba8"))
    Dependency(PackageSpec(name="Cairo_jll", uuid="83423d85-b0ee-5818-9007-b63ccbeb887a"); compat=string(version))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
export AbstractSpaceSymmetryOperation

export Permutation
export generate_group
export isidentity

"""
    Permutation(perms; max_order=2048)

Create a permutation of integers from 1 to n.
`perms` should be a permutation of `1:n`.

# Arguments
- `perms`: an integer vector containing a permutation of integers from 1 to n
- `max_order`: maximum order

# Note
The convention for the permutation is that map[i] gets mapped to i.
In other words, map tells you where each element is from.
"""
struct Permutation
    map ::Vector{Int}
    order ::Int
    function Permutation(perms::AbstractVector{<:Integer}; max_order=2048)
        n = length(perms)
        map = Vector{Int}(perms)
        let # check for duplicates
            duplicates = Set{Int}()
            for j in perms
                if !(1 <= j <= n)
                    throw(ArgumentError("argument not a proper permutation (target != universe)"))
                elseif j in duplicates
                    throw(ArgumentError("argument not a proper permutation (contains duplicates)"))
                end
                push!(duplicates, j)
            end
        end

        order = 1
        let # compute cycle length
            current = Int[map[x] for x in 1:n]
            while !issorted(current) && order <= max_order
                current = Int[map[x] for x in current]
                order += 1
            end
            if order > max_order
                throw(OverflowError("cycle length exceeds maximum value (max = $max_order)"))
            end
        end
        return new(map, order)
    end
end

"""
    *(p1 ::Permutation, p2 ::Permutation)

Multiply the two permutation.
NOT THIS: (Return `[p2.map[x] for x in p1.map]`.)
BUT THIS: (Return `[p1.map[x] for x in p2.map]`.)

# Examples
```jldoctest
julia> using TightBindingLattice

julia> Permutation([2,1,3]) * Permutation([1,3,2])
Permutation([2, 3, 1], 3)

julia> Permutation([1,3,2]) * Permutation([2,1,3])
Permutation([3, 1, 2], 3)
```
"""
function Base.:(*)(p1 ::Permutation, p2 ::Permutation)
    if length(p1.map) != length(p2.map)
        throw(ArgumentError("permutations of different universes"))
    end
    #  return Permutation(Int[p2.map[x] for x in p1.map]) original version
    return Permutation(Int[p1.map[x] for x in p2.map])
end

#=
function Base.:(*)(lhs ::AbstractSet{Permutation}, rhs::Permutation)
    return generate_group(lhs..., rhs)
end

function Base.:(*)(lhs ::Permutation, rhs::AbstractSet{Permutation})
    return generate_group(lhs, rhs...)
end

function Base.:(*)(lhs ::AbstractSet{Permutation}, rhs::AbstractSet{Permutation})
    return generate_group(lhs..., rhs...)
end
=#

"""
    ^(perm ::Permutation, pow ::Integer)

Exponentiate the permutation.

# Examples
```jldoctest
julia> using TightBindingLattice

julia> Permutation([2,3,4,1])^2
Permutation([3, 4, 1, 2], 2)
```
"""
function Base.:(^)(perm::Permutation, pow::Integer)
    p = mod(pow, perm.order)
    out = collect(1:length(perm.map))
    for i in 1:p
        out = collect(perm.map[x] for x in out)
    end
    return Permutation(out)
end

function Base.inv(perm::Permutation)
    out = zeros(Int, length(perm.map))
    for (i, x) in enumerate(perm.map)
        out[x] = i
    end
    return Permutation(out)
end

Base.:(==)(p1 ::Permutation, p2::Permutation) = p1.map == p2.map

(p::Permutation)(i::Integer) = p.map[i]

function Base.isless(p1 ::Permutation, p2::Permutation)
    return Base.isless(p1.order, p2.order) || ((p1.order == p2.order) && Base.isless(p1.map, p2.map))
end

function Base.hash(p::Permutation, h::UInt=UInt(0x0))
    return Base.hash(p.map, hash(Permutation, h))
end

"""
    generate_group(generators...)

Return a FiniteGroup generated by the `generators`.

# Arguments
* `generators::Permutation...`: generating permutations
"""
function generate_group(generators::Permutation...)
    change = true
    group = Set{Permutation}([generators...])
    while change
        change = false
        for g1 in generators, g2 in group
            g3 = g1 * g2
            if !(g3 in group)
                change = true
                push!(group, g3)
            end
        end
    end
    return group
end

"""
    isidentity(perm::Permutation)

Test whether the permutation is an identity.
"""
isidentity(perm::Permutation) = perm.map == 1:length(perm.map)


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	thermal_no_commit(EP::Model, inputs::Dict, Reserves::Int)

This function defines the operating constraints for thermal power plants NOT subject to unit commitment constraints on power plant start-ups and shut-down decisions ($y \in H \setminus UC$).

**Ramping limits**

Thermal resources not subject to unit commitment ($y \in H \setminus UC$) adhere instead to the following ramping limits on hourly changes in power output:

```math
\begin{aligned}
	\Theta_{y,z,t-1} - \Theta_{y,z,t} \leq \kappa_{y,z}^{down} \Delta^{\text{total}}_{y,z} \hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```

```math
\begin{aligned}
	\Theta_{y,z,t} - \Theta_{y,z,t-1} \leq \kappa_{y,z}^{up} \Delta^{\text{total}}_{y,z} \hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```
(See Constraints 1-2 in the code)

This set of time-coupling constraints wrap around to ensure the power output in the first time step of each year (or each representative period), $t \in \mathcal{T}^{start}$, is within the eligible ramp of the power output in the final time step of the year (or each representative period), $t+\tau^{period}-1$.

**Minimum and maximum power output**

When not modeling regulation and reserves, thermal units not subject to unit commitment decisions are bound by the following limits on maximum and minimum power output:

```math
\begin{aligned}
	\Theta_{y,z,t} \geq \rho^{min}_{y,z} \times \Delta^{total}_{y,z}
	\hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```

```math
\begin{aligned}
	\Theta_{y,z,t} \leq \rho^{max}_{y,z,t} \times \Delta^{total}_{y,z}
	\hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```
(See Constraints 3-4 in the code)
"""
function thermal_no_commit(EP::Model, inputs::Dict, Reserves::Int)

	println("Thermal (No Unit Commitment) Resources Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	hours_per_subperiod = inputs["hours_per_subperiod"] #total number of hours per subperiod

	START_SUBPERIODS = inputs["START_SUBPERIODS"]
	INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]

	THERM_NO_COMMIT = inputs["THERM_NO_COMMIT"]

	### Expressions ###

	## Power Balance Expressions ##
	@expression(EP, ePowerBalanceThermNoCommit[t=1:T, z=1:Z],
		sum(EP[:vP][y,t] for y in intersect(THERM_NO_COMMIT, dfGen[dfGen[!,:Zone].==z,:][!,:R_ID])))

	EP[:ePowerBalance] += ePowerBalanceThermNoCommit

	### Constraints ###

	### Maximum ramp up and down between consecutive hours (Constraints #1-2)
	@constraints(EP, begin

		## Maximum ramp up between consecutive hours
		# Start Hours: Links last time step with first time step, ensuring position in hour 1 is within eligible ramp of final hour position
		# NOTE: We should make wrap-around a configurable option
		[y in THERM_NO_COMMIT, t in START_SUBPERIODS], EP[:vP][y,t]-EP[:vP][y,(t+hours_per_subperiod-1)] <= dfGen[!,:Ramp_Up_Percentage][y]*EP[:eTotalCap][y]

		# Interior Hours
		[y in THERM_NO_COMMIT, t in INTERIOR_SUBPERIODS], EP[:vP][y,t]-EP[:vP][y,t-1] <= dfGen[!,:Ramp_Up_Percentage][y]*EP[:eTotalCap][y]

		## Maximum ramp down between consecutive hours
		# Start Hours: Links last time step with first time step, ensuring position in hour 1 is within eligible ramp of final hour position
		[y in THERM_NO_COMMIT, t in START_SUBPERIODS], EP[:vP][y,(t+hours_per_subperiod-1)] - EP[:vP][y,t] <= dfGen[!,:Ramp_Dn_Percentage][y]*EP[:eTotalCap][y]

		# Interior Hours
		[y in THERM_NO_COMMIT, t in INTERIOR_SUBPERIODS], EP[:vP][y,t-1] - EP[:vP][y,t] <= dfGen[!,:Ramp_Dn_Percentage][y]*EP[:eTotalCap][y]
	end)

	### Minimum and maximum power output constraints (Constraints #3-4)
	if Reserves == 1
		# If modeling with regulation and reserves, constraints are established by thermal_no_commit_reserves() function below
		EP = thermal_no_commit_reserves(EP, inputs)
	else
		@constraints(EP, begin
			# Minimum stable power generated per technology "y" at hour "t" Min_Power
			[y in THERM_NO_COMMIT, t=1:T], EP[:vP][y,t] >= dfGen[!,:Min_Power][y]*EP[:eTotalCap][y]

			# Maximum power generated per technology "y" at hour "t"
			[y in THERM_NO_COMMIT, t=1:T], EP[:vP][y,t] <= inputs["pP_Max"][y,t]*EP[:eTotalCap][y]
		end)

	end
	# END Constraints for thermal resources not subject to unit commitment

	return EP
end

@doc raw"""
	thermal_no_commit_reservesEP::Model, inputs::Dict)

This function is called by the ```thermal_no_commit()``` function when regulation and reserves constraints are active and defines reserve related constraints for thermal power plants not subject to unit commitment constraints on power plant start-ups and shut-down decisions.

**Maximum contributions to frequency regulation and reserves**

Thermal units not subject to unit commitment adhere instead to the following constraints on maximum reserve and regulation contributions:

```math
\begin{aligned}
	f_{y,z,t} \leq \upsilon^{reg}_{y,z} \times \rho^{max}_{y,z,t} \Delta^{\text{total}}_{y,z} \hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```

```math
\begin{aligned}
	r_{y,z,t} \leq \upsilon^{rsv}_{y,z} \times \rho^{max}_{y,z,t} \Delta^{\text{total}}_{y,z} \hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```

where $f_{y,z,t}$ is the frequency regulation contribution limited by the maximum regulation contribution $\upsilon^{reg}_{y,z}$, and $r_{y,z,t}$ is the reserves contribution limited by the maximum reserves contribution $\upsilon^{rsv}_{y,z}$. Limits on reserve contributions reflect the maximum ramp rate for the thermal resource in whatever time interval defines the requisite response time for the regulation or reserve products (e.g., 5 mins or 15 mins or 30 mins). These response times differ by system operator and reserve product, and so the user should define these parameters in a self-consistent way for whatever system context they are modeling.

**Minimum and maximum power output**

When modeling regulation and spinning reserves, thermal units not subject to unit commitment are bound by the following limits on maximum and minimum power output:

```math
\begin{aligned}
	\Theta_{y,z,t} - f_{y,z,t} \geq \rho^{min}_{y,z} \times \Delta^{\text{total}}_{y,z}
	\hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```

```math
\begin{aligned}
	\Theta_{y,z,t} + f_{y,z,t} + r_{y,z,t} \leq \rho^{max}_{y,z,t} \times \Delta^{\text{total}}_{y,z}
	\hspace{1cm} \forall y \in \mathcal{H \setminus UC}, \forall z \in \mathcal{Z}, \forall t \in \mathcal{T}
\end{aligned}
```

Note there are multiple versions of these constraints in the code in order to avoid creation of unecessary constraints and decision variables for thermal units unable to provide regulation and/or reserves contributions due to input parameters (e.g. ```Reg_Max=0``` and/or ```RSV_Max=0```).
"""
function thermal_no_commit_reserves(EP::Model, inputs::Dict)

	println("Thermal (No Unit Commitment) Reserves Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)

	THERM_NO_COMMIT = setdiff(inputs["THERM_ALL"], inputs["COMMIT"])

	# Constraints on contribution to regulation and reserves
	# Thermal units without commitment constraints assumed to be quick-start units, so can contribute a specified fraction of their
	# total installed capacity to reserves

	THERM_NO_COMMIT_REG_RSV = intersect(THERM_NO_COMMIT, inputs["REG"], inputs["RSV"]) # Set of thermal resources with both regulation and spinning reserves

	THERM_NO_COMMIT_REG = intersect(THERM_NO_COMMIT, inputs["REG"]) # Set of thermal resources with regulation reserves
	THERM_NO_COMMIT_RSV = intersect(THERM_NO_COMMIT, inputs["RSV"]) # Set of thermal resources with spinning reserves

	THERM_NO_COMMIT_NO_RES = setdiff(THERM_NO_COMMIT, THERM_NO_COMMIT_REG, THERM_NO_COMMIT_RSV) # Set of thermal resources with no reserves

	THERM_NO_COMMIT_REG_ONLY = setdiff(THERM_NO_COMMIT_REG, THERM_NO_COMMIT_RSV) # Set of thermal resources only with regulation reserves
	THERM_NO_COMMIT_RSV_ONLY = setdiff(THERM_NO_COMMIT_RSV, THERM_NO_COMMIT_REG) # Set of thermal resources only with spinning reserves

	if !isempty(THERM_NO_COMMIT_REG_RSV)
		@constraints(EP, begin

			[y in THERM_NO_COMMIT_REG_RSV, t=1:T], EP[:vREG][y,t] <= inputs["pP_Max"][y,t]*dfGen[!,:Reg_Max][y]*EP[:eTotalCap][y]
			[y in THERM_NO_COMMIT_REG_RSV, t=1:T], EP[:vRSV][y,t] <= inputs["pP_Max"][y,t]*dfGen[!,:Rsv_Max][y]*EP[:eTotalCap][y]

			# Minimum stable power generated per technology "y" at hour "t" and contribution to regulation down must be > min power
			[y in THERM_NO_COMMIT_REG_RSV, t=1:T], EP[:vP][y,t]-EP[:vREG][y,t] >= dfGen[!,:Min_Power][y]*EP[:eTotalCap][y]

			# Maximum power generated per technology "y" at hour "t" and contribution to regulation and reserves must be < max power
			[y in THERM_NO_COMMIT_REG_RSV, t=1:T], EP[:vP][y,t]+EP[:vREG][y,t]+EP[:vRSV][y,t] <= inputs["pP_Max"][y,t]*EP[:eTotalCap][y]
		end)
	end

	if !isempty(THERM_NO_COMMIT_REG)
		@constraints(EP, begin

			[y in THERM_NO_COMMIT_REG, t=1:T], EP[:vREG][y,t] <= inputs["pP_Max"][y,t]*dfGen[!,:Reg_Max][y]*EP[:eTotalCap][y]

			# Minimum stable power generated per technology "y" at hour "t" and contribution to regulation down must be > min power
			[y in THERM_NO_COMMIT_REG, t=1:T], EP[:vP][y,t]-EP[:vREG][y,t] >= dfGen[!,:Min_Power][y]*EP[:eTotalCap][y]

			# Maximum power generated per technology "y" at hour "t" and contribution to regulation and reserves must be < max power
			[y in THERM_NO_COMMIT_REG, t=1:T], EP[:vP][y,t]+EP[:vREG][y,t] <= inputs["pP_Max"][y,t]*EP[:eTotalCap][y]
		end)
	end

	if !isempty(THERM_NO_COMMIT_RSV)
		@constraints(EP, begin

			[y in THERM_NO_COMMIT_RSV, t=1:T], EP[:vRSV][y,t] <= inputs["pP_Max"][y,t]*dfGen[!,:Rsv_Max][y]*EP[:eTotalCap][y]

			# Minimum stable power generated per technology "y" at hour "t" and contribution to regulation down must be > min power
			[y in THERM_NO_COMMIT_RSV, t=1:T], EP[:vP][y,t] >= dfGen[!,:Min_Power][y]*EP[:eTotalCap][y]

			# Maximum power generated per technology "y" at hour "t" and contribution to regulation and reserves must be < max power
			[y in THERM_NO_COMMIT_RSV, t=1:T], EP[:vP][y,t]+EP[:vRSV][y,t] <= inputs["pP_Max"][y,t]*EP[:eTotalCap][y]
		end)
	end

	if !isempty(THERM_NO_COMMIT_NO_RES)
		@constraints(EP, begin
			# Minimum stable power generated per technology "y" at hour "t" Min_Power
			[y in THERM_NO_COMMIT_NO_RES, t=1:T], EP[:vP][y,t] >= dfGen[!,:Min_Power][y]*EP[:eTotalCap][y]

			# Maximum power generated per technology "y" at hour "t"
			[y in THERM_NO_COMMIT_NO_RES, t=1:T], EP[:vP][y,t] <= inputs["pP_Max"][y,t]*EP[:eTotalCap][y]
		end)
	end

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM