text:
export pseudo_matrix, pseudo_hnf, PseudoMatrix, pseudo_hnf_with_transform, coefficient_ideals, matrix
import Base.vcat, Base.hcat

add_verbose_scope(:PseudoHnf)
add_assert_scope(:PseudoHnf)

function _det_bound(M::Generic.Mat{NfOrdElem})
  n = nrows(M)
  O = base_ring(M)
  d = degree(O)
  c1, c2 = Hecke.norm_change_const(O)

  return fmpz(BigInt(ceil(sqrt(c2)*c1^(n/2)*BigInt(n)^n*BigInt(d)^n*BigInt(_max_max(M))^n)))
end

function _max_max(M::Generic.Mat{NfOrdElem})
  d = FlintZZ(1)
  for i in 1:nrows(M)
    for j in 1:ncols(M)
      if !iszero(M[i, j])
        v = coordinates(M[i, j])
        for k in degree(base_ring(M))
          d = max(d, abs(v[k]))
        end
      end
    end
  end
  return d
end

#function det(M::Generic.Mat{NfOrdElem})
#  O = base_ring(M)::NfOrd
#  B = _det_bound(M)
#  p = next_prime(2^60) # magic numbers
#  P = fmpz(1)
#  i = 1
#  res = zero(O)
#  t = 0.0
#  while P < 2*B
#    # reject some bad primes
#    if true
#      #println("$p");
#      Omodp, pi_p = quo(O, ideal(O, p))
#      Mmodp = MatrixSpace(Omodp, nrows(M), ncols(M))(M)
#      t += @elapsed detmodp = pi_p\Hecke.det(Mmodp)
#      if i == 1
#        res = detmodp
#      elseif i >= 2
#        g, e, f = gcdx(P, fmpz(p))
#        @assert isone(g)
#        res = f*p*res + e*P*detmodp
#        res = mod_sym(res, P*p)
#        #@assert _basis_coord_nonneg(res)
#      end
#      P = P*p
#      p = next_prime(p)
#      i += 1
#    end
#  end
#  #println("Modular determinant time: $t");
#  return res
#end

function _get_coeff_raw(x::nmod_poly, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, libflint), UInt, (Ref{nmod_poly}, Int), x, i)
  return u
end

function _get_coeff_raw(x::gfp_poly, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, libflint), UInt, (Ref{gfp_poly}, Int), x, i)
  return u
end

function det(M::Generic.Mat{NfOrdElem})
  O = base_ring(M)::NfOrd
  K = nf(O)
  B = _det_bound(M)
  p = next_prime(p_start) #global
  P = fmpz(1)
  i = 1
  res = O()
  t_det = 0.0
  t_reducing = 0.0
  t_splitting = 0.0
  used_primes = Vector{UInt}()
  tmp_polys = Vector{nf_elem}()

  while P < 2*B
    # reject some bad primes
    if isindex_divisor(O, p)
      continue
    end

    ttt = @elapsed me = modular_init(K, p)

    push!(used_primes, UInt(p))

    t_splitting += ttt

    detmodp = nmod_poly(UInt(p))

    t_reducing += @elapsed Mp = modular_proj(M, me)

    t_det += @elapsed detmodQ = [det(x) for x = Mp]

      # now the CRT step
    detmodp = modular_lift(detmodQ, me)

    push!(tmp_polys, detmodp)

    P = P*p
    p = next_prime(p)
    i += 1
  end

  # now the CRT on each coefficient

  fc = crt_env(fmpz[x for x = used_primes])
  fv = Array{fmpz}(undef, length(used_primes))
  for i=1:length(used_primes)
    fv[i] = fmpz(0)
  end
  zz = fmpz()

  @assert length(used_primes) == length(tmp_polys)

  tmp_fmpz_poly = PolynomialRing(FlintZZ)[1]()

  for i in 0:degree(O)
    for j=1:length(used_primes)
      Nemo.num_coeff!(fv[j], tmp_polys[j], i)
    end
    crt!(zz, fv, fc)
    setcoeff!(tmp_fmpz_poly, i, zz)
  end

  tut = fmpq_poly(tmp_fmpz_poly)
  tut.parent = parent(nf(O).pol)
  res = mod_sym(O(nf(O)(tut)), P)

  #println("Modular determinant time: $t_det");
  #println("Time for reducing: $t_reducing");
  #println("Time for splitting: $t_splitting");
  #println("Used $(length(used_primes)) primes")

  return res
end

# s, t are auxillary variables, r1, r2 are the residues, m1, m2 are the moduli
# aliasing is not allowed (?)
function crt!(z::nmod_poly, r1::nmod_poly, r2::Union{nmod_poly, fq_nmod}, m1::nmod_poly, m2::nmod_poly, s::nmod_poly, t::nmod_poly)
  ccall((:nmod_poly_xgcd, libflint), Nothing, (Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}, Ref{nmod_poly}), z, s, t, m1, m2)
  @assert Bool(ccall((:nmod_poly_is_one, libflint), Cint, (Ref{nmod_poly}, ), z))
  # z = s*m1*r2 + t*m2*r1
  mul!(z, s, m1)
  ccall((:nmod_poly_mul, libflint), Nothing, (Ref{nmod_poly}, Ref{nmod_poly}, Ref{fq_nmod}), z, z, r2)
  mul!(t, t, m2)
  mul!(t, t, r1)
  add!(z, z, t)
  mul!(s, m1, m2)
  rem!(z, z, s)
end

function set!(z::nmod_poly, x::nmod_poly)
  ccall((:nmod_poly_set, libflint), Nothing, (Ref{nmod_poly}, Ref{nmod_poly}), z, x)
end

function set!(z::gfp_poly, x::gfp_poly)
  ccall((:nmod_poly_set, libflint), Nothing, (Ref{gfp_poly}, Ref{gfp_poly}), z, x)
end

function __helper!(z, mF, entries)
  s = size(entries)
  for i in 1:s[2]
    for j in 1:s[1]
      z[j, i] = mF(entries[j, i])
    end
  end
  return z
end

function mod_sym(x::NfOrdElem, m::fmpz)
  z = coordinates(x)
  for i in 1:length(z)
    z[i] = mod(z[i], m)
    if 2*z[i] > m
      z[i] = z[i] - m
    end
  end
  return parent(x)(z)
end
mod_sym(x::NfOrdElem, m::Integer) = mod_sym(x, fmpz(m))

function _basis_coord_nonneg(x::NfOrdElem)
  b = coordinates(x)
  for i in 1:length(b)
    if b[i] < 0
      return false
    end
  end
  return true
end

function rand(M::Generic.MatSpace{NfOrdElem}, B::Union{Int, fmpz})
  z = M()
  for i in 1:nrows(z)
    for j in 1:ncols(z)
      z[i, j] = Hecke.rand(Hecke.base_ring(M), B)
    end
  end
  return z
end

==(P::PMat, Q::PMat) = P.matrix == Q.matrix && P.coeffs == Q.coeffs

function Base.deepcopy_internal(P::PMat{T, S}, dict::IdDict) where {T, S}
  z = PMat{T, S}()
  for x in fieldnames(typeof(P))
    if x != :parent && isdefined(P, x)
      setfield!(z, x, Base.deepcopy_internal(getfield(P, x), dict))
    end
  end
  if isdefined(P, :parent)
    z.parent = P.parent
  end
  return z
end

function show(io::IO, P::PMat)
  compact = get(io, :compact, false)
  if compact
    for i in 1:nrows(P.matrix)
      i == 1 || print(io, "\n")
      print(io, "(")
      print(IOContext(io, :compact => true), P.coeffs[i])
      print(io, ") * ")
      print(io, sub(P.matrix, i:i, 1:ncols(P.matrix)))
    end
  else
    print(io, "Pseudo-matrix over $(parent(P.matrix[1, 1]))")
    for i in 1:nrows(P.matrix)
      print(io, "\n")
      show(IOContext(io, :compact => true), P.coeffs[i])
      print(io, " with row $(sub(P.matrix, i:i, 1:ncols(P.matrix)))")
    end
  end
end

@doc Markdown.doc"""
    coefficient_ideals(M::PMat)

Returns the vector of coefficient ideals.
"""
coefficient_ideals(M::PMat) = M.coeffs

@doc Markdown.doc"""
    matrix(M::PMat)

Returns the matrix part of the `PMat`.
"""
matrix(M::PMat) = M.matrix

@doc Markdown.doc"""
    base_ring(M::PMat)

The `PMat` $M$ defines an $R$-module for soem maximal order $R$.
This function returns the $R$ that was used to defined $M$.
"""
base_ring(M::PMat) = order(M.coeffs[1])

function PseudoMatrix(m::AbstractAlgebra.MatElem{T}, c::Vector{S}) where {T, S}
  # sanity checks
  @assert nrows(m) == length(c)
  return PMat{T, S}(m ,c)
end

@doc Markdown.doc"""
    PseudoMatrix(m::Generic.Mat{nf_elem}, c::Vector{NfOrdIdl}) -> PMat{nf_elem, NfOrdFracIdl}

Returns the (row) pseudo matrix representing the $Z_k$-module
 $$\sum c_i m_i$$
 where $c_i$ are the ideals in $c$ and $m_i$ the rows of $M$.
"""
function PseudoMatrix(m::AbstractAlgebra.MatElem{nf_elem}, c::Vector{NfOrdIdl})
  @assert nrows(m) == length(c)
  cc = map(z -> NfOrdFracIdl(z, fmpz(1)), c)
  return PMat{nf_elem, NfOrdFracIdl}(m, cc)
end

@doc Markdown.doc"""
    PseudoMatrix(m::Generic.Mat{NfOrdElem}, c::Vector{NfOrdIdl}) -> PMat{nf_elem, NfOrdFracIdl}

Returns the (row) pseudo matrix representing the $Z_k$-module
 $$\sum c_i m_i$$
 where $c_i$ are the ideals in $c$ and $m_i$ the rows of $M$.
"""
function PseudoMatrix(m::Generic.Mat{NfOrdElem}, c::Vector{NfOrdIdl})
  @assert nrows(m) == length(c)
  mm = change_base_ring(nf(base_ring(m)), m)
  cc = map(z -> NfOrdFracIdl(z, fmpz(1)), c)
  return PMat{nf_elem, NfOrdFracIdl}(mm, cc)
end

@doc Markdown.doc"""
    PseudoMatrix(m::Generic.Mat{NfOrdElem}) -> PMat{nf_elem, NfOrdFracIdl}

Returns the free (row) pseudo matrix representing the $Z_k$-module
 $$\sum Z_k m_i$$
 where $m_i$ are the rows of $M$.
"""
function PseudoMatrix(m::Generic.Mat{nf_elem})
   K = base_ring(m)
   O = maximal_order(K)
   return PseudoMatrix(m, [ideal(O, K(1)) for i = 1:nrows(m)])
end

function PseudoMatrix(m::MatElem{S}) where S <: NumFieldElem
  L = base_ring(m)
  OL = maximal_order(L)
  K = base_field(L)
  return PseudoMatrix(m, fractional_ideal_type(OL)[ fractional_ideal(OL, ideal(OL, 1)) for i = 1:nrows(m) ])
end

function PseudoMatrix(m::MatElem{S}, c::Vector{T}) where {S <: NumFieldElem, T <: NfRelOrdIdl}
  @assert nrows(m) == length(c)
  cc = [ fractional_ideal(order(c[i]), basis_pmatrix(c[i]), true) for i = 1:length(c) ]
  return PMat{S, typeof(cc[1])}(m, cc)
end

PseudoMatrix(m::MatElem{NfOrdElem}) = PseudoMatrix(change_base_ring(nf(base_ring(m)), m))

function PseudoMatrix(c::Vector{S}) where S
   K = nf(order(c[1]))
   m = identity_matrix(K, length(c))
   return PseudoMatrix(m, c)
end

PseudoMatrix(c::Vector{NfOrdIdl}) = PseudoMatrix(map(z -> NfOrdFracIdl(z, fmpz(1)), c))

function nrows(m::PMat)
  return nrows(m.matrix)
end

function ncols(m::PMat)
  return ncols(m.matrix)
end

function det(m::PMat)
  z = m.coeffs[1]
  for i in 2:nrows(m)
    z = z*m.coeffs[i]
  end
  return det(m.matrix)*z
end

function *(P::PMat{T, S}, x::U) where { T, S, U <: Union{ Int, RingElem } }
  if nrows(P) == 0 || ncols(P) == 0
    return P
  end

  K = parent(P.matrix[1, 1])
  x = K(x)

  PP = deepcopy(P)
  for i = 1:nrows(PP)
    PP.coeffs[i] = PP.coeffs[i]*x
    PP.coeffs[i] = simplify(PP.coeffs[i])
  end
  return PP
end

*(x::U, P::PMat{T, S}) where { T, S, U <: Union{ Int, RingElem } } = P*x

# this is slow
function _coprime_integral_ideal_class(x::Union{NfOrdFracIdl, NfOrdIdl}, y::NfOrdIdl)
  O = order(y)
  x_inv = inv(x)
  check = true
  z = ideal(O, O(1))
  a = nf(O)()
  i = 0
  while check
    i += 1
    a = rand(x_inv, 10)
    if iszero(a)
      continue
    end
    b = x*a
    z = divexact(numerator(b), denominator(b))
    norm(z + y) == 1 ? (check = false) : (check = true)
  end
  return z, a
end

function _coprime_integral_ideal_class(x::FacElem{NfOrdIdl, NfOrdIdlSet}, y::NfOrdIdl)
  D = typeof(x.fac)()
  D2 = Dict{nf_elem, fmpz}()
  O = order(y)
  for (I, e) in x
    II, a = _coprime_integral_ideal_class(I, y)
    D[II] = e
    D2[a] = e
  end
  return FacElem(D), FacElem(D2)
end

# this is slow
function _coprime_norm_integral_ideal_class(x, y) #x::NfOrdFracIdl, y::NfOrdIdl)
  # x must be nonzero
  O = order(y)
  if iscoprime(norm(numerator(x, copy = false), copy = false), norm(y, copy = false))
    return numerator(x, copy = false), nf(O)(denominator(x))
  end
  x_inv = inv(x)
  check = true
  z = ideal(O, O(1))
  a = nf(O)()
  i = 0
  while check && i < 20
    a = rand(x_inv, 10)
    if iszero(a)
      continue
    end
    i += 1
    b = x*a
    simplify(b)
    @assert isone(denominator(b, copy = false))
    z = numerator(b, copy = false)
    check = !(gcd(norm(z, copy = false), norm(y, copy = false)) == 1)
  end
  if !check
    return z, a
  end
  a = nf(O)(denominator(x, copy = false))
  lp = factor(ideal(O, gcd(minimum(numerator(x, copy = false), copy = false), minimum(y, copy = false))))
  J, b = coprime_deterministic(numerator(x, copy = false), y, lp)
  res2 = b*a
  @hassert :PseudoHnf 1 res2*x == J
  @hassert :PseudoHnf 1 iscoprime(norm(J, copy = false), norm(y, copy = false))
  return J, res2
end

RandomExtensions.maketype(I::NfOrdIdl, ::Int) = NfOrdElem

function rand(rng::AbstractRNG, sp::SamplerTrivial{<:Make2{NfOrdElem,NfOrdIdl,Int}})
  I, B = sp[][1:end]
  r = rand(rng, -B:B, degree(order(I)))
  b = basis(I)
  z::Random.gentype(sp) = r[1]*b[1] # type assert to help inference in Julia 1.0
  for i in 2:degree(order(I))
    z = z + r[i]*b[i]
  end
  return z
end

rand(I::NfOrdIdl, B::Int) = rand(GLOBAL_RNG, I, B)
rand(rng::AbstractRNG, I::NfOrdIdl, B::Int) = rand(rng, make(I, B))

RandomExtensions.maketype(I::NfOrdFracIdl, ::Int) = nf_elem

function rand(rng::AbstractRNG, sp::SamplerTrivial{<:Make2{nf_elem,NfOrdFracIdl,Int}})
  I, B = sp[][1:2]
  z = rand(rng, make(numerator(I), B))
  return divexact(elem_in_nf(z), denominator(I))::nf_elem
  # type assert to help inference in Julia 1.0
end

rand(I::NfOrdFracIdl, B::Int) = rand(GLOBAL_RNG, I, B)
rand(rng::AbstractRNG, I::NfOrdFracIdl, B::Int) = rand(rng, make(I, B))

@doc Markdown.doc"""
    pseudo_hnf(P::PMat)

Transforms $P$ into pseudo-Hermite form as defined by Cohen. Essentially the
matrix part of $P$ will be upper triangular with some technical normalisation
for the off-diagonal elements. This operation preserves the module.

A optional second argument can be secified as a symbols, indicating the desiered
shape of the echelon form. Possible are
`:upperright` (the default) and `:lowerleft`
"""
function pseudo_hnf(P::PMat{nf_elem, NfOrdFracIdl}, shape::Symbol = :upperright, full_rank::Bool = false)
  if full_rank
    return pseudo_hnf_full_rank(P, shape)
  else
    # TODO: If P is not of full rank and nrows(P) > ncols(P)
    # find_pseudo_hnf_modulus (called by pseudo_hnf_full_rank)
    # starts an infinite loop.
    Q = try pseudo_hnf_full_rank(P, shape)
    catch e
      pseudo_hnf_kb(P, shape)
    end
    return Q
  end
end

@doc Markdown.doc"""
    pseudo_hnf_with_transform(P::PMat)

Transforms $P$ into pseudo-Hermite form as defined by Cohen. Essentially the
matrix part of $P$ will be upper triangular with some technical normalisation
for the off-diagonal elements. This operation preserves the module.
The used transformation is returned as a second return value.

A optional second argument can be secified as a symbols, indicating the desiered
shape of the echelon form. Possible are
`:upperright` (the default) and `:lowerleft`
"""
pseudo_hnf_with_transform(P::PMat{nf_elem, NfOrdFracIdl}, shape::Symbol = :upperright, full_rank::Bool = false) = pseudo_hnf_kb_with_transform(P, shape)

pseudo_hnf(P::PMat{T, S}, shape::Symbol = :upperright, full_rank::Bool = false) where {T <: NumFieldElem, S} = pseudo_hnf_kb(P, shape)

pseudo_hnf_with_transform(P::PMat{T, S}, shape::Symbol = :upperright, full_rank::Bool = false) where {T <: NumFieldElem, S} = pseudo_hnf_kb_with_transform(P, shape)

function pseudo_hnf_full_rank(P::PMat, shape::Symbol = :upperright)
  PP = deepcopy(P)
  K = parent(PP.matrix[1, 1])
  integralizer = _make_integral!(PP)
  m = find_pseudo_hnf_modulus(PP)
  PPhnf = pseudo_hnf_mod(PP, m, shape)
  for i in 1:nrows(PP)
    PPhnf.coeffs[i] = PPhnf.coeffs[i]*inv(K(integralizer))
    simplify(PPhnf.coeffs[i])
  end
  return PPhnf
end

function pseudo_hnf_full_rank_with_modulus!(P::PMat, m::NfOrdIdl, shape::Symbol = :upperright)
  K = parent(P.matrix[1, 1])
  integralizer = _make_integral!(P)
  m = integralizer*m
  PPhnf = pseudo_hnf_mod(P, m, shape)
  for i in 1:nrows(PPhnf)
    PPhnf.coeffs[i] = PPhnf.coeffs[i]*inv(K(integralizer))
    simplify(PPhnf.coeffs[i])
  end
  return PPhnf
end

function pseudo_hnf_full_rank_with_modulus(P::PMat, m::NfOrdIdl, shape::Symbol = :upperright)
  PP = deepcopy(P)
  return pseudo_hnf_full_rank_with_modulus!(PP, m, shape)
end

pseudo_hnf_full_rank_with_modulus(P::PMat, m::NfRelOrdIdl, shape::Symbol = :upperright) = pseudo_hnf_kb(P, shape)

function find_pseudo_hnf_modulus(P::PMat{T, S}) where {T, S}
  K = parent(P.matrix[1, 1])
  O = order(P.coeffs[1])
  if nrows(P) == ncols(P)
    m = det(P)
    simplify(m)
    return numerator(m)
  end
  p = next_prime(2^61)
  permGroup = SymmetricGroup(nrows(P))
  rowPerms = elem_type(permGroup)[]
  cnt = 0
  to_remove = Int[]
  while length(rowPerms) < 2 && cnt < min(5, ncols(P))
    cnt += 1
    lp = prime_ideals_over(O, p)
    for t in lp
      F, mF = ResidueField(O, t)
      mFF = extend(mF, K)
      Pt = zero_matrix(F, nrows(P), ncols(P))
      nextIdeal = false
      for i = 1:nrows(P)
        if i in to_remove
          continue
        end
        for j = 1:ncols(P)
          try Pt[i, j] = mFF(P.matrix[i, j])
          catch
            nextIdeal = true
            break
          end
        end
        if nextIdeal
          break
        end
      end
      if nextIdeal
        continue
      end
      rowPerm = one(permGroup)
      rank = lu!(rowPerm, Pt)
      if rank == ncols(P) && !(rowPerm in rowPerms)
        j = 1
        while rowPerm[j] > ncols(P)
          j += 1
        end
        push!(to_remove, j)
        push!(rowPerms, rowPerm)
      end
    end
    p = next_prime(p)
  end
  dets = Vector{NfOrdIdl}()
  for s = 1:length(rowPerms)
    rowPerm = rowPerms[s]
    Minor = zero_matrix(K, ncols(P), ncols(P))
    C = Vector{S}(undef, ncols(P))
    for i = 1:nrows(P)
      if rowPerm[i] > ncols(P)
        continue
      end
      for j = 1:ncols(P)
        Minor[rowPerm[i], j] = P.matrix[i, j]
      end
      C[rowPerm[i]] = P.coeffs[i]
    end
    PMinor = PseudoMatrix(Minor, C)
    m1 = det(PMinor)
    simplify(m1)
    push!(dets, numerator(m1))
  end
  return sum(dets)
end

#TODO: das kann man besser machen
function _make_integral!(P::PMat{T, S}) where {T, S}
  K = parent(P.matrix[1, 1])
  O = order(P.coeffs[1])
  integralizer = one(FlintZZ)

  for i = 1:nrows(P)
    divide_row!(P.matrix, i, K(denominator(P.coeffs[i])))
    P.coeffs[i] = P.coeffs[i]*K(denominator(P.coeffs[i]))
    simplify(P.coeffs[i])
  end

  z = one(FlintZZ)
  for i in 1:nrows(P)
    for j in 1:ncols(P)
      z = lcm(z, denominator(P.matrix[i, j], O))
    end
  end

  for i in 1:nrows(P)
    mul_row!(P.matrix, i, K(z))
  end

  return z
end

function pseudo_hnf_mod(P::PMat, m, shape::Symbol = :upperright, strategy = :split)
  O = order(m)

  t_comp_red = 0.0
  t_mod_comp = 0.0
  t_sum = 0.0
  t_div = 0.0
  t_idem = 0.0

  t_comp_red += @elapsed z = _matrix_for_reduced_span(P, m)
  @vprint :PseudoHnf 1 "Computation of reduction: $t_comp_red\n"
  t_mod_comp += @elapsed zz = strong_echelon_form(z, shape, strategy)
  @vprint :PseudoHnf 1 "Modular computation: $t_mod_comp\n"

  res_mat = zero_matrix(nf(O), nrows(P), ncols(P))
  for i in 1:nrows(P)
    for j in 1:ncols(P)
      res_mat[i, j] = lift(zz[i, j]).elem_in_nf
    end
  end

  res = PMat{nf_elem, NfOrdFracIdl}(res_mat, P.coeffs)

  shift = 0
  if shape == :lowerleft
    shift = nrows(P) - ncols(P)
  end

  for i in 1:ncols(P)
    if iszero(zz[i + shift, i])
      res.matrix[i + shift, i] = one(nf(O))
      res.coeffs[i + shift] = NfOrdFracIdl(m, fmpz(1))
    else
      o = ideal(O, lift(zz[i + shift, i]))
      t_sum += @elapsed g = o + m
      if isone(g)
        oo = o
        mm = m
      else
        t_div += @elapsed oo = divexact(o, g)
        t_div += @elapsed mm = divexact(m, g)
      end
      t_idem += @elapsed e, f = idempotents(oo, mm)
      res.coeffs[i + shift] = NfOrdFracIdl(g, fmpz(1))
      t = divexact(elem_in_nf(e), elem_in_nf(zz[i + shift, i].elem))
      mul_row!(res.matrix, i + shift, t)
      res.matrix[i + shift, i] = one(nf(O))
    end
  end

  if shape == :upperright
    t = nf(O)()
    for i = (ncols(res) - 1):-1:1
      for j = (i + 1):ncols(res)
        if iszero(res.matrix[i, j])
          continue
        end
        d = res.coeffs[j]//res.coeffs[i]
        q = mod(res.matrix[i, j], d)
        q = q - res.matrix[i, j]
        for c = j:ncols(res)
          mul!(t, q, res.matrix[j, c])
          addeq!(res.matrix[i, c], t)
        end
      end
    end
  elseif shape == :lowerleft
    t = nf(O)()
    for i = (shift + 2):nrows(res)
      for j = (i - shift - 1):-1:1
        if iszero(res.matrix[i, j])
          continue
        end
        d = res.coeffs[j + shift]//res.coeffs[i]
        q = mod(res.matrix[i, j], d)
        q = q - res.matrix[i, j]
        for c = 1:j
          mul!(t, q, res.matrix[j + shift, c])
          addeq!(res.matrix[i, c], t)
        end
      end
    end
  end

  #println("computation of ideal sum : $t_sum")
  #println("computation of ideal div : $t_div")
  #println("computation of idems     : $t_idem")

  return res
end

#this is Algorithm 4 of FH2014
# we assume that span(P) \subseteq O^r
function _matrix_for_reduced_span(P::PMat, m::NfAbsOrdIdl)
  O = order(m)
  Om, OtoOm = quo(O, m)
  z = zero_matrix(Om, nrows(P), ncols(P))
  if isone(m)
    return z
  end

  for i in 1:nrows(z)
    @vprint :PseudoHnf 4 "New row\n"
    @vtime :PseudoHnf 4 I, a = _coprime_norm_integral_ideal_class(P.coeffs[i], m)
    n = norm(I, copy = false)
    qq = Om(invmod(n, minimum(m, copy = false)))
    for j in 1:ncols(z)
      q = OtoOm(O(n*divexact(P.matrix[i, j], a)))
      z[i, j] = mul!(z[i, j], q, qq)
    end
  end
  return z
end

function _matrix_for_reduced_span(P::PMat, m::NfRelOrdIdl)
  O = order(m)
  Om, OtoOm = quo(O, m)
  z = zero_matrix(Om, nrows(P), ncols(P))
  if isone(m)
    return z
  end

  for i in 1:nrows(z)
    I, a = _coprime_norm_integral_ideal_class(P.coeffs[i], m)
    n = absolute_norm(I)
    Omn = OtoOm(O(n))
    qq = inv(Omn)
    for j in 1:ncols(z)
      @assert euclid(Omn) == 1
      q = OtoOm(O(n*divexact(P.matrix[i, j], a)))
      z[i, j] = mul!(z[i, j], q, qq)
    end
  end
  return z
end
(::NfOrdQuoRing)(a::NfOrdQuoRingElem) = a

_check(a) = a.has_coord ? dot(a.coordinates, basis(parent(a))) == a : true

function _check(m::Generic.Mat{NfOrdElem})
  for i in 1:nrows(m)
    for j in 1:ncols(m)
      if !_check(m[i, j].elem)
        println("$i $j not consistent: $(m[i, j]) : $(m[i, j].coordinates)")
        error("dasdsad")
      end
    end
  end
end

function _check(m::Generic.Mat{NfOrdQuoRingElem})
  for i in 1:nrows(m)
    for j in 1:ncols(m)
      if !_check(m[i, j].elem)
        println("$i $j not consistent: $(m[i, j].elem) : $(m[i, j].elem.coordinates)")
        error("dasdsad")
      end
    end
  end
end

function divide_row!(M::Generic.Mat{T}, i::Int, r::T) where T
  for j in 1:ncols(M)
    M[i, j] = divexact(M[i, j], r)
  end
  return nothing
end

function mul_row!(M::Generic.Mat{T}, i::Int, r::T) where T
  for j in 1:ncols(M)
    M[i, j] = mul!(M[i, j], M[i, j], r)
  end
  return nothing
end

function mul_col!(M::Generic.Mat{T}, i::Int, r::T) where T
  for j in 1:nrows(M)
    M[j, i] = mul!(M[j, i], M[j, i], r)
  end
  return nothing
end

function _contained_in_span_of_pseudohnf(v::Generic.Mat{T}, P::PMat{T, S}, shape::Symbol = :upperright) where {T, S}
  # accepts :upperright or :lowerleft for the shape of P
  (shape != :upperright && shape != :lowerleft) && error("Not yet implemented.")
  start = 1
  stop = nrows(P)
  step = 1
  if shape == :lowerleft
    start, stop = stop, start
    step = -1
  end
  w = deepcopy(v)
  for i = start:step:stop
    if !(w[1, i]//P.matrix[i, i] in P.coeffs[i])
      return false
    end
    e = w[1, i]//P.matrix[i, i]
    if shape == :upperright
      for j = i:ncols(P)
        w[1, j] = w[1, j] - e*P.matrix[i, j]
      end
    elseif shape == :lowerleft
      for j = 1:i
        w[1, j] = w[1, j] - e*P.matrix[i, j]
      end
    end
  end
  if !iszero(w)
    return false
  end
  return true
end

function _contained_in_span_of_pseudohnf(v::Generic.Mat{T}, a::S, P::PMat{T, S}, shape::Symbol = :upperright) where {T, S}
  # accepts :upperright or :lowerleft for the shape of P
  (shape != :upperright && shape != :lowerleft) && error("Not yet implemented.")
  start = 1
  stop = nrows(P)
  step = 1
  if shape == :lowerleft
    start, stop = stop, start
    step = -1
  end
  w = deepcopy(v)
  for i = start:step:stop
    if !(w[1, i]//P.matrix[i, i] in P.coeffs[i]*inv(a))
      return false
    end
    e = w[1, i]//P.matrix[i, i]
    if shape == :upperright
      for j = i:ncols(P)
        w[1, j] = w[1, j] - e*P.matrix[i, j]
      end
    elseif shape == :lowerleft
      for j = 1:i
        w[1, j] = w[1, j] - e*P.matrix[i, j]
      end
    end
  end
  if !iszero(w)
    return false
  end
  return true
end

function _spans_subset_of_pseudohnf(M::PMat{nf_elem, NfOrdFracIdl}, P::PMat{nf_elem, NfOrdFracIdl}, shape::Symbol = :upperright)
  # accepts :upperright or :lowerleft for the shape of P
  (shape != :upperright && shape != :lowerleft) && error("Not yet implemented.")
  for i in 1:nrows(M)
    A = M.coeffs[i]
    v = sub(M.matrix, i:i, 1:ncols(P))
    for b in basis(A.num)
      bb = divexact(elem_in_nf(b), A.den)
      if !Hecke._contained_in_span_of_pseudohnf(bb*v, P, shape)
        return false
      end
    end
  end
  return true
end

function _spans_subset_of_pseudohnf(M::PMat{T, S}, P::PMat{T, S}, shape::Symbol = :upperright) where {T, S}
  # accepts :upperright or :lowerleft for the shape of P
  (shape != :upperright && shape != :lowerleft) && error("Not yet implemented.")
  for i in 1:nrows(M)
    v = sub(M.matrix, i:i, 1:ncols(P))
    if !Hecke._contained_in_span_of_pseudohnf(v, M.coeffs[i], P, shape)
      return false
    end
  end
  return true
end

function sub(M::Generic.Mat, rows::UnitRange{Int}, cols::UnitRange{Int})
  @assert step(rows) == 1 && step(cols) == 1
  z = zero_matrix(base_ring(M), length(rows), length(cols))
  for i in rows
    for j in cols
      z[i - first(rows) + 1, j - first(cols) + 1] = M[i, j]
    end
  end
  return z
end

function sub(P::PMat, rows::UnitRange{Int}, cols::UnitRange{Int})
  m = sub(P.matrix, rows, cols)
  return PseudoMatrix(m, P.coeffs[rows])
end

pseudo_matrix(x...) = PseudoMatrix(x...)

function pseudo_hnf_cohen(P::PMat)
   return _pseudo_hnf_cohen(P, Val{false})
end

function pseudo_hnf_cohen_with_transform(P::PMat)
   return _pseudo_hnf_cohen(P, Val{true})
end

function _pseudo_hnf_cohen(P::PMat, trafo::Type{Val{T}} = Val{false}) where T
   H = deepcopy(P)
   m = nrows(H)
   if trafo == Val{true}
      U = identity_matrix(base_ring(H.matrix), m)
      pseudo_hnf_cohen!(H, U, true)
      return H, U
   else
      U = similar(H.matrix, 0, 0)
      pseudo_hnf_cohen!(H, U, false)
      return H
   end
end

#=
Algorithm 2.6 in "Hermite and Smith normal form algorithms over Dedekind domains"
by H. Cohen.
=#
function pseudo_hnf_cohen!(H::PMat, U::Generic.Mat{T}, with_transform::Bool = false) where T <: nf_elem
   m = nrows(H)
   n = ncols(H)
   A = H.matrix
   K = base_ring(A)
   t = K()
   t1 = K()
   t2 = K()
   k = 1
   for i = 1:n
      j = k
      while j <= m && A[j, i] == 0
         j += 1
      end
      if j > m
         continue
      end
      if j > k
         swap_rows!(H, j, k)
         with_transform ? swap_rows!(U, j, k) : nothing
      end
      H.coeffs[k] = H.coeffs[k]*A[k, i]
      simplify(H.coeffs[k])
      with_transform ? divide_row!(U, k, A[k, i]) : nothing
      divide_row!(A, k, A[k, i])
      for j = k+1:m
         if iszero(A[j, i])
            continue
         end
         Aji = deepcopy(A[j, i])
         a = H.coeffs[j]
         aa = Aji*a
         b = H.coeffs[k]
         d = aa + b
         ad = aa//d
         simplify(ad)
         bd = b//d
         simplify(bd)
         if ad.den != 1 || bd.den != 1
            error("Ideals are not integral.")
         end
         u, v = map(K, idempotents(ad.num, bd.num))
         u = divexact(u, Aji)
         for c = i:n
            t = deepcopy(A[j, c])
            mul!(t1, A[k, c], -Aji)
            addeq!(A[j, c], t1)
            mul!(t1, t, u)
            mul!(t2, A[k, c], v)
            add!(A[k, c], t1, t2)
         end
         if with_transform
            for c = 1:m
               t = deepcopy(U[j, c])
               mul!(t1, U[k, c], -Aji)
               addeq!(U[j, c], t1)
               mul!(t1, t, u)
               mul!(t2, U[k, c], v)
               add!(U[k, c], t1, t2)
            end
         end
         H.coeffs[j] = a*b//d
         simplify(H.coeffs[j])
         H.coeffs[k] = d
         simplify(H.coeffs[k])
      end
      if iszero(A[k, i])
         continue
      end
      for j = 1:k-1
         if iszero(A[j, i])
            continue
         end
         d = H.coeffs[k]//H.coeffs[j]
         q = mod(A[j, k], d)
         q = q - A[j, k]
         for c = k:n
            mul!(t, q, A[k, c])
            addeq!(A[j, c], t)
         end
         if with_transform
            for c = 1:m
               mul!(t, q, U[k, c])
               addeq!(U[j, c], t)
            end
         end
      end
      k += 1
   end
   return nothing
end

function swap_rows!(P::PMat, i::Int, j::Int)
   swap_rows!(P.matrix, i, j)
   P.coeffs[i], P.coeffs[j] = P.coeffs[j], P.coeffs[i]
   return nothing
end

function _in_span(v::Vector{nf_elem}, P::PMat)
   @assert length(v) == ncols(P)
   m = nrows(P)
   n = ncols(P)
   K = base_ring(P.matrix)
   x = zeros(K, m)
   t = K()
   k = 1
   for i = 1:n
      s = K()
      for j = 1:k-1
         mul!(t, P.matrix[j, i], x[j])
         addeq!(s, t)
      end
      s = v[i] - s
      if iszero(P.matrix[k, i])
         if !iszero(s)
            return false, x
         end
         continue
      end
      x[k] = divexact(s, P.matrix[k, i])
      if !(x[k] in P.coeffs[k])
         return false, x
      end
      if k == min(m, n)
         break
      end
      k += 1
   end
   return true, x
end

function pseudo_hnf_kb(P::PMat, shape::Symbol = :upperright)
  if shape == :lowerleft
    H = _pseudo_hnf_kb(PseudoMatrix(reverse_cols(P.matrix), P.coeffs), Val{false})
    reverse_cols!(H.matrix)
    reverse_rows!(H.matrix)
    reverse!(H.coeffs)
    return H
  elseif shape == :upperright
    return _pseudo_hnf_kb(P, Val{false})
  else
    error("Not yet implemented")
  end
end

function pseudo_hnf_kb_with_transform(P::PMat, shape::Symbol = :upperright)
  if shape == :lowerleft
    H, U = _pseudo_hnf_kb(PseudoMatrix(reverse_cols(P.matrix), P.coeffs), Val{true})
    reverse_cols!(H.matrix)
    reverse_rows!(H.matrix)
    reverse!(H.coeffs)
    reverse_rows!(U)
    return H, U
  elseif shape == :upperright
    return _pseudo_hnf_kb(P, Val{true})
  else
    error("Not yet implemented")
  end
end

function _pseudo_hnf_kb(P::PMat, trafo::Type{Val{T}} = Val{false}) where T
   H = deepcopy(P)
   m = nrows(H)
   if trafo == Val{true}
      U = identity_matrix(base_ring(H.matrix), m)
      pseudo_hnf_kb!(H, U, true)
      return H, U
   else
      U = similar(H.matrix, 0, 0)
      pseudo_hnf_kb!(H, U, false)
      return H
   end
end

function kb_search_first_pivot(H::PMat, start_element::Int = 1)
   for r = start_element:nrows(H)
      for c = start_element:ncols(H)
         if !iszero(H.matrix[r,c])
            return r, c
         end
      end
   end
   return 0, 0
end

function kb_reduce_row!(H::PMat{T, S}, U::Generic.Mat{T}, pivot::Vector{Int}, c::Int, with_transform::Bool) where {T <: NumFieldElem, S}
   r = pivot[c]
   A = H.matrix
   t = base_ring(A)()
   for i = c+1:ncols(A)
      p = pivot[i]
      if p == 0 || iszero(A[r, i])
         continue
      end
      d = H.coeffs[p]//H.coeffs[r]
      q = mod(A[r, i], d)
      q = q - A[r, i]
      for j = i:ncols(A)
         mul!(t, q, A[p,j])
         addeq!(A[r,j], t)
      end
      if with_transform
         for j = 1:ncols(U)
            mul!(t, q, U[p,j])
            addeq!(U[r,j], t)
         end
      end
   end
   return nothing
end

function kb_reduce_column!(H::PMat{T, S}, U::Generic.Mat{T}, pivot::Vector{Int}, c::Int, with_transform::Bool, start_element::Int = 1) where {T <: NumFieldElem, S}
   r = pivot[c]
   A = H.matrix
   t = base_ring(A)()
   for i = start_element:c-1
      p = pivot[i]
      if p == 0 || iszero(A[p, c])
         continue
      end
      d = H.coeffs[r]//H.coeffs[p]
      q = mod(A[p, c], d)
      q = q - A[p, c]
      for j = c:ncols(A)
         mul!(t, q, A[r,j])
         addeq!(A[p,j], t)
      end
      if with_transform
         for j = 1:ncols(U)
            mul!(t, q, U[r,j])
            addeq!(U[p,j], t)
         end
      end
   end
   return nothing
end

function kb_sort_rows!(H::PMat{T, S}, U::Generic.Mat{T}, pivot::Vector{Int}, with_transform::Bool, start_element::Int = 1) where {T <: NumFieldElem, S}
   m = nrows(H)
   n = ncols(H)
   pivot2 = zeros(Int, m)
   for i = 1:n
      if pivot[i] == 0
         continue
      end
      pivot2[pivot[i]] = i
   end

   r1 = start_element
   for i = start_element:n
      r2 = pivot[i]
      if r2 == 0
         continue
      end
      if r1 != r2
         swap_rows!(H, r1, r2)
         with_transform ? swap_rows!(U, r1, r2) : nothing
         p = pivot2[r1]
         pivot[i] = r1
         if p != 0
            pivot[p] = r2
         end
         pivot2[r1] = i
         pivot2[r2] = p
      end
      r1 += 1
      if r1 == m
         break
      end
   end
   return nothing
end

function pseudo_hnf_kb!(H::PMat{T, S}, U::Generic.Mat{T}, with_transform::Bool = false, start_element::Int = 1) where {T <: NumFieldElem, S}
   m = nrows(H)
   n = ncols(H)
   A = H.matrix
   K = base_ring(A)
   pivot = zeros(Int, n)
   row1, col1 = kb_search_first_pivot(H, start_element)
   if row1 == 0
      return nothing
   end
   pivot[col1] = row1
   pivot_max = col1
   H.coeffs[row1] = H.coeffs[row1]*A[row1, col1]
   simplify(H.coeffs[row1])
   with_transform ? divide_row!(U, row1, A[row1, col1]) : nothing
   divide_row!(A, row1, A[row1, col1])
   t = K()
   t1 = K()
   t2 = K()
   for i=row1:m-1
      new_pivot = false
      for j = start_element:pivot_max
         if iszero(A[i+1,j])
            continue
         end
         if pivot[j] == 0
            pivot[j] = i+1
            pivot_max = max(pivot_max, j)
            new_pivot = true
            H.coeffs[i+1] = H.coeffs[i+1]*A[i+1, j]
            simplify(H.coeffs[i+1])
            with_transform ? divide_row!(U, i+1, A[i+1, j]) : nothing
            divide_row!(A, i+1, A[i+1, j])
            kb_reduce_row!(H, U, pivot, j, with_transform)
         else
            p = pivot[j]
            Aij = deepcopy(A[i+1, j])
            a = H.coeffs[i+1]
            aa = Aij*a
            b = H.coeffs[p]
            d = aa + b
            ad = aa//d
            simplify(ad)
            bd = b//d
            simplify(bd)
            if typeof(ad) == NfOrdFracIdl
              if ad.den != 1 || bd.den != 1
                error("Ideals are not integral.")
              end
              u, v = map(K, idempotents(ad.num, bd.num))
            else
              if !isintegral(ad) || !isintegral(bd)
                error("Ideals are not integral.")
              end
              # numerator(ad) would make a deepcopy...
              adint = ideal_type(order(ad))(order(ad), basis_pmatrix(ad, copy = false))
              bdint = ideal_type(order(bd))(order(bd), basis_pmatrix(bd, copy = false))
              u, v = map(K, idempotents(adint, bdint))
            end
            u = divexact(u, Aij)
            for c = j:n
               t = deepcopy(A[i+1, c])
               #t1 = mul!(t1, A[p, c], -Aij)
               mul!(t1, A[p, c], -Aij)
               #A[i+1, c] = addeq!(A[i+1, c], t1)
               addeq!(A[i+1, c], t1)
               #t1 = mul!(t1, t, u)
               mul!(t1, t, u)
               #t2 = mul!(t2, A[p, c], v)
               mul!(t2, A[p, c], v)
               #A[p, c] = add!(A[p, c], t1, t2)
               add!(A[p, c], t1, t2)
            end
            if with_transform
               for c = 1:m
                  t = deepcopy(U[i+1, c])
                  #t1 = mul!(t1, U[p, c], -Aij)
                  mul!(t1, U[p, c], -Aij)
                  #U[i+1, c] = addeq!(U[i+1, c], t1)
                  addeq!(U[i+1, c], t1)
                  #t1 = mul!(t1, t, u)
                  mul!(t1, t, u)
                  #t2 = mul!(t2, U[p, c], v)
                  mul!(t2, U[p, c], v)
                  #U[p, c] = add!(U[p, c], t1, t2)
                  add!(U[p, c], t1, t2)
               end
            end
            H.coeffs[i+1] = a*b//d
            simplify(H.coeffs[i+1])
            H.coeffs[p] = d
            simplify(H.coeffs[p])
         end
         kb_reduce_column!(H, U, pivot, j, with_transform, start_element)
         if new_pivot
            break
         end
      end
      if !new_pivot
         for c = pivot_max+1:n
            if !iszero(A[i+1,c])
               pivot[c] = i+1
               pivot_max = max(pivot_max, c)
               H.coeffs[i+1] = H.coeffs[i+1]*A[i+1, c]
               simplify(H.coeffs[i+1])
               with_transform ? divide_row!(U, i+1, A[i+1, c]) : nothing
               divide_row!(A, i+1, A[i+1, c])
               kb_reduce_column!(H, U, pivot, c, with_transform, start_element)
               break
            end
         end
      end
   end
   kb_sort_rows!(H, U, pivot, with_transform, start_element)
   return nothing
end

mutable struct PMat2
   parent
   matrix::Generic.MatSpaceElem{nf_elem}
   row_coeffs::Vector{NfOrdFracIdl}
   col_coeffs::Vector{NfOrdFracIdl}

   function PMat2(m::Generic.MatSpaceElem{nf_elem}, r::Vector{NfOrdFracIdl}, c::Vector{NfOrdFracIdl})
      z = new()
      z.matrix = m
      z.row_coeffs = r
      z.col_coeffs = c
      return z
   end
end

function show(io::IO, P::PMat2)
   print(io, "Pseudo-matrix over $(parent(P.matrix[1, 1]))\n")
   print(io, "$(P.matrix)\n")
   print(io, "\nwith row ideals\n")
   for I in P.row_coeffs
      show(IOContext(io, :compact => true), I)
      print(io, "\n")
   end
   print(io, "\nand column ideals")
   for I in P.col_coeffs
      print(io, "\n")
      show(IOContext(io, :compact => true), I)
   end
end

function PseudoMatrix2(m::Generic.Mat{nf_elem}, r::Vector{NfOrdFracIdl}, c::Vector{NfOrdFracIdl})
   ( nrows(m) != length(r) || ncols(m) != length(c) ) && error("Dimensions do not match.")
   return PMat2(m, r, c)
end

function PseudoMatrix2(m::Generic.Mat{NfOrdElem}, r::Vector{NfOrdFracIdl}, c::Vector{NfOrdIdl})
   mm = change_base_ring(nf(base_ring(m)), m)
   rr = map(z -> NfOrdFracIdl(z, fmpz(1)), r)
   cc = map(z -> NfOrdFracIdl(z, fmpz(1)), c)
   return PMat(mm, rr, cc)
end

function nrows(m::PMat2)
   return nrows(m.matrix)
end

function ncols(m::PMat2)
   return ncols(m.matrix)
end

function pseudo_snf_kb(P::PMat2)
   return _pseudo_snf_kb(P, Val{false})
end

function pseudo_snf_kb_with_transform(P::PMat2)
   return _pseudo_snf_kb(P, Val{true})
end

function _pseudo_snf_kb(P::PMat2, trafo::Type{Val{T}} = Val{false}) where T
   S = deepcopy(P)
   m = nrows(S)
   n = ncols(S)
   if trafo == Val{true}
      U = identity_matrix(base_ring(S.matrix), m)
      K = identity_matrix(base_ring(S.matrix), m)
      pseudo_snf_kb!(S, U, K, true)
      return S, U, K
   else
      U = similar(S.matrix, 0, 0)
      K = U
      pseudo_snf_kb!(S, U, K, false)
      return S
   end
end

function kb_clear_row!(S::PMat2, K::Generic.Mat{nf_elem}, i::Int, with_transform::Bool)
   m = nrows(S)
   n = ncols(S)
   A = S.matrix
   t = base_ring(A)()
   t1 = base_ring(A)()
   t2 = base_ring(A)()
   for j = i+1:n
      if A[i, j] == 0
         continue
      end
      #Just for debugging, this should not happen.
      if A[i, i] != 1
         error("Pivot is not 1.")
      end
      Aij = deepcopy(A[i, j])
      a = S.col_coeffs[j]
      aa = Aij*a
      b = S.col_coeffs[i]
      d = aa + b
      ad = aa//d
      simplify(ad)
      bd = b//d
      simplify(bd)
      if ad.den != 1 || bd.den != 1
         error("Ideals are not integral.")
      end
      u, v = map(base_ring(A), idempotents(ad.num, bd.num))
      u = divexact(u, Aij)
      for r = i:m
         t = deepcopy(A[r, j])
         mul!(t1, A[r, i], -Aij)
         addeq!(A[r,j], t1)
         mul!(t1, t, u)
         mul!(t2, A[r, i], v)
         add!(A[r, i], t1, t2)
      end
      if with_transform
         for r = 1:n
            t = deepcopy(K[r, j])
            mul!(t1, K[r, i], -Aij)
            addeq!(K[r,j], t1)
            mul!(t1, t, u)
            mul!(t2, K[r, i], v)
            add!(K[r, i], t1, t2)
         end
      end
      S.col_coeffs[j] = a*b//d
      simplify(S.col_coeffs[j])
      S.col_coeffs[i] = d
      simplify(S.col_coeffs[i])
   end
   return nothing
end

function pseudo_snf_kb!(S::PMat2, U::Generic.Mat{nf_elem}, K::Generic.Mat{nf_elem}, with_transform::Bool = false)
   m = nrows(S)
   n = ncols(S)
   A = S.matrix
   l = min(m,n)
   i = 1
   t = base_ring(A)()
   t1 = base_ring(A)()
   t2 = base_ring(A)()
   H = PseudoMatrix(A, S.row_coeffs)
   if !iszero(A[1, 1])
      S.row_coeffs[1] = S.row_coeffs[1]*A[1, 1]
      simplify(S.row_coeffs[1])
      with_transform ? divide_row!(U, 1, A[1, 1]) : nothing
      divide_row!(A, 1, A[1, 1])
   end
   while i<=l
      !iszero(A[i, i]) ? kb_clear_row!(S, K, i, with_transform) : nothing
      pseudo_hnf_kb!(H, U, with_transform, i)
      c = i + 1
      while c <= n && iszero(A[i, c])
         c += 1
      end
      if c != n + 1
         continue
      end
      i += 1
   end
   return nothing
end

mutable struct ModDed
   pmatrix::PMat
   base_ring::NfOrd
   is_triu::Bool
   function ModDed(P::PMat, is_triu::Bool = false; check::Bool = true)
      if check
         is_triu = istriu(P.matrix)
      end
      z = new()
      z.pmatrix = P
      z.base_ring = order(P.coeffs[1])
      z.is_triu = is_triu
      return z
   end
end

base_ring(M::ModDed) = M.base_ring

function istriu(A::Generic.Mat)
   m = nrows(A)
   n = ncols(A)
   d = 0
   for r = 1:m
      for c = 1:n
         if !iszero(A[r, c])
            if c <= d
               return false
            end
            d = c
            break
         end
      end
   end
   return true
end

function show(io::IO, M::ModDed)
   print(io, "Module over $(M.base_ring) with defining pseudo-matrix")
   for i in 1:nrows(M.pmatrix.matrix)
      print(io, "\n")
      show(IOContext(io, :compact => true), M.pmatrix.coeffs[i])
      print(io, " with row $(sub(M.pmatrix.matrix, i:i, 1:ncols(M.pmatrix.matrix)))")
   end
end

ModDed(m::Generic.Mat{nf_elem}, is_triu::Bool = false; check::Bool = true) = ModDed(PseudoMatrix(m), is_triu; check = check)

ModDed(m::Generic.Mat{NfOrdElem}, is_triu::Bool = false; check::Bool = true) = ModDed(PseudoMatrix(m), is_triu; check = check)

ModDed(c::Vector{NfOrdFracIdl}) = ModDed(PseudoMatrix(c), true; check = false)

ModDed(c::Vector{NfOrdIdl}) = ModDed(PseudoMatrix(c), true; check = false)

function in(v::Vector{nf_elem}, M::ModDed)
   P = M.pmatrix
   if !M.is_triu
      P = pseudo_hnf_kb(M.pmatrix)
   end
   return _in_span(v, P)[1]
end

function simplify_basis!(M::ModDed)
   if !M.is_triu
      pseudo_hnf_kb!(M.pmatrix, M.pmatrix.matrix, false)
   end
   M.is_triu = true
   P = M.pmatrix
   r = nrows(P)
   for i = nrows(P):-1:1
      if !iszero_row(P.matrix, i)
         break
      end
      r -= 1
   end
   deleteat!(P.coeffs, r + 1:nrows(P))
   t = reshape(P.matrix.entries, nrows(P)*ncols(P))
   for i = nrows(P):-1:r + 1
      deleteat!(t, [ i*j for j in 1:ncols(P)])
   end
   P.matrix.entries = reshape(t, r, ncols(P))
   return nothing
end

function simplify_basis(M::ModDed)
   N = deepcopy(M)
   simplify_basis!(N)
   return N
end

function vcat(P::PMat, Q::PMat)
   @assert base_ring(P.matrix) == base_ring(Q.matrix)
   m = vcat(P.matrix, Q.matrix)
   c = vcat(P.coeffs, Q.coeffs)
   return PseudoMatrix(m, c)
end

function vcat(A::Vector{ <: PMat })
  m = vcat([ P.matrix for P in A ])
  c = vcat([ P.coeffs for P in A ]...)
  return PseudoMatrix(m, c)
end

function vcat(A::PMat...)
  m = vcat([ P.matrix for P in A ])
  c = vcat([ P.coeffs for P in A ]...)
  return PseudoMatrix(m, c)
end

function hcat(P::PMat, Q::PMat)
   @assert base_ring(P.matrix) == base_ring(Q.matrix)
   @assert P.coeffs == Q.coeffs
   m = hcat(P.matrix, Q.matrix)
   return PseudoMatrix(m, P.coeffs)
end

function hcat(A::Vector{ <: PMat })
  @assert all( P -> P.coeffs == A[1].coeffs, A)
  m = hcat([ P.matrix for P in A ])
  return PseudoMatrix(m, A[1].coeffs)
end

function hcat(A::PMat...)
  @assert all( P -> P.coeffs == A[1].coeffs, A)
  m = hcat([ P.matrix for P in A ])
  return PseudoMatrix(m, A[1].coeffs)
end

function +(M::ModDed, N::ModDed)
   @assert base_ring(M) == base_ring(N)
   m = deepcopy(M.pmatrix)
   n = deepcopy(N.pmatrix)
   mn = vcat(m, n)
   MN = ModDed(mn; check = false)
   simplify_basis!(MN)
   return MN
end

function intersect(M::ModDed, N::ModDed)
   @assert base_ring(M) == base_ring(N)
   @assert ncols(M.pmatrix) == ncols(N.pmatrix)
   MM = simplify_basis(M)
   NN = simplify_basis(N)
   A = deepcopy(MM.pmatrix)
   B = deepcopy(NN.pmatrix)
   if nrows(B) > nrows(A)
      A, B = B, A
   end
   C1 = hcat(A, deepcopy(A))
   m = zero_matrix(base_ring(B.matrix), nrows(B), ncols(B))
   C2 = hcat(B, PseudoMatrix(m, B.coeffs))
   C = vcat(C1, C2)
   # C = [A A; B 0]
   pseudo_hnf_kb!(C, C.matrix, false)
   for i = 1:nrows(B)
      for j = 1:ncols(B)
         m[i, j] = C.matrix[i + nrows(A), j + ncols(A)]
      end
   end
   D = PseudoMatrix(m, C.coeffs[nrows(A) + 1:nrows(C)])
   MN = ModDed(D, true; check = false)
   simplify_basis!(MN)
   return MN
end

function mod(M::ModDed, p::NfOrdIdl)
   O = base_ring(M)
   Op = ResidueRing(O, p)
   N = zero_matrix(Op, nrows(M.pmatrix), ncols(M.pmatrix))
   MM = M.pmatrix.matrix
   ideals = M.pmatrix.coeffs
   for i = 1:nrows(N)
      a = numerator(ideals[i]).gen_one
      if valuation(a, p) > valuation(numerator(ideals[i]).gen_two, p)
         a = numerator(ideals[i]).gen_two
      end
      for j = 1:ncols(N)
         N[i, j] = Op(O(MM[i, j]*a))
      end
   end
   return N
end

################################################################################
#
#  Print the size of a pseudo matrix
#
################################################################################

# Just for debugging
# Prints the size of the ideals/entries of a pseudo matrix
# The first column is nbits(norm(numerator)) + nbits(denominator) of the ideal
# The rest of entries are nbits(max(numerator)) + nbits(denominator)
# (The size of the entries is with respect the equation order
function size(A::PMat)
  K = parent(A.matrix[1, 1])

  println("Size is:")
  size = Array{String}(nrows(A), ncols(A) + 2)
  for i in 1:nrows(A)
    size[i, 1] = "$i"
    size[i, 2] = "$(nbits(norm(numerator(A.coeffs[i])))) $(nbits(denominator(A.coeffs[i])))"
  end
  for i in 1:nrows(A)
    for j in 1:ncols(A)
      if iszero(A.matrix[i, j])
        size[i, j + 2] = "0"
      else
        a = numerator(A.matrix[i, j])
        b = denominator(A.matrix[i, j])
        size[i, j + 2] = "$(nbits(maximum([ZZ(coeff(a, i)) for i in 0:degree(K) - 1]))) $(nbits(b))"
      end
    end
  end
  display(size)
end

function ispseudo_hnf(M, shape::Symbol = :lowerleft)
  return istriangular(M.matrix, shape)
end

function test_triangular()
  M = zero_matrix(FlintZZ, 3, 3)

  M = FlintZZ[1 0 0;
              0 1 0;
              0 0 1]

  @assert istriangular(M)

  M = FlintZZ[0 0 0;
              0 1 0;
              0 0 1]

  @assert istriangular(M)

  M = FlintZZ[1 0 0;
              0 0 0;
              0 0 1]

  @assert !istriangular(M)

  M = FlintZZ[0 1 0;
              0 0 1;
              0 0 0]

  @assert !istriangular(M)

  M = FlintZZ[1 0 0;
              0 1 0;
              0 0 0]

  @assert !istriangular(M)

  M = FlintZZ[0 1 0;
              1 0 0;
              0 0 1]

  @assert !istriangular(M)
end

function istriangular(M::MatElem, shape::Symbol = :lowerleft)
  r = nrows(M)
  c = ncols(M)

  if shape == :lowerleft

    piv = 0

    k = 1
    while iszero_row(M, k) && k <= r
      k = k + 1
    end
    if k == r + 1
      return true # Zero matrix
    end

    for i in k:r
      for j in c:-1:1
        if !iszero(M[i, j])
          if j <= piv
            return false
          else
            piv = j
            break
          end
        end
        if j == 1 && piv != 0
          piv = 0
        end
      end
      if piv > 0
        continue
      end
      return false # There should not be a zero row
    end
    return true
  elseif shape == :upperright
    return istriangular(transpose(M), :lowerleft)
  end
end

function Base.hash(P::PMat, h::UInt)
  h = Base.hash(P.matrix, h)
  return Base.hash(P.coeffs, h)
end

# Returns x in K with xa integral and coprime to m

function integral_and_coprime_to(a::NfOrdFracIdl, m::NfAbsOrdIdl)
  O = order(m)
  b = inv(a)
  B = absolute_basis(b)
  while true
    z = rand(B, -10:10)
    if iszero(z)
      continue
    end
    I = z * a
    I = simplify(I)
    @assert denominator(I) == 1
    if iscoprime(I.num, m)
      return z
    end
  end
end

function integral_and_coprime_to(a::Union{ NfOrdFracIdl, NfRelOrdFracIdl }, m::Union{ NfAbsOrdIdl, NfRelOrdIdl })
  O = order(m)

  facm = factor(m)
  faca = factor(a)

  primes = Vector{ideal_type(O)}()
  v = Vector{Int}()
  for (p, e) in faca
    if e < 0
      push!(primes, p)
      push!(v, -e)
    end
  end

  for (p, e) in facm
    if haskey(faca, p)
      if faca[p] >= 0
        push!(primes, p)
        push!(v, -faca[p])
      end
    else
      push!(primes, p)
      push!(v, fmpz(0))
    end
  end

  if isempty(primes)
    return one(nf(O))
  end

  return approximate(v, primes)
end

function steinitz_form(P::PMat)
  return _steinitz_form(P, Val{false})
end

function steinitz_form_with_transform(P::PMat)
  return _steinitz_form(P, Val{true})
end

function _steinitz_form(P::PMat, trafo::Type{Val{T}} = Val{false}) where T
  if trafo == Val{true}
    S, U = pseudo_hnf_with_transform(P, :lowerleft)
  else
    S = pseudo_hnf(P, :lowerleft)
  end

  K = base_ring(S.matrix)
  oneK = one(K)
  O = order(S.coeffs[1])
  start_row = 1
  for i = 1:nrows(S)
    if !iszero(S.matrix[i, 1])
      start_row = i
      break
    end
    S.coeffs[i] = oneK*O
  end
  if trafo == Val{true}
    steinitz_form!(S, U, true, start_row)
    return S, U
  else
    U = similar(S.matrix, 0, 0)
    steinitz_form!(S, U, false, start_row)
    return S
  end
end

# Algorithm 4.6.2 in Hoppe: Normal forms over Dedekind domains
function steinitz_form!(M::PMat{T, S}, U::Generic.Mat{T}, with_transform::Bool = false, start_row::Int = 1) where { T <: NumFieldElem, S }
  if nrows(M) < start_row
    return nothing
  end

  A = M.matrix
  K = base_ring(A)
  oneK = one(K)
  t = K()
  t1 = K()
  O = order(M.coeffs[1])
  for r = start_row:nrows(M) - 1
    a = M.coeffs[r]

    if isone(a)
      continue
    end

    if isone(M.coeffs[r + 1])
      swap_rows!(M, r, r + 1)
      continue
    end

    if a isa NfOrdFracIdl && a.num.is_principal == 1
      x = divexact(K(a.num.princ_gen), K(a.den))
      mul_row!(A, r, x)
      with_transform ? divide_row!(U, r, x) : nothing
      M.coeffs[r] = oneK*O
      continue
    end

    b = M.coeffs[r + 1]
    # Hoppe, Algorithm 1.8.5
    d = lcm(denominator(a), denominator(b))
    ad = simplify(a*d)
    bd = simplify(b*d)
    @assert denominator(ad) == 1 && denominator(bd) == 1
    ad = numerator(ad)
    bd = numerator(bd)
    iad = inv(ad)
    m1 = integral_and_coprime_to(iad, bd)
    mad = simplify(m1*iad)
    @assert denominator(mad) == 1
    x, m2 = idempotents(numerator(mad), bd)
    @assert x + m2 == 1
    m1 = divexact(m1, d)
    m2 = divexact(elem_in_nf(m2), d)
    @assert x + d * m2 == 1
    n1 = divexact(x, m1)
    @assert n1 * m1 + d * m2 == 1
    n2 = K(-1)*d
    # We now have m1 in a, m2 in b and n1 in a^-1, n2 in b^-1 with m1n1 - m2n2 = 1
    @assert m1 in a
    @assert m2 in b
    @assert n1 in inv(a)
    @assert n2 in inv(b)
    @assert m1 * n1 - m2*n2 == 1

    for c = 1:ncols(M)
      t = deepcopy(A[r, c])

      A[r, c] = mul!(A[r, c], A[r, c], m1)
      t1 = mul!(t1, A[r + 1, c], m2)
      A[r, c] = add!(A[r, c], A[r, c], t1)

      t1 = mul!(t1, t, n2)
      A[r + 1, c] = mul!(A[r + 1, c], A[r + 1, c], n1)
      A[r + 1, c] = add!(A[r + 1, c], A[r + 1, c], t1)
    end

    if with_transform
      for c = 1:ncols(U)
        t = deepcopy(U[r, c])

        U[r, c] = mul!(U[r, c], U[r, c], m1)
        t1 = mul!(t1, U[r + 1, c], m2)
        U[r, c] = add!(U[r, c], U[r, c], t1)

        t1 = mul!(t1, t, n2)
        U[r + 1, c] = mul!(U[r + 1, c], U[r + 1, c], n1)
        U[r + 1, c] = add!(U[r + 1, c], U[r + 1, c], t1)
      end
    end

    M.coeffs[r] = oneK*O
    M.coeffs[r + 1] = a*b
  end
  return nothing
end

################################################################################
#
#  Inverse
#
################################################################################

function inv(x::Generic.MatSpaceElem{NfOrdElem})
  R = base_ring(x)
  K = nf(R)
  return change_base_ring(R, inv(change_base_ring(K, x)))
end


text:
using Logging
using Revise
# https://github.com/timholy/CodeTracking.jl
using CodeTracking
# see also Sugar.jl
# https://discourse.julialang.org/t/source-code-of-a-function/4499/4

"""
    call_1(func, x...)

Call the argument function.  This is for testing delaying calling a function.

...
# Arguments
- `func`: function
...
"""
function call_1(func, x...)
    @info "call_1"
    func(x...)
end

"""
Use a Dict to pass array of functions to call later.

...
# Arguments
- `db`: Dict of functions
...
"""
function call_2(db, x...)
    @info "call_2"
    func = db["func_1"]
    func(x...)
end

function func_1()
    @info "func_1"
end

function func_2(a,b)
    # Does comment included in code_string?
    # Yes, it does.
    @info "func_2: $a, $b"
end

# call_1(func_1)
# call_1(func_2, 13, 17)

db = Dict()
db["func_1"] = func_1
call_2(db)
#db["func_1"]()

code = @code_string(func_2(1,2))
@info code
@which call_2(Dict())

ext: jl


text:
type Network
  pairs::Vector{Tuple{Int,Int}}
  pair_idx::Dict{Tuple{Int,Int}, Int}
  np::Int

  buses::Vector{Int}
  bus_idx::Dict{Int, Int}
  nb::Int

  from::Vector{Int}
  to::Vector{Int}

  function Network(pairs::Vector{Tuple{Int,Int}}, buses::Vector{Int})
    pair_idx = Dict(p=>j for (j,p) in enumerate(pairs))
    np = length(pairs)

    bus_idx = Dict(b=>j for (j,b) in enumerate(buses))
    nb = length(buses)

    from  = [bus_idx[p[1]] for p in pairs]
    to  =   [bus_idx[p[2]] for p in pairs]

    new(pairs, pair_idx, np, buses, bus_idx, nb, from, to)
  end
end

function Network(pm::GenericPowerModel)
  pairs = [k for k in keys(pm.ref[:buspairs])]
  buses = [k for k in keys(pm.ref[:bus])]

  Network(pairs, buses)
end

function calc_YYY(pm::GenericPowerModel, net::Network)
  Yᵈ = zeros(Complex, net.nb);
  Yᶠ = zeros(Complex, net.np);
  Yᵗ = zeros(Complex, net.np);

  for (br, br_data) in pm.ref[:branch]
    g, b = PowerModels.calc_branch_y(br_data)
    tr, ti = PowerModels.calc_branch_t(br_data)

    Ys = br_data["br_status"]*(g + im*b)
    Ytt = Ys + im*br_data["br_status"]*br_data["br_b"]/2

    br_ft = (br_data["f_bus"], br_data["t_bus"])
    f_idx = net.bus_idx[br_ft[1]]
    t_idx = net.bus_idx[br_ft[2]]

    Yᵈ[f_idx] = Yᵈ[f_idx] + Ytt/br_data["tap"]^2
    Yᵈ[t_idx] = Yᵈ[t_idx] + Ytt

    if haskey(net.pair_idx, br_ft)
      br_idx = net.pair_idx[br_ft]
      Yᶠ[br_idx] = Yᶠ[br_idx] - Ys/(tr - im*ti)
      Yᵗ[br_idx] = Yᵗ[br_idx] - Ys/(tr + im*ti)
    elseif haskey(net.pair_idx, reverse(br_ft))
      br_idx = net.pair_idx[reverse(br_ft)]
      Yᵗ[br_idx] = Yᵗ[br_idx] - Ys/(tr - im*ti)
      Yᶠ[br_idx] = Yᶠ[br_idx] - Ys/(tr + im*ti)
    else
      ErrorException("Inconsistent pair and branch data")
    end
  end

  for (bus, bus_data) in pm.ref[:bus]
    Yᵈ[net.bus_idx[bus]] = Yᵈ[net.bus_idx[bus]] + bus_data["gs"] + im*bus_data["bs"]
  end

  Yᵈ = spdiagm(Yᵈ, (0))
  Yᶠ = sparse(net.from, 1:net.np, Yᶠ, net.nb, net.np)
  Yᵗ = sparse(net.to,   1:net.np, Yᵗ, net.nb, net.np)

  return (Yᵈ, Yᶠ, Yᵗ)
end

function get_bus_types(pm::GenericPowerModel, net::Network)
  pq    = find(b -> pm.ref[:bus][b]["bus_type"] == 1, net.buses)
  pv    = find(b -> pm.ref[:bus][b]["bus_type"] == 2, net.buses)
  slack = find(b -> pm.ref[:bus][b]["bus_type"] == 3, net.buses)

  return (pq, pv, slack)
end

function get_pq(pm::GenericPowerModel, net::Network)
  p = [-pm.ref[:bus][idx]["pd"]  for idx in net.buses]
  q = [-pm.ref[:bus][idx]["qd"]  for idx in net.buses]

  for data in values(pm.ref[:gen])
    idx = net.bus_idx[data["gen_bus"]]
    p[idx] = p[idx] + data["gen_status"]*data["pg"]
    q[idx] = q[idx] + data["gen_status"]*data["qg"]
  end
  return (p, q)
end

function get_vϑ(pm::GenericPowerModel, net::Network)
  v = [pm.ref[:bus][idx]["vm"]        for idx in net.buses]
  ϑ = [pm.ref[:bus][idx]["va"]*pi/180 for idx in net.buses]
  return (v, ϑ)
end

function get_sol_vϑ(sol::Dict, net::Network)
  v = [sol["solution"]["bus"][string(idx)]["vm"]        for idx in net.buses]
  ϑ = [sol["solution"]["bus"][string(idx)]["va"]*pi/180 for idx in net.buses]
  return (v, ϑ)
end

# function zero_loading(pm::GenericPowerModel)
#     pm0 = pm;
#     allbuses = keys(pm.data["bus"])
#     allgens = keys(pm.data["gen"])
#
#     for idx in allbuses
#       pm0.data["bus"][idx]["pd"] = 0.01*pm.data["bus"][idx]["pd"]
#       pm0.data["bus"][idx]["qd"] = 0.01*pm.data["bus"][idx]["qd"]
#     end
#     for idx in allgens
#       pm0.data["gen"][idx]["pg"] = 0.01*pm.data["gen"][idx]["pg"]
#     end
#    return pm0
# end

function zero_loading(case::Dict)
    case0 = case;
    allbuses = keys(case["bus"])
    allgens = keys(case["gen"])

    for idx in allbuses
      case0["bus"][idx]["pd"] = 0.01*case["bus"][idx]["pd"]
      case0["bus"][idx]["qd"] = 0.01*case["bus"][idx]["qd"]
    end
    for idx in allgens
      case0["gen"][idx]["pg"] = 0.01*case["gen"][idx]["pg"]
    end
   return case0
end

export Network, calc_YYY, get_bus_types, get_pq, get_vϑ, get_sol_vϑ, zero_loading

max_stars_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_issues_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_forks_repo_name: hunghtd/PowerFlowAnalysis-tutorial


text:
using Colors
using Interact

myslider = slider(0:0.01:1, label="myslider")
cmap=map(x->RGB(x,0.5,0.5), observe(myslider))
map(display,[myslider,cmap])

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using DataFrames
using Gadfly
using CSV
using StatsBase
using ColorTypes

# Code taken from https://julialang.org/benchmarks/

benchmarks = CSV.read("benchmarks.csv", DataFrame; header=["language", "benchmark", "time"])

# Capitalize and decorate language names from datafile
dict = Dict("c"=>"C", "julia"=>"Julia", "lua"=>"LuaJIT", "fortran"=>"Fortran", "java"=>"Java",
    "javascript"=>"JavaScript", "matlab"=>"Matlab", "mathematica"=>"Mathematica",
    "python"=>"Python", "octave"=>"Octave", "r"=>"R", "rust"=>"Rust", "go"=>"Go");
benchmarks[!,:language] = [dict[lang] for lang in benchmarks[!,:language]]

# Normalize benchmark times by C times
ctime = benchmarks[benchmarks[!,:language].== "C", :]
benchmarks = innerjoin(benchmarks, ctime, on=:benchmark, makeunique=true)
select!(benchmarks, Not(:language_1)) # delete column language_1
rename!(benchmarks, :time_1 =>:ctime)
benchmarks[!,:normtime] = benchmarks[!,:time] ./ benchmarks[!,:ctime];

# Compute the geometric mean for each language
langs = [];
means = [];
priorities = [];
for lang in values(dict)
    data = benchmarks[benchmarks[!,:language].== lang, :]
    gmean = geomean(data[!,:normtime])
    push!(langs, lang)
    push!(means, gmean)
    if (lang == "C")
        push!(priorities, 1)
    elseif (lang == "Julia")
        push!(priorities, 2)
    else
        push!(priorities, 3)
    end
end

# Add the geometric means back into the benchmarks dataframe
langmean = DataFrame(language=langs, geomean = means, priority = priorities)
benchmarks = innerjoin(benchmarks, langmean, on=:language)

# Put C first, Julia second, and sort the rest by geometric mean
sort!(benchmarks, [:priority, :geomean]);
sort!(langmean,   [:priority, :geomean]);

p = plot(benchmarks,
    x = :language,
    y = :normtime,
    color = :benchmark,
    Scale.y_log10,
    Guide.ylabel(nothing),
    Guide.xlabel(nothing, :vertical),
    Coord.Cartesian(xmin=1,xmax=13.3,ymin=-0.5,ymax=4.2),
    Theme(
        guide_title_position = :left,
        colorkey_swatch_shape = :circle,
        minor_label_color = colorant"white",
        major_label_color = colorant"white",
        key_label_color = colorant"white",
        key_title_color = colorant"white",
        minor_label_font_size = 12pt,
        major_label_font_size = 12pt,
        key_label_font_size = 12pt,
        key_title_font_size = 13pt,
        highlight_width = 0pt
    ),
)

golden = (1+√5)/2
draw(SVG(8inch,8inch/golden), p)
draw(SVG("benchmarks_dark.svg", 8inch, 8inch/golden), p)
draw(PNG("../images/benchmark.png", 9inch, 9inch/golden), p)

max_stars_repo_name: csimal/Julia-Unamur

max_issues_repo_name: csimal/Julia-Unamur

max_forks_repo_name: csimal/Julia-Unamur


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DifferentialEquations
using DynamicalSystems
using DelayEmbeddings
using DelimitedFiles
using Random

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we look at a chaotic system: The Duffing attractor in the funnel regime

## Integrate system and determine theiler window
Random.seed!(234)
# set time interval for integration
N = 5000 # number of samples
transients = 2000
dt_tra = 0.1
tspan = (dt_tra*N) + (dt_tra*transients)

function duffing!(du,u,p,t)
    x,y,z = u
    μ,α,β,γ,ω = p

    du[1] = y # dx
    du[2] = μ*(1-x^2)*y - α*x - β*x^3 + z
    du[3] = γ*ω*cos(ω*t)# dz
end

function prob_func_duf(prob,i,repeat)
  @. prob.u0 = [.2*rand(1); .2*rand(1); .2*rand(1)]
  prob
end

# set parameters for regular Duffing system
μ = .1
α = 1
β = 0
γ = .5
ω = 2
p = [μ,α,β,γ,ω]

# random initial condition
u0 = [.2*rand(1); .2*rand(1); .2*rand(1)]
duf = ODEProblem(duffing!,u0,tspan,p)

duf_ensemble = EnsembleProblem(duf; prob_func = prob_func_duf)
sim = solve(duf_ensemble; saveat = dt_tra, trajectories = 1000)

## Analysis of ensemble

# method for standard τ-estimation
method = "mi_min"
#method = "ac_zero"

# noise levels
σs = [0, .1]

#preallocation
L_ref = zeros(length(σs),1000)
L_tde = zeros(length(σs),1000)
L_mdop = zeros(length(σs),1000)
L_GA = zeros(length(σs),1000)
L_pec = zeros(length(σs),1000)

dim_tde = zeros(length(σs),1000)
dim_mdop = zeros(length(σs),1000)
dim_GA = zeros(length(σs),1000)
dim_pec = zeros(length(σs),1000)

mfnn_tde = zeros(length(σs),1000)
mfnn_mdop = zeros(length(σs),1000)
mfnn_GA = zeros(length(σs),1000)
mfnn_pec = zeros(length(σs),1000)

jrrf_tde = zeros(length(σs),1000)
jrrf_mdop = zeros(length(σs),1000)
jrrf_GA = zeros(length(σs),1000)
jrrf_pec = zeros(length(σs),1000)

DET_ref = zeros(length(σs),1000)
DET_tde = zeros(length(σs),1000)
DET_mdop = zeros(length(σs),1000)
DET_GA = zeros(length(σs),1000)
DET_pec = zeros(length(σs),1000)

ENTR_ref = zeros(length(σs),1000)
ENTR_tde = zeros(length(σs),1000)
ENTR_mdop = zeros(length(σs),1000)
ENTR_GA = zeros(length(σs),1000)
ENTR_pec = zeros(length(σs),1000)

RTE_ref = zeros(length(σs),1000)
RTE_tde = zeros(length(σs),1000)
RTE_mdop = zeros(length(σs),1000)
RTE_GA = zeros(length(σs),1000)
RTE_pec = zeros(length(σs),1000)

LAM_ref = zeros(length(σs),1000)
LAM_tde = zeros(length(σs),1000)
LAM_mdop = zeros(length(σs),1000)
LAM_GA = zeros(length(σs),1000)
LAM_pec = zeros(length(σs),1000)

TRANS_ref = zeros(length(σs),1000)
TRANS_tde = zeros(length(σs),1000)
TRANS_mdop = zeros(length(σs),1000)
TRANS_GA = zeros(length(σs),1000)
TRANS_pec = zeros(length(σs),1000)

@time for (cnt,σ) in enumerate(σs)
    for i = 1:1000
        display("σ=$σ, i=$i")
        tr = Array(sim[i])'[transients+1:end,:].+σ*randn(5001,3)
        tr = regularize(Dataset(tr))

        if i == 1
            w1 = estimate_delay(tr[:,1], method)
            w2 = estimate_delay(tr[:,2], method)
            global w = maximum(hcat(w1,w2))
            global taus = 0:100

            tw1 = mdop_maximum_delay(tr[:,1])
            tw2 = mdop_maximum_delay(tr[:,2])
            lm1 = DelayEmbeddings.findlocalminima(tw1[2])
            if lm1[1] ≤ 2
                lmm1 = try lm1[2]
                catch
                    lm1[1]
                end
            else
                lmm1 = lm1[1]
            end
            lm2 = DelayEmbeddings.findlocalminima(tw2[2])
            if lm2[1] ≤ 2
                lmm2 = try lm2[2]
                catch
                    lm2[1]
                end
            else
                lmm2 = lm2[1]
            end
            global tau_max = maximum(hcat(lmm1, lmm2))
            global taus_mdop = 0:tau_max
        end

        ## Perform reconstructions and compute according L-value

        #Standard TDE
        Y_tde, τ_tde, _ = optimal_traditional_de(tr[:,1]; τs = taus)
        dim_tde[cnt,i] = size(Y_tde,2)
        τ_tdes = [(i-1)*τ_tde for i = 1:(size(Y_tde,2))]
        L_tde[cnt,i] = compute_delta_L(tr[:,1], τ_tdes, taus[end]; w = w)

        #MDOP
        Y_mdop, τ_vals_mdop, ts_vals_mdop, FNNs_mdop , βs_mdop = mdop_embedding(tr[:,1];
                                                            τs = taus_mdop , w = w)
        dim_mdop[cnt,i] = size(Y_mdop,2)
        L_mdop[cnt,i] = compute_delta_L(tr[:,1], τ_vals_mdop, taus_mdop[end]; w = w)

        #Garcia & Almeida
        Y_GA, τ_vals_GA, ts_vals_GA, FNNs_GA , ns_GA = garcia_almeida_embedding(tr[:,1];
                                                            τs = taus , w = w, T = w)
        dim_GA[cnt,i] = size(Y_GA,2)
        L_GA[cnt,i] = compute_delta_L(tr[:,1], τ_vals_GA, taus[end]; w = w)

        #Pecuzal
        Y_pec, τ_vals_pec, ts_vals_pec, Ls_pec , εs_pec = pecuzal_embedding_update(tr[:,1];
                                                                    τs = taus , w = w)
        dim_pec[cnt,i] = size(Y_pec,2)
        L_pec[cnt,i] = sum(Ls_pec)

        ## compute other evaluation statistics

        mfnn_tde[cnt,i], mfnn_mdop[cnt,i], mfnn_GA[cnt,i], mfnn_pec[cnt,i],
        jrrf_tde[cnt,i], jrrf_mdop[cnt,i], jrrf_GA[cnt,i], jrrf_pec[cnt,i],
        RQA_ref, RQA1, RQA2, RQA3, RQA4, _, _, _, _, _ =
                perform_recurrence_analysis(tr, Dataset(Y_tde), Dataset(Y_mdop),
                            Dataset(Y_GA), Dataset(Y_pec); ε = 0.08, w = w, kNN = 10)

        DET_ref[cnt,i] = RQA_ref.DET
        RTE_ref[cnt,i] = RQA_ref.RTE
        ENTR_ref[cnt,i] = RQA_ref.ENTR
        LAM_ref[cnt,i] = RQA_ref.LAM
        TRANS_ref[cnt,i] = RQA_ref.TRANS

        DET_tde[cnt,i] = RQA1.DET
        RTE_tde[cnt,i] = RQA1.RTE
        ENTR_tde[cnt,i] = RQA1.ENTR
        LAM_tde[cnt,i] = RQA1.LAM
        TRANS_tde[cnt,i] = RQA1.TRANS

        DET_mdop[cnt,i] = RQA2.DET
        RTE_mdop[cnt,i] = RQA2.RTE
        ENTR_mdop[cnt,i] = RQA2.ENTR
        LAM_mdop[cnt,i] = RQA2.LAM
        TRANS_mdop[cnt,i] = RQA2.TRANS

        DET_GA[cnt,i] = RQA3.DET
        RTE_GA[cnt,i] = RQA3.RTE
        ENTR_GA[cnt,i] = RQA3.ENTR
        LAM_GA[cnt,i] = RQA3.LAM
        TRANS_GA[cnt,i] = RQA3.TRANS

        DET_pec[cnt,i] = RQA4.DET
        RTE_pec[cnt,i] = RQA4.RTE
        ENTR_pec[cnt,i] = RQA4.ENTR
        LAM_pec[cnt,i] = RQA4.LAM
        TRANS_pec[cnt,i] = RQA4.TRANS
    end
end

# save results
writedlm("./scripts/Duffing/ensemble results uni/L_ref.csv", L_ref)
writedlm("./scripts/Duffing/ensemble results uni/L_tde.csv", L_tde)
writedlm("./scripts/Duffing/ensemble results uni/L_mdop.csv", L_mdop)
writedlm("./scripts/Duffing/ensemble results uni/L_GA.csv", L_GA)
writedlm("./scripts/Duffing/ensemble results uni/L_pec.csv", L_pec)

writedlm("./scripts/Duffing/ensemble results uni/dim_tde.csv", dim_tde)
writedlm("./scripts/Duffing/ensemble results uni/dim_mdop.csv", dim_mdop)
writedlm("./scripts/Duffing/ensemble results uni/dim_GA.csv", dim_GA)
writedlm("./scripts/Duffing/ensemble results uni/dim_pec.csv", dim_pec)

writedlm("./scripts/Duffing/ensemble results uni/mfnn_tde.csv", mfnn_tde)
writedlm("./scripts/Duffing/ensemble results uni/mfnn_mdop.csv", mfnn_mdop)
writedlm("./scripts/Duffing/ensemble results uni/mfnn_GA.csv", mfnn_GA)
writedlm("./scripts/Duffing/ensemble results uni/mfnn_pec.csv", mfnn_pec)

writedlm("./scripts/Duffing/ensemble results uni/jrrf_tde.csv", jrrf_tde)
writedlm("./scripts/Duffing/ensemble results uni/jrrf_mdop.csv", jrrf_mdop)
writedlm("./scripts/Duffing/ensemble results uni/jrrf_GA.csv", jrrf_GA)
writedlm("./scripts/Duffing/ensemble results uni/jrrf_pec.csv", jrrf_pec)

writedlm("./scripts/Duffing/ensemble results uni/DET_ref.csv", DET_ref)
writedlm("./scripts/Duffing/ensemble results uni/DET_tde.csv", DET_tde)
writedlm("./scripts/Duffing/ensemble results uni/DET_mdop.csv", DET_mdop)
writedlm("./scripts/Duffing/ensemble results uni/DET_GA.csv", DET_GA)
writedlm("./scripts/Duffing/ensemble results uni/DET_pec.csv", DET_pec)

writedlm("./scripts/Duffing/ensemble results uni/ENTR_ref.csv", ENTR_ref)
writedlm("./scripts/Duffing/ensemble results uni/ENTR_tde.csv", ENTR_tde)
writedlm("./scripts/Duffing/ensemble results uni/ENTR_mdop.csv", ENTR_mdop)
writedlm("./scripts/Duffing/ensemble results uni/ENTR_GA.csv", ENTR_GA)
writedlm("./scripts/Duffing/ensemble results uni/ENTR_pec.csv", ENTR_pec)

writedlm("./scripts/Duffing/ensemble results uni/RTE_ref.csv", RTE_ref)
writedlm("./scripts/Duffing/ensemble results uni/RTE_tde.csv", RTE_tde)
writedlm("./scripts/Duffing/ensemble results uni/RTE_mdop.csv", RTE_mdop)
writedlm("./scripts/Duffing/ensemble results uni/RTE_GA.csv", RTE_GA)
writedlm("./scripts/Duffing/ensemble results uni/RTE_pec.csv", RTE_pec)

writedlm("./scripts/Duffing/ensemble results uni/LAM_ref.csv", LAM_ref)
writedlm("./scripts/Duffing/ensemble results uni/LAM_tde.csv", LAM_tde)
writedlm("./scripts/Duffing/ensemble results uni/LAM_mdop.csv", LAM_mdop)
writedlm("./scripts/Duffing/ensemble results uni/LAM_GA.csv", LAM_GA)
writedlm("./scripts/Duffing/ensemble results uni/LAM_pec.csv", LAM_pec)

writedlm("./scripts/Duffing/ensemble results uni/TRANS_ref.csv", TRANS_ref)
writedlm("./scripts/Duffing/ensemble results uni/TRANS_tde.csv", TRANS_tde)
writedlm("./scripts/Duffing/ensemble results uni/TRANS_mdop.csv", TRANS_mdop)
writedlm("./scripts/Duffing/ensemble results uni/TRANS_GA.csv", TRANS_GA)
writedlm("./scripts/Duffing/ensemble results uni/TRANS_pec.csv", TRANS_pec)


text:
"""
    MomentsAttribute(N)
    MomentsAttribute()

A constraint attribute for the vector of moments corresponding to the
constraint that a polynomial is zero in the full space in result `N`. If the
polynomial is constrained to be zero in an algebraic set, it is the moments
for the constraint once it is rewritten into an constraint on the full space.
If `N` is omitted, it is 1 by default.

## Examples

Consider the following program:
```julia
@variable(model, α)
@variable(model, β ≤ 1)

using DynamicPolynomials
@polyvar x y
cref = @constraint(model, α * x - β * y == 0, domain = @set x == y)
```
The constraint is equivalent to
```julia
@constraint(model, (α - β) * y == 0)
```
for which the dual is the 1-element vector with the moment of `y` of value `-1`.
This is the result of `moments(cref)`. However, the dual of `cref` obtained by
`dual(cref)` is the 2-elements vector with both the moments of `x` and `y`
of value `-1`.
"""
struct MomentsAttribute <: MOI.AbstractConstraintAttribute
    N::Int
end
MomentsAttribute() = MomentsAttribute(1)

function MOI.is_set_by_optimize(::MomentsAttribute)
    return true
end

# If a variable is bridged, the `VectorOfVariables`-in-`ZeroPolynomialSet` is
# bridged by `MOI.Bridges.Constraint.VectorFunctionizeBridge` and it has
# to pass the constraint to the PolyJuMP bridge.
MOI.Bridges.Constraint.invariant_under_function_conversion(::MomentsAttribute) = true

# This is type piracy but we tolerate it.
const ObjectWithoutIndex = Union{
    AbstractMonomial, AbstractTerm{<:MOI.Utilities.ObjectWithoutIndex},
    AbstractPolynomial{<:MOI.Utilities.ObjectWithoutIndex},
    MB.AbstractPolynomialBasis,
    MultivariateMoments.AbstractMomentMatrix{<:MOI.Utilities.ObjectWithoutIndex},
    MultivariateMoments.AbstractMeasure{<:MOI.Utilities.ObjectWithoutIndex}
}
const ObjectOrTupleWithoutIndex = Union{ObjectWithoutIndex, Tuple{Vararg{ObjectWithoutIndex}}}
const ObjectOrTupleOrArrayWithoutIndex = Union{ObjectOrTupleWithoutIndex, AbstractArray{<:ObjectOrTupleWithoutIndex}}
MOI.Utilities.map_indices(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x
MOI.Utilities.substitute_variables(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x


text:
"""
    eom_linearmap1(x, p, n) -> SVector{2}

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended
Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function eom_linearmap1(x, p, t)
    c = p[1]
    x, y = (x...,)
    dx = 3.4*x*(t - 1)* (1 - x^2*(t - 1)) * exp((-x^2)*(t - 1)) + 0.8*x*(t - 2)
    dy = 3.4*y*(t - 1)* (1 - y^2*(t - 1)) * exp((-y^2)*(t - 1)) + 0.5*y*(t - 2) + c*x*(t - 2)
    return SVector{2}(dx, dy)
end

"""
    linearmap1(u₀, c) -> DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function linearmap1(u₀, c)
    p = [c]
    logistic_system = DiscreteDynamicalSystem(eom_linearmap1, u₀, p)
    return logistic_system
end

"""
    linearmap1(;u₀ = [1, rand(2)], c = 2.0) -> DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
linearmap1(;u₀ = rand(2), c = 0.5) = linearmap1(u₀, c)

"""
    eom_linearmap2(x, p, n) -> SVector{2}

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended
Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function eom_chaoticmaps(x, p, t)
    c = p[1]
    x, y, z = (x...)
    t = t + 2
    dx = 3.4*x*(t - 1) * (1 - x^2*(t - 1)) * exp((-x^2)*(t - 1))
    dy = 3.4*y*(t - 1) * (1 - y^2*(t - 1)) * exp((-y^2)*(t - 1)) + 0.5*x*(t - 1)
    dz = 3.4*z*(t - 1) * (1 - z^2*(t - 1)) * exp((-z^2)*(t - 1)) + c*x*(t - 1) + 0.3*y*(t - 1)

    return SVector{3}(dx, dy, dz)
end

doc"""
    linearmap2(u₀, c) -> DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function chaoticmaps(u₀, c)
    p = [c]
    cmap = DiscreteDynamicalSystem(eom_chaoticmaps, u₀, p)
    return cmap
end

doc"""
    linearmap2(;u₀ = rand(3), c = 2.0) -> DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
chaoticmaps(;u₀ = rand(3), c = 0.5) = chaoticmaps(u₀, c)


text:
abstract type AbstractContainer end
abstract type AbstractNode end
abstract type AbstractNeurons <: AbstractNode end
abstract type AbstractSynapses <: AbstractNode end


text:
"""
This file is inspired from:
https://github.com/adap/flower/blob/main/src/py/flwr/server/client_manager.py
"""

import Base.length
using StatsBase: sample

"""
    ClientManager()

Structure which manages the clients registration and syncrhonizes the server to
make it work with them.
"""
struct ClientManager
    clients::Set{String}
    cond::Threads.Condition

    ClientManager() = begin
        clients = Set{String}()
        cond = Threads.Condition()
        return new(clients, cond)
    end
end

"""
    length(cm::ClientManager)

Return the number of client managed by the `ClientManager`.
"""
length(cm::ClientManager) = length(cm.clients)

"""
    wait_for(cm::ClientManager, num_clients::Int)

Block until at least `num_clients` are available.
"""
function wait_for(cm::ClientManager, num_clients::Int)
    lock(cm.cond) do
        while length(cm.clients) < num_clients
            wait(cm.cond)
        end
    end
end

"""
    num_available_clients(cm::ClientManager)

Return the number of available clients.
"""
num_available_clients(cm::ClientManager) = length(cm)

"""
    register!(cm::ClientManager, client_url::String)

Register a new `client_url`.

**Example**
```julia
register!(cm, "http://127.0.0.1:8081")
```
"""
function register!(cm::ClientManager, client_url::String)
    push!(cm.clients, client_url)
    lock(cm.cond) do
        notify(cm.cond)
    end
end

"""
    sample_clients(cm::ClientManager, fraction::Float32)

Sample randomly `fraction`% of the clients, without replacement. Return at
least one client if the number of available clients > 0.
"""
function sample_clients(cm::ClientManager, fraction::Float32)::Vector{String}
    length(cm) == 0 && return Vector{String}(undef, 0)
    num_clients = max(round(Int, fraction * num_available_clients(cm)), 1)
    return sample(collect(cm.clients), num_clients, replace=false)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 98601dc8-fd46-11ea-2560-e762bfd97ed7
using Pkg, DrWatson

# ╔═╡ 986050d6-fd46-11ea-26b6-7f618638f1ab
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanQuap
	using StatisticalRethinking
end

# ╔═╡ f7be0558-fd45-11ea-2cb4-c9f411d2e55e
md"## Clip-05-25-27s.jl"

# ╔═╡ 986d9d7c-fd46-11ea-305a-915f690d446a
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame)
	scale!(df, [:Marriage, :MedianAgeMarriage, :Divorce])
end;

# ╔═╡ 946d9a5a-833e-11eb-2c8e-452d18e7c2db
stan5_3_A = "
data {
  int N;
  vector[N] divorce_s;
  vector[N] marriage_s;
  vector[N] medianagemarriage_s;
}
parameters {
  real a;
  real bA;
  real bM;
  real aM;
  real bAM;
  real<lower=0> sigma;
  real<lower=0> sigma_M;
}
model {
  // A -> D <- M
  vector[N] mu = a + bA * medianagemarriage_s + bM * marriage_s;
  a ~ normal( 0 , 0.2 );
  bA ~ normal( 0 , 0.5 );
  bM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  divorce_s ~ normal( mu , sigma );
  // A -> M
  vector[N] mu_M = aM + bAM * medianagemarriage_s;
  aM ~ normal( 0 , 0.2 );
  bAM ~ normal( 0 , 0.5 );
  sigma_M ~ exponential( 1 );
  marriage_s ~ normal( mu_M , sigma_M );
}
";

# ╔═╡ 946dcbb0-833e-11eb-25c9-af83df7acab1
begin
	data = (N = size(df, 1), divorce_s = df.Divorce_s,
	  marriage_s = df.Marriage_s, medianagemarriage_s = df.MedianAgeMarriage_s)
	init = (a = 0.0, bM = 0.0, bA = -1.0, sigma = 1.0,
		aM = 0.0, bAM = -1.0, sigma_m = 1.0)
	q5_3_As, m5_3_As, o5_3_As = stan_quap("m5.3_A", stan5_3_A; data, init)
	if !isnothing(q5_3_As)
		quap5_3_As_df = sample(q5_3_As)
	end
	if !isnothing(m5_3_As)
	  post5_3_As_df = read_samples(m5_3_As, :dataframe)
	  PRECIS(post5_3_As_df)
	end
end

# ╔═╡ 946eb598-833e-11eb-2df5-5f66fa9a18d7
# Rethinking results
rethinking_results = "
           mean   sd  5.5% 94.5%
  a        0.00 0.10 -0.16  0.16
  bM      -0.07 0.15 -0.31  0.18
  bA      -0.61 0.15 -0.85 -0.37
  sigma    0.79 0.08  0.66  0.91
  aM       0.00 0.09 -0.14  0.14
  bAM     -0.69 0.10 -0.85 -0.54
  sigma_M  0.68 0.07  0.57  0.79
";

# ╔═╡ 948f1e3c-833e-11eb-391d-6b4ec346d46c
if !isnothing(q5_3_As)
		PRECIS(quap5_3_As_df)
	end

# ╔═╡ 9877c52c-fd46-11ea-293f-3d07c0cd6734
md"##### Rethinking results"

# ╔═╡ b9af5a04-fd49-11ea-07ba-fb61574aed90
part5_3_As = read_samples(m5_3_As, :particles)

# ╔═╡ 988ca05a-fd46-11ea-2ae3-910f7baa2b3f
md"## Snippet 5.25"

# ╔═╡ c79668d4-fd48-11ea-0cea-838eb6be744c
a_seq = range(-2, stop=2, length=100)

# ╔═╡ 988d33a8-fd46-11ea-27e5-7ba54e7b04fa
md"## Snippet 5.26"

# ╔═╡ e46cd1dc-fd48-11ea-0802-4d13a4981a23
begin
	m_sim = zeros(size(post5_3_As_df, 1), length(a_seq))
end;

# ╔═╡ 9899e134-fd46-11ea-0499-b94859cad8d1
for j in 1:size(post5_3_As_df, 1)
  for i in 1:length(a_seq)
    d = Normal(part5_3_As.aM[j] + part5_3_As.bAM[j]*a_seq[i], part5_3_As.sigma_M[j])
    m_sim[j, i] = rand(d, 1)[1]
  end
end

# ╔═╡ 98a1fc3e-fd46-11ea-12f3-81a21baa6353
md"## Snippet 5.27"

# ╔═╡ eee2e318-fd48-11ea-2433-e1f6e65a082a
d_sim = zeros(size(post5_3_As_df, 1), length(a_seq));

# ╔═╡ 98a9de04-fd46-11ea-1a1b-b7512b456dc6
for j in 1:size(post5_3_As_df, 1)
  for i in 1:length(a_seq)
    d = Normal(part5_3_As.a[j] + part5_3_As.bA[j]*a_seq[i] + part5_3_As.bM[j]*m_sim[j, i], part5_3_As.sigma[j])
    d_sim[j, i] = rand(d, 1)[1]
  end
end

# ╔═╡ 98ac248e-fd46-11ea-37ca-fbce7e1a8203
begin
	plot(xlab="Manipulated A", ylab="Counterfactual D",
		title="Total counterfactual effect of A on D")
	plot!(a_seq, mean(d_sim, dims=1)[1, :], leg=false)
	hpdi_array = zeros(length(a_seq), 2)
	for i in 1:length(a_seq)
		hpdi_array[i, :] =  hpdi(d_sim[i, :])
	end
	plot!(a_seq, mean(d_sim, dims=1)[1, :]; ribbon=(hpdi_array[:, 1], -hpdi_array[:, 2]))
end

# ╔═╡ 98d0d7e8-fd46-11ea-02d8-191f0f81da25
md"## End of clip-05-25-27s.jl"

# ╔═╡ Cell order:
# ╟─f7be0558-fd45-11ea-2cb4-c9f411d2e55e
# ╠═98601dc8-fd46-11ea-2560-e762bfd97ed7
# ╠═986050d6-fd46-11ea-26b6-7f618638f1ab
# ╠═986d9d7c-fd46-11ea-305a-915f690d446a
# ╠═946d9a5a-833e-11eb-2c8e-452d18e7c2db
# ╠═946dcbb0-833e-11eb-25c9-af83df7acab1
# ╠═946eb598-833e-11eb-2df5-5f66fa9a18d7
# ╠═948f1e3c-833e-11eb-391d-6b4ec346d46c
# ╟─9877c52c-fd46-11ea-293f-3d07c0cd6734
# ╠═b9af5a04-fd49-11ea-07ba-fb61574aed90
# ╟─988ca05a-fd46-11ea-2ae3-910f7baa2b3f
# ╠═c79668d4-fd48-11ea-0cea-838eb6be744c
# ╠═988d33a8-fd46-11ea-27e5-7ba54e7b04fa
# ╠═e46cd1dc-fd48-11ea-0802-4d13a4981a23
# ╠═9899e134-fd46-11ea-0499-b94859cad8d1
# ╟─98a1fc3e-fd46-11ea-12f3-81a21baa6353
# ╠═eee2e318-fd48-11ea-2433-e1f6e65a082a
# ╠═98a9de04-fd46-11ea-1a1b-b7512b456dc6
# ╠═98ac248e-fd46-11ea-37ca-fbce7e1a8203
# ╟─98d0d7e8-fd46-11ea-02d8-191f0f81da25


text:
# Autogenerated wrapper script for libnabo_jll for i686-w64-mingw32-cxx11
export libnabo

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libnabo")
JLLWrappers.@declare_library_product(libnabo, "libnabo.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libnabo,
        "bin\\libnabo.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Flux: calc_padding, expand, @nograd
import Flux: SamePad, DenseConvDims, PoolDims
import Flux: Conv, DepthwiseConv, CrossCor, MaxPool, MeanPool
import Flux: depthwiseconv, crosscor, maxpool, meanpool

export CausalPad

struct CausalPad end

function calc_padding(lt, pad::CausalPad, k::NTuple{N,T}, dilation, stride) where {N,T}
  @assert length(k) == 1
  # Effective kernel size, including dilation
  k_eff = @. k + (k - 1) * (dilation - 1)
  # How much total padding needs to be applied?
  pad_amt = @. k_eff - 1
  return Tuple(mapfoldl(i -> [i, 0], vcat, pad_amt))
end

need_manual_padding(x, pad) = unwrap(x) isa Flux.CUDA.CuArray && pad[1:2:end] != pad[2:2:end]

@nograd need_manual_padding

function get_paddings(x, pad)
  xl = fill!(similar(x, pad[1:2:end]..., size(x)[end-1:end]...), 0)
  xr = fill!(similar(x, pad[2:2:end]..., size(x)[end-1:end]...), 0)
  return xl, xr
end

@nograd get_paddings

function add_padding(x, pad::NTuple{N, Int}) where N
  all(iszero, pad) && return x
  if sum(pad) > 0
    xl, xr = get_paddings(x, pad)
    cat(xl, x, xr, dims = 1:(N ÷ 2))
  else
    ids = ntuple(Val(N ÷ 2)) do d
      l, r = pad[2d - 1], pad[2d]
      (1 - l):(size(x, d) + r)
    end
    x[ids..., :, :]
  end
end

macro fixpad(ex)
  var = filter(x -> x isa Symbol, ex.args)[2]
  pad = nothing
  for ex′ in ex.args[2].args
    if ex′ isa Expr && ex′.head == :kw && ex′.args[1] == :padding
      pad = ex′.args[2]
      ex′.args[2] = :(needpad ? zero.($pad) : $pad)
      break
    end
  end
  isnothing(pad) && return esc(ex)
  quote
    needpad = need_manual_padding($var, $pad)
    $var = needpad ? add_padding($var, $pad) : $var
    $ex
  end |> esc
end

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

function (c::Conv)(x::AbstractArray)
  # TODO: breaks gpu broadcast :(
  # ndims(x) == ndims(c.weight)-1 && return squeezebatch(c(reshape(x, size(x)..., 1)))
  σ, b = c.σ, reshape(c.bias, ntuple(_->1, length(c.stride))..., :, 1)
  cdims = @fixpad DenseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  σ.(conv(x, c.weight, cdims) .+ b)
end

function (c::DepthwiseConv)(x)
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  cdims = @fixpad DepthwiseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  σ.(depthwiseconv(x, c.weight, cdims) .+ b)
end

function (c::CrossCor)(x::AbstractArray)
  # TODO: breaks gpu broadcast :(
  # ndims(x) == ndims(c.weight)-1 && return squeezebatch(c(reshape(x, size(x)..., 1)))
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  cdims = @fixpad DenseConvDims(x, c.weight; stride=c.stride, padding=c.pad, dilation=c.dilation)
  σ.(crosscor(x, c.weight, cdims) .+ b)
end

function (m::MaxPool)(x)
  pdims = @fixpad PoolDims(x, m.k; padding=m.pad, stride=m.stride)
  return maxpool(x, pdims)
end

function (m::MeanPool)(x)
  pdims = @fixpad PoolDims(x, m.k; padding=m.pad, stride=m.stride)
  return meanpool(x, pdims)
end

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

function (c::ConvTranspose)(x::AbstractArray)
  # ndims(x) == ndims(c.weight)-1 && return squeezebatch(c(reshape(x, size(x)..., 1)))
  σ, b = c.σ, reshape(c.bias, map(_->1, c.stride)..., :, 1)
  if need_manual_padding(x, c.pad)
    c0 = ConvTranspose(c.σ, c.weight, c.bias, c.stride, zero.(c.pad), c.dilation)
    cimds = conv_transpose_dims(c0, x)
    add_padding(σ.(∇conv_data(x, c.weight, cdims) .+ b), .-c.pad)
  else
    cdims = conv_transpose_dims(c, x)
    σ.(∇conv_data(x, c.weight, cdims) .+ b)
  end
end


text:
module ArrayLayouts
using Base, Base.Broadcast, LinearAlgebra, FillArrays
import LinearAlgebra.BLAS

import Base: AbstractArray, AbstractMatrix, AbstractVector,
        ReinterpretArray, ReshapedArray, AbstractCartesianIndex, Slice,
             RangeIndex, BroadcastStyle, copyto!, length, broadcastable, axes,
             getindex, eltype, tail, IndexStyle, IndexLinear, getproperty,
             *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
         fld, cld, div, min, max, minimum, maximum, mod,
         <, ≤, >, ≥, promote_rule, convert, copy,
         size, step, isempty, length, first, last, ndims,
         getindex, setindex!, intersect, @_inline_meta, inv,
         sort, sort!, issorted, sortperm, diff, cumsum, sum, in, broadcast,
         eltype, parent, real, imag,
         conj, transpose, adjoint, permutedims, vec,
         exp, log, sqrt, cos, sin, tan, csc, sec, cot,
                   cosh, sinh, tanh, csch, sech, coth,
                   acos, asin, atan, acsc, asec, acot,
                   acosh, asinh, atanh, acsch, asech, acoth, (:),
         AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo, one, zero,
        to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
      checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
      _in_range, _range, _rangestyle, Ordered,
      ArithmeticWraps, floatrange, reverse, unitrange_last,
      AbstractArray, AbstractVector, axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
      similar, @_gc_preserve_end, @_gc_preserve_begin,
      @nexprs, @ncall, @ntuple, tuple_type_tail,
      all, any, isbitsunion, issubset, replace_in_print_matrix, replace_with_centered_mark

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                        combine_eltypes, DefaultArrayStyle, instantiate, materialize,
                        materialize!, eltypes

import LinearAlgebra: AbstractTriangular, AbstractQ, checksquare, pinv, fill!, tilebufsize, Abuf, Bbuf, Cbuf, dot, factorize, qr, lu, cholesky,
                        norm2, norm1, normInf, normMinusInf, qr, lu, qr!, lu!

import LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex

import FillArrays: AbstractFill, getindex_value

if VERSION < v"1.2-"
    import Base: has_offset_axes
    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
else
    import Base: require_one_based_indexing
end

export materialize, materialize!, MulAdd, muladd!, Ldiv, Rdiv, Lmul, Rmul, lmul, rmul, ldiv, rdiv, mul, MemoryLayout, AbstractStridedLayout,
        DenseColumnMajor, ColumnMajor, ZerosLayout, FillLayout, AbstractColumnMajor, RowMajor, AbstractRowMajor,
        DiagonalLayout, ScalarLayout, SymTridiagonalLayout, HermitianLayout, SymmetricLayout, TriangularLayout,
        UnknownLayout, AbstractBandedLayout, ApplyBroadcastStyle, ConjLayout, AbstractFillLayout,
        colsupport, rowsupport, layout_getindex, QLayout, LayoutArray, LayoutMatrix, LayoutVector

struct ApplyBroadcastStyle <: BroadcastStyle end
@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyBroadcastStyle})
    @assert length(bc.args) == 1
    copyto!(dest, first(bc.args))
end

# Subtypes of LayoutArray default to
# ArrayLayouts routines
abstract type LayoutArray{T,N} <: AbstractArray{T,N} end
const LayoutMatrix{T} = LayoutArray{T,2}
const LayoutVector{T} = LayoutArray{T,1}

include("memorylayout.jl")
include("muladd.jl")
include("lmul.jl")
include("ldiv.jl")
include("diagonal.jl")
include("triangular.jl")
include("factorizations.jl")

@inline sub_materialize(_, V, _) = Array(V)
@inline sub_materialize(L, V) = sub_materialize(L, V, axes(V))
@inline sub_materialize(V::SubArray) = sub_materialize(MemoryLayout(typeof(V)), V)

@inline layout_getindex(A, I...) = sub_materialize(view(A, I...))

macro layoutmatrix(Typ)
    esc(quote
        ArrayLayouts.@layoutldiv $Typ
        ArrayLayouts.@layoutmul $Typ
        ArrayLayouts.@layoutlmul $Typ
        ArrayLayouts.@layoutfactorizations $Typ

        @inline Base.getindex(A::$Typ, kr::Colon, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
        @inline Base.getindex(A::$Typ, kr::Colon, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr)
        @inline Base.getindex(A::$Typ, kr::AbstractUnitRange, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
        @inline Base.getindex(A::$Typ, kr::AbstractUnitRange, jr::AbstractUnitRange) = ArrayLayouts.layout_getindex(A, kr, jr)
        @inline Base.getindex(A::$Typ, kr::AbstractVector, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
        @inline Base.getindex(A::$Typ, kr::Colon, jr::AbstractVector) = ArrayLayouts.layout_getindex(A, kr, jr)
        @inline Base.getindex(A::$Typ, kr::AbstractVector, jr::Colon) = ArrayLayouts.layout_getindex(A, kr, jr)
    end)
end

@layoutmatrix LayoutMatrix

_copyto!(_, _, dest::AbstractArray{T,N}, src::AbstractArray{V,N}) where {T,V,N} =
    Base.invoke(copyto!, Tuple{AbstractArray{T,N},AbstractArray{V,N}}, dest, src)

copyto!(dest::LayoutArray{<:Any,N}, src::LayoutArray{<:Any,N}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)
copyto!(dest::AbstractArray{<:Any,N}, src::LayoutArray{<:Any,N}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)
copyto!(dest::LayoutArray{<:Any,N}, src::AbstractArray{<:Any,N}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)

copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::SubArray{<:Any,N,<:LayoutArray}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)
copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::LayoutArray{<:Any,N}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)
copyto!(dest::LayoutArray{<:Any,N}, src::SubArray{<:Any,N,<:LayoutArray}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)
copyto!(dest::SubArray{<:Any,N,<:LayoutArray}, src::AbstractArray{<:Any,N}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)
copyto!(dest::AbstractArray{<:Any,N}, src::SubArray{<:Any,N,<:LayoutArray}) where N =
    _copyto!(MemoryLayout(typeof(dest)), MemoryLayout(typeof(src)), dest, src)

zero!(A::AbstractArray{T}) where T = fill!(A,zero(T))
function zero!(A::AbstractArray{<:AbstractArray})
    for a in A
        zero!(a)
    end
    A
end

_fill_lmul!(β, A::AbstractArray{T}) where T = iszero(β) ? zero!(A) : lmul!(β, A)

# Elementary reflection similar to LAPACK. The reflector is not Hermitian but
# ensures that tridiagonalization of Hermitian matrices become real. See lawn72
@inline function reflector!(x::AbstractVector)
    require_one_based_indexing(x)
    n = length(x)
    n == 0 && return zero(eltype(x))
    @inbounds begin
        ξ1 = x[1]
        normu = abs2(ξ1)
        for i = 2:n
            normu += abs2(x[i])
        end
        if iszero(normu)
            return zero(ξ1/normu)
        end
        normu = sqrt(normu)
        ν = copysign(normu, real(ξ1))
        ξ1 += ν
        x[1] = -ν
        for i = 2:n
            x[i] /= ξ1
        end
    end
    ξ1/ν
end

# apply reflector from left
@inline function reflectorApply!(x::AbstractVector, τ::Number, A::AbstractVecOrMat)
    m,n = size(A,1),size(A,2)
    if length(x) != m
        throw(DimensionMismatch("reflector has length $(length(x)), which must match the first dimension of matrix A, $m"))
    end
    m == 0 && return A
    @inbounds begin
        for j = 1:n
            # dot
            vAj = A[1, j]
            for i = 2:m
                vAj += x[i]'*A[i, j]
            end

            vAj = conj(τ)*vAj

            # ger
            A[1, j] -= vAj
            for i = 2:m
                A[i, j] -= x[i]*vAj
            end
        end
    end
    return A
end

end


text:
module H3

include("../gen/Lib.jl")
include("API.jl")

end # module H3