text:
export rres, rresx

################################################################################
#
#  Computation of a generator of the ideal of the resultant
#
################################################################################
@doc Markdown.doc"""
    resultant_ideal(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
A generator for the ideal of the resultant of $f$ anf $g$ using a quadratic-time algorithm.
One of the two polynomials must be monic.
"""
function resultant_ideal(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  #The algorithm is the same as the resultant. We assume that one fo the 2 polynomials is monic. Under this assumption, at every
  #step the same is true and we can discard the unti obtained from the fun_factor function
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  easy = isprime_power(m)

  if easy
    return resultant_ideal_pp(f,g)
  end

  #Some initial checks
  res = R(1)
  if degree(f) < 1 && degree(g) < 1
    if iszero(f) || iszero(g)
      return R(0)
    end
    return res
  end

  if degree(f) < 1
    res = mul!(res, res, lead(f)^degree(g))
    return res
  end

  c, f = primsplit(f)
  if !isone(c)
    res = mul!(res, res, R(c)^degree(g))
  end

  c, g = primsplit(g)
  if !isone(c)
    res = mul!(res, res, R(c)^degree(f))
  end

  if degree(f) < degree(g)
    f, g = g, f
  end

  if iszero(res)
    return res
  end

  #Now, I can safely assume that the degree of f is always greater than the degree of g
  while true

    if degree(g) < 1
      res = mul!(res, res, lead(g)^degree(f))
      return res
    end

    c, g = primsplit(g)
    if !isone(c)
      res = mul!(res, res, R(c)^degree(f))
    end

    if iszero(res)
      return res
    end
    #want f % g which works iff lead(g) | lead(f)

    if isunit(lead(g)) #accelerate a bit...possibly.
      rem!(f, f, g)
      f, g = g, f
      continue
    end
    break
  end

  # factoring case, need to split the ring as well.
  # we need a coprime factorisation and then we go recursively
  #If res is not coprime to the modulus, I can continue the computations modulo a smaller one.
  s = gcd(m, lift(res))
  if !isone(s)
    m = divexact(m, s)
  end
  cp = fmpz[gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
  push!(cp, m)
  cp = fmpz[x for x = cp if !iszero(x)]
  cp = coprime_base(cp)
  cp = fmpz[x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
  resp = fmpz[]
  pg = fmpz[]
  for p = cp
    lg = p^valuation(m, p)
    push!(pg, lg)
    R1 = ResidueRing(FlintZZ, S(lg), cached = false)
    R1t = PolynomialRing(R1, cached = false)[1]
    #g is bad in R1, so factor it
    gR1 = R1t(T[R1(lift(coeff(g, i))) for i = 0:degree(g)])
    fR1 = R1t(T[R1(lift(coeff(f, i))) for i = 0:degree(f)])

    if degree(fR1) < degree(f) && degree(gR1) < degree(g)
      res1 = R1(0)
    elseif degree(fR1) < degree(f)
      res1 = lead(gR1)^(degree(f) - degree(fR1))
    else
      res1 = R1(1)
    end

    if !isunit(lead(gR1))
      g1, g2 = fun_factor(gR1)
      res1 = mul!(res1, res1, resultant_ideal(fR1, g2))
      push!(resp, lift(res1))
    else
      #gR1 has a invertible leading coeff
      res1 = mul!(res1, res1, resultant_ideal(fR1, gR1))
      push!(resp, lift(res1))
    end
  end
  if length(resp) == 1
    res = mul!(res, res, R(resp[1]))
  else
    res = mul!(res, res, R(crt(resp, pg)))
  end
  return res
end

function resultant_ideal_pp(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  #The algorithm is the same as the resultant. We assume that one fo the 2 polynomials is monic. Under this assumption, at every
  #step the same is true and we can discard the unti obtained from the fun_factor function
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  #pn = fmpz(modulus(R))
  pn = modulus(R)

  #Some initial checks
  if degree(f) < 1 && degree(g) < 1
    if iszero(f) || iszero(g)
      return R(0)
    else
      return R(1)
    end
  end

  res = R(1)

  if degree(f) < 1
    res = mul!(res, res, lead(f)^degree(g))
    return res
  end

  if degree(g) < 1
    res = mul!(res, res, lead(g)^degree(f))
    return res
  end

  c, f = primsplit(f)
  if !isone(c)
    res = mul!(res, res, c^degree(g))
  end

  c, g = primsplit(g)
  if !isone(c)
    res = mul!(res, res, c^degree(f))
  end

  if degree(f) < degree(g)
    f, g = g, f
  end

  if iszero(res)
    return res
  end

  while true
    #want f % g which works iff lead(g) | lead(f)

    if isunit(lead(g)) #accelerate a bit...possibly.
      rem!(f, f, g)
      if degree(f) < 1
        res = mul!(res, res, lead(f)^degree(g))
        return res
      end
      c, f = primsplit!(f)
      if !isone(c)
        res = mul!(res, res, c^degree(g))
      end
      f, g = g, f
    else
      s = gcd(lift(res), pn)
      if !isone(s)
        new_pn = divexact(pn, s)
        R1 = ResidueRing(FlintZZ, S(new_pn), cached = false)
        R1t = PolynomialRing(R1, "y", cached = false)[1]
        f2 = R1t(T[R1(lift(coeff(f, i))) for i = 0:degree(f)])
        g2 = R1t(T[R1(lift(coeff(g, i))) for i = 0:degree(g)])
        g2 = fun_factor(g2)[2]
        z = resultant_ideal_pp(f2, g2)
        return mul!(res, res, R(lift(z)))
      end
      g = fun_factor(g)[2]
      if degree(g) < 1
        res = mul!(res, res, lead(g)^degree(f))
        return res
      end
    end
  end

end

################################################################################
#
#  Extended gcd
#
################################################################################

@doc Markdown.doc"""
    xxgcd(a::ResElem{fmpz}, b::ResElem{fmpz}) -> g, e, f, u, v
    xxgcd(a::ResElem{Integer}, b::ResElem{Integer}) -> g, e, f, u, v
Coputes $g = \gcd(a, b)$, the Bezout coefficients, $e$, $f$ s.th.
$g = ea+fb$ and $u$, $v$ s.th. $ev-fu = 1$, $gu = a$ and $gv = b$
"""
function xxgcd(a::ResElem{S}, b::ResElem{S}) where S <: Union{fmpz, Integer}
  g, e, f = gcdx(a, b)
  if iszero(g)
    return g, e, f, f, e
  else
    return g, e, f, divexact(a, g), divexact(b, g)
  end
end

#for testing: the obvious(?) naive method(s)
function rres_hnf(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  R = base_ring(f)
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  s = Nemo.Generic.sylvester_matrix(lift(Zx, f), lift(Zx, g))
  h = hnf(s)
  return gcd(R(0), R(h[nrows(h), ncols(h)]))
end

function rres_bez(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  R = base_ring(f)
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  Qx = PolynomialRing(FlintQQ, "x", cached = false)[1]
  g, q, w = gcdx(Qx(lift(Zx, f)), Qx(lift(Zx, g)))
  return gcd(R(0), R(lcm(denominator(q), denominator(w))))
end

#polynomial remainder sequence - almost
#=
function prs_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}

  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)
  @assert easy

  Zx = PolynomialRing(FlintZZ, cached = false)[1]

  rs = []

  while degree(g) >0
    g *= inv(canonical_unit(lead(g)))
    c, gg = primsplit(g)
    @show f, (g, mu) = gg, my_divrem(f, gg)
    push!(rs, (c, gg, mu))
  end
  return rs, g
end
=#

function Nemo.gcd(a::ResElem{T}, b::ResElem{T}) where T <: Union{Integer, fmpz}
  m = modulus(a)
  return parent(a)(gcd(gcd(a.data, m), b.data))
end

function Nemo.gcdx(a::ResElem{T}, b::ResElem{T}) where T <: Union{Integer, fmpz}
  m = modulus(a)
  R = parent(a)
  g, u, v = gcdx(fmpz(a.data), fmpz(b.data))
  G, U, V = gcdx(g, fmpz(m))
  return R(G), R(U)*R(u), R(U)*R(v)
end

@doc Markdown.doc"""
    annihilator(a::ResElem{fmpz}) -> r
    annihilator(a::ResElem{Integer}) -> r
The annihilator of $a$, ie. a generator for the annihilator ideal
$\{x | xa = 0\} = \langle r\rangle$
"""
function annihilator(a::ResElem{T}) where T <: Union{Integer, fmpz}
  R = parent(a)
  m = modulus(R)
  return R(divexact(m, gcd(m, a.data)))
end

@doc Markdown.doc"""
    isunit(f::Union{fmpz_mod_poly,nmod_poly}) -> Bool
Tests if $f$ is a unit in the polynomial ring, ie. if
$f = u + n$ where $u$ is a unit in the coeff. ring
and $n$ is nilpotent.
"""
function Nemo.isunit(f::T) where T <: Union{fmpz_mod_poly,nmod_poly}
  if !isunit(constant_coefficient(f))
    return false
  end
  for i=1:degree(f)
    if !isnilpotent(coeff(f, i))
      return false
    end
  end
  return true
end

@doc Markdown.doc"""
    isnilpotent(a::ResElem{fmpz}) -> Bool
    isnilpotent(a::ResElem{Integer}) -> Bool
Tests if $a$ is nilpotent.
"""
function isnilpotent(a::ResElem{T}) where T <: Union{Integer, fmpz}
  #a is nilpontent if it is divisible by all primes divising the modulus
  # the largest exponent a prime can divide is nbits(m)
  l = nbits(modulus(a))
  return iszero(a^l)
end

function iszerodivisor(f::T) where T <: Union{fmpz_mod_poly,nmod_poly}
  c = content(f)
  return isnilpotent(c)
end

function Nemo.inv(f::T) where T <: Union{fmpz_mod_poly,nmod_poly}
  if !isunit(f)
    error("impossible inverse")
  end
  Rx = parent(f)
  g = Rx(inv(trailing_coefficient(f)))
  #lifting: to invert a, start with an inverse b mod m, then
  # then b -> b*(2-ab) is an inverse mod m^2
  # starting with this g, and using the fact that all coeffs are nilpotent
  # we have an invers modulo s.th. nilpotent. Hence it works
  c = Rx()
  mul!(c, f, g)
  while !isone(c)
    mul!(g, g, 2-c)
    mul!(c, f, g)
  end
  return g
end

function Nemo.invmod(f::fmpz_mod_poly, M::fmpz_mod_poly)
  if !isunit(f)
    error("impossible inverse")
  end
  if !isunit(lead(M))
    error("not yet implemented")
  end
  g = parent(f)(inv(trailing_coefficient(f)))
  #lifting: to invert a, start with an inverse b mod m, then
  # then b -> b*(2-ab) is an inverse mod m^2
  # starting with this g, and using the fact that all coeffs are nilpotent
  # we have an invers modulo s.th. nilpotent. Hence it works
  c = f*g
  rem!(c, c, M)
  while !isone(c)
    mul!(g, g, 2-c)
    rem!(g, g, M)
    mul!(c, f, g)
    rem!(c, c, M)
  end
  return g
end

function rres_sircana_pp(f1::PolyElem{T}, g1::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f1, g1)
  @assert typeof(f1) == typeof(g1)
  Rt = parent(f1)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  f = deepcopy(f1)
  g = deepcopy(g1)

  res = R(1)
  while true
    #First, some trivial check.
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) && iszero(g)
        res = zero(R)
      else
        res = mul!(res, res, R(gcd(lift(lead(f)), lift(lead(g)))))
      end
      return res
    end

    if degree(f) < degree(g)
      f, g = g, f
    end
    if degree(g) < 1

      if divisible(lift(coeff(f, degree(f))), p)
        #need the constant coeff * the annihilator of the others...
        a = coeff(f, 1)
        for i = 2:degree(f)
          a = gcd(a, coeff(f, i))
          if isone(a)
            break
          end
        end
        a = annihilator(a)
        if iszero(a)
          return res*coeff(g, 0)
        else
          res = mul!(res, res, R(gcd(lift(coeff(g, 0)), lift(a)*lift(coeff(f, 0)))))
          return res
        end
      else
        res = mul!(res, res, constant_coefficient(g))
        return res
      end
    end

    c, g = primsplit!(g)
    res = mul!(res, res, c)
    if iszero(res)
      return res
    end

    if divisible(lift(coeff(g, degree(g))), p)
      #one of the coefficient will now be invertible (at least after the splitting)
      g = fun_factor(g)[2]
    end
    rem!(f, f, g)
  end
end

# Based on these formulas:
# rres(f, g) = rres(f - kg, g), so I can divide.
# rres(f, g) = rres(g, f)
# rres(uf, g) = rres(f, g)
# rres(pf, g) mod p^n = p*(rres(f, g) mod p^(n-1)) (under right hypotheses)
function rres_sircana(f1::PolyElem{T}, g1::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f1, g1)
  @assert typeof(f1) == typeof(g1)
  Rt = parent(f1)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  easy = isprime_power(m)
  if easy
    return rres_sircana_pp(f1, g1)
  end
  f = deepcopy(f1)
  g = deepcopy(g1)

  Zx = PolynomialRing(FlintZZ, cached = false)[1]

  res = R(1)
  while true
    #First, some trivial check.
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) && iszero(g)
        res = zero(R)
      else
        res *= R(gcd(lift(lead(f)), lift(lead(g))))
      end
      return res
    end

    if degree(f) < degree(g)
      f, g = g, f
    end

    if degree(g) < 1
      if !isunit(lead(f))
        #need the constant coeff * the annihilator of the others...
        a = coeff(f, 1)
        for i = 2:degree(f)
          a = gcd(a, coeff(f, i))
          if isone(a)
            break
          end
        end
        a = annihilator(a)
        if iszero(a)
          return lead(g)
        else
          res *= gcd(lead(g), a*constant_coefficient(f))
          return res
        end
      else
        res *= constant_coefficient(g)
        return res
      end
    end

    c, g = primsplit!(g)
    res *= R(c)
    if iszero(res)
      return res
    end

    if !isunit(lead(g))
      #one of the coefficient will now be invertible (at least after the splitting)
      s = gcd(m, lift(res))
      if !isone(s)
        m = divexact(m, s)
      end
      if easy
        cp = S[m]
      else
        cp = S[gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
        push!(cp, m)
        cp = S[x for x = cp if !iszero(x)]
        cp = coprime_base(cp)
        cp = S[x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
      end
      resp = fmpz[]
      pg = fmpz[]
      for p = cp
        lg = p^valuation(m, p)
        push!(pg, lg)
        R1 = ResidueRing(FlintZZ, S(lg), cached=false)
        R1t = PolynomialRing(R1, cached=false)[1]
        #g is bad in R1, so factor it
        gR1 = R1t(lift(Zx, g))
        fR1 = R1t(lift(Zx, f))
        res1 = one(R1)
        if isunit(lead(gR1))
          g2 = gR1
        else
          if iszero(gR1)
            push!(resp, fmpz(0))
            continue
          end
          c, gR1 = primsplit!(gR1)
          res1 *= R1(c)
          g1, g2 = fun_factor(gR1)
        end
        push!(resp, lift(res1)*lift(rres_sircana(fR1, g2)))
      end
      if length(resp) == 1
        res *= R(resp[1])
      else
        res *= R(crt(resp, pg))
      end
      return res
    end

    f = rem(f, g)
  end
end

@doc Markdown.doc"""
    rresx(f::PolyElem{ResElem{fmpz}}, g::PolyElem{ResElem{fmpz}}) -> r, u, v
The reduced resultant $r$ and polynomials $u$ and $v$ s.th.
$r = uf+vg$ and $\langle r\rangle = \langle f, g\rangle\cap \mathbb Z$.
"""
function rresx(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  return rresx_sircana(f, g)
end

function rresx_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  @assert isunit(lead(f)) || isunit(lead(g)) #can be weakened to invertable lead
  res, u, v = _rresx_sircana(f, g)
  if !iszero(res)
    cu = canonical_unit(res)
    cu = inv(cu)
    res = mul!(res, res, cu)
    u *= cu
    v *= cu
  end
  if isunit(lead(g))
    q, r = divrem(u, g)
    @hassert :NfOrd 1 res == r*f + (v+q*f)*g
    return res, r, v+q*f
  else
    q, r = divrem(v, f)
    @hassert :NfOrd 1 res == (u+q*g)*f + r*g
    return res, u+q*g, r
  end
end

function rresx_sircana_pp(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  @assert isunit(lead(f)) || isunit(lead(g)) #can be weakened to invertable lead
  res, u, v = _rresx_sircana_pp(f, g)
  if !iszero(res)
    cu = canonical_unit(res)
    cu = inv(cu)
    res = mul!(res, res, cu)
    u *= cu
    v *= cu
  end
  if isunit(lead(g))
    q, r = divrem(u, g)
    @hassert :NfOrd 1 res == r*f + (v+q*f)*g
    return res, r, v+q*f
  else
    q, r = divrem(v, f)
    @hassert :NfOrd 1 res == (u+q*g)*f + r*g
    return res, u+q*g, r
  end
end

function _rresx_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)

  Rt = parent(f)
  R = base_ring(Rt)
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  m = fmpz(modulus(R))
  easy = isprime_power(m)
  if easy
    return _rresx_sircana_pp(f, g)
  end

  u, v = Rt(0), Rt(1)
  U, V = Rt(1), Rt(0)

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        res = R(0)
        u = Rt(0)
        v = Rt(0)
      else
        res, uu, vv = gcdx(lead(f), lead(g))
        #res = uu*f + vv*g = uu*(U f_in + V g_in) + vv*(u f_in + v g_in)
        #    = uu*U + vv*u  uu*V + vv*v
        u, v = (uu*U + vv*u), (uu*V + vv*v)
      end
      return res, u, v
    end

    if degree(f) < degree(g)
      f, g = g, f
      U, u = u, U
      V, v = v, V
    end

    if degree(g) < 1
      if !isunit(lead(f))
        #need the constant coeff * the annihilator of the others...
        a = coeff(f, 1)
        for i = 2:degree(f)
          a = gcd(a, coeff(f, i))
          if isone(a)
            break
          end
        end
        a = annihilator(a)
        if iszero(a)
          return lead(g), u, v
        else
          res, uu, vv = gcdx(a*constant_coefficient(f), lead(g))
          return res, uu*a*U + vv*u, uu*a*V + vv*v
        end
      else
        return constant_coefficient(g), u, v
      end
    end

    if !isunit(lead(g))
      c, g = primsplit(g)
      cp = fmpz[gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
      push!(cp, m)
      cp = coprime_base(cp)
      cp = fmpz[x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
      resp = fmpz[]
      resB = Tuple{typeof(f), typeof(f)}[]
      pg = fmpz[]
      for p = cp
        lg = p^valuation(m, p)
        push!(pg, lg)
        R1 = ResidueRing(FlintZZ, S(lg), cached=false)
        R1t = PolynomialRing(R1, cached=false)[1]
        #g is bad in R1, so factor it
        gR1 = R1t(lift(Zx, g))
        fR1 = R1t(lift(Zx, f))

        if iszero(R1(lift(c)))
          push!(resp, fmpz(0))
          push!(resB, (R1t(0), R1t(1))) #relation need to be primitive
        else
          if isunit(lead(gR1))
            g2 = gR1
            g1 = R1t(1)
          else
            g1, g2 = fun_factor(gR1)
          end
          #@assert isunit(g1)
          rr, uu, vv = rresx_sircana(fR1, g2)
          push!(resp, lift(lift(c)*rr))
          push!(resB, (uu*lift(c), inv(g1)*vv))
        end
      end
      if length(cp) == 1
        res, u_, v_ = R(resp[1]), Rt(lift(Zx, resB[1][1])), Rt(lift(Zx, resB[1][2]))
      else
        ce = crt_env(pg)
        res = R(crt(resp, ce))
        u_ = Rt(induce_crt(typeof(f)[x[1] for x = resB], ce))
        v_ = Rt(induce_crt(typeof(f)[x[2] for x = resB], ce))
      end
      # f = U*f_in + V*g_in
      # g = u*f_in + v*g_in
      # r = u_ * f + v_ * g
      return res, (u_*U + v_*u), (u_*V + v_*v)
    end

    q, f = divrem(f, g)
    #f -> f-qg, so U*f_in + V * g_in -> ...
    U = U - q*u
    V = V - q*v
  end
end

function _rresx_sircana_pp(f1::PolyElem{T}, g1::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f1, g1)

  Rt = parent(f1)
  R = base_ring(Rt)
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  m = fmpz(modulus(R))
  f = f1
  g = g1

  u, v = Rt(0), Rt(1)
  U, V = Rt(1), Rt(0)

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) && iszero(g)
        res = R(0)
        u = Rt(0)
        v = Rt(0)
      else
        res, uu, vv = gcdx(lead(f), lead(g))
        #res = uu*f + vv*g = uu*(U f_in + V g_in) + vv*(u f_in + v g_in)
        #    = uu*U + vv*u  uu*V + vv*v
        u, v = (uu*U + vv*u), (uu*V + vv*v)
      end
      return res, u, v
    end

    if degree(f) < degree(g)
      f, g = g, f
      U, u = u, U
      V, v = v, V
    end

    if degree(g) < 1
      if !isunit(lead(f))
        #need the constant coeff * the annihilator of the others...
        a = coeff(f, 1)
        for i = 2:degree(f)
          a = gcd(a, coeff(f, i))
          if isone(a)
            break
          end
        end
        a = annihilator(a)
        if iszero(a)
          return lead(g), u, v
        else
          res, uu, vv = gcdx(a*constant_coefficient(f), lead(g))
          u, v = (uu*a*U + vv*u), (uu*a*V + vv*v)

          return res, u, v
        end
      else
        return constant_coefficient(g), u, v
      end
    end

    if !isunit(lead(g))
      c, g = primsplit(g)
      g1, g2 = fun_factor(g)
      rr, uu, vv = _rresx_sircana_pp(f, g2)
      res = rr*c
      u_ = uu*c
      v_ = inv(g1)*vv
      return res, (u_*U + v_*u), (u_*V + v_*v)
    end

    q, f = divrem(f, g)
    #f -> f-qg, so U*f_in + V * g_in -> ...
    U = U - q*u
    V = V - q*v
  end
end

function my_divrem(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  g1, g2 = fun_factor(g)
  if degree(g2) < 1 # g2 is monic, so it's 1
    return parent(f)(0), g2
  end
  u = invmod(g1, g2)
  q, r = divrem(f, g2)
  q, r = divrem(r*u, g2)
  return g1*r, g2
end

#key idea (Carlo): if g = ab and a is a unit mod p, then it is actually a unit
# in Z/p^kZ, hence the ideal (f, g) = (f, b) where b is now monic.
#Thus rres(f,g ) = rres(f, b).... and the division can continue
@doc Markdown.doc"""
    rres(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
The reduced resultant of $f$ and $g$ using a quadratic-time algorithm.
That is a generator for the $(f, g) \cap Z$
"""
function rres(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  return rres_sircana(f, g)
end

################################################################################
#
#  Gcd
#
################################################################################

@doc Markdown.doc"""
    gcd_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
The 'gcd' of $f$ anf $g$ using a quadratic-time algorithm.
"""
function gcd_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)
  @assert easy #for now...

  while !iszero(g)
    @show cg = content(g)
    if !isunit(cg)
      for i=0:degree(g)
        setcoeff!(g, i, divexact(coeff(g, i), cg))
      end
    end
    if !isunit(lead(g))
      u, g = fun_factor(g)
    end
    @show f, g = g, (f%g)
  end
  c = canonical_unit(lead(f))
  f = divexact(f, c)
  return f
end

################################################################################
#
#  Resultant
#
################################################################################

@doc Markdown.doc"""
    resultant_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
The resultant of $f$ anf $g$ using a quadratic-time algorithm.
"""
function resultant_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)

  Zx = PolynomialRing(FlintZZ, cached = false)[1]

  res = R(1)

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        res *= R(0)
      else
        res *= R(1)
      end
      return res
    end

    if degree(f) < 1
      res *= lead(f)^degree(g)
      return res
    end

    if degree(g) < 1
      res *= lead(g)^degree(f)
      return res
    end

    c, f = primsplit(f)
    if !isone(c)
      res *= c^degree(g)
    end

    c, g = primsplit(g)
    if !isone(c)
      res *= c^degree(f)
    end

    if degree(f) < degree(g)
      res *= (-1)^(degree(f)*degree(g))
      f, g = g, f
    end

    #want f % g which works iff lead(g) | lead(f)

    if isunit(lead(g)) #accelerate a bit...possibly.
      q, r = divrem(f, g)
      res *= lead(g)^(degree(f) - degree(r))
      res *= R(-1)^(degree(g)*(degree(f) - degree(r)))
      f = r
      continue
    end

    if gcd(lift(lead(f)), m)  % gcd(lift(lead(g)), m) == 0
      q = divexact(lead(f), lead(g))
      ff = f - q*g*gen(Rt)^(degree(f) - degree(g))
      @assert degree(f) > degree(ff)
      res *= lead(g)^(degree(f) - degree(ff))
      res *= R(-1)^(degree(g)*(degree(f) - degree(ff)))
      f = ff
      continue
    end
    break
  end

  #factoring case, need to split the ring as well.
  #merde : we need a coprime factorisation: take
  # 6t^2+2t+3 mod 36....
  if easy
    cp = [m]
  else
    cp = [gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
    push!(cp, m)
    cp = [x for x = cp if !iszero(x)]
    cp = coprime_base(cp)
    cp = [x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
  end
  resp = fmpz[]
  pg = fmpz[]
  for p = cp
    lg = p^valuation(m, p)
    push!(pg, lg)

    if lg != m
      R1 = ResidueRing(FlintZZ, S(lg), cached=false)
      R1t = PolynomialRing(R1, cached=false)[1]
      #g is bad in R1, so factor it
      gR1 = R1t(lift(Zx, g))
      fR1 = R1t(lift(Zx, f))
    else
      gR1 = g
      fR1 = f
      R1 = R
      R1t = Rt
    end

    if degree(fR1) < degree(f) &&
       degree(gR1) < degree(g)
       res1 = R1(0)
    elseif degree(fR1) < degree(f)
        res1 = R1(-1)^(degree(g) * (degree(f) - degree(fR1))) *
               lead(gR1)^(degree(f) - degree(fR1))
    elseif degree(gR1) < degree(g)
        res1 = lead(fR1)^(degree(g) - degree(gR1))
    else
        res1 = R1(1)
    end

    if !isunit(lead(gR1))
      g1, g2 = fun_factor(gR1)

      #careful: lead(g) = 0 mod lg is possible, so the degree drops!
      #from Wiki:
      # phi: R -> S
      # phi(res(f, g)) = res(phi(f), phi(g)) if the degrees are the same
      #                = 0                   if both degrees drop (zero column in
      #                                         sylvester matrix)
      #                = (-1)^(delta(f)*deg(g))*
      #                  phi(l(g))^delta(f)  if only degree(f) drops
      #                = phi(l(f))^delta(g)  if only degree(g) drops
      #
      # we use res(f, gh) = res(f, g)res(f, h) which seems to be true in general
      #next: res(fR1, gR1) = res(phi(f), g1)
      #g1 is = 1 mod p
      #however, reverse(fR1) can have a smaller degree then fR1 (eg. x^2 -> 1)
      # res(f, g) = t(g)^delta(f) * (-1)^(deg(g)*delta(f)) * res(rev(f), rev(g))
      # there is a (-1)^deg(f)*deg(g) from Wiki for the reverse on top.

      res1 *= resultant_sircana(reverse(fR1), reverse(g1))
      v = valuation(fR1, gen(parent(fR1)))
           #should be delta(f) = degree(rev(f)) - degree(f)
      res1 *= constant_coefficient(g1)^v*R1(-1)^(v*degree(g1))
      res1 *= R1(-1)^(degree(g1)*degree(fR1))

      res1 *= resultant_sircana(fR1, g2)
      push!(resp, lift(res1))
    else
      #gR1 has a invertible leading coeff
      res1 *= resultant_sircana(fR1, gR1)
      push!(resp, lift(res1))
    end
  end
  res *= length(cp)==1 ? R(resp[1]) : R(crt(resp, pg))
  return res
end

################################################################################
#
#  Fun Factor
#
################################################################################
# factors f as unit * monic
# requires some coefficient of f to be a unit

function fun_factor(f::T) where T <: Union{fmpz_mod_poly, nmod_poly}
  Rx = parent(f)
  if isunit(lead(f))
    l = lead(f)
    return Rx(l), f*inv(l)
  end
  R = base_ring(Rx)
  smod = lift(coeff(f, degree(f)))
  ind = degree(f)-1
  while !isunit(coeff(f, ind))
    smod = gcd(smod, lift(coeff(f, ind)))
    ind -= 1
  end
  if iszero(ind)
    return f, one(Rx)
  end

  u0 = zero(Rx)
  for i = 0:degree(f)-ind
    setcoeff!(u0, i, coeff(f, ind+i))
  end
  g0 = zero(Rx)
  invc = inv(coeff(f, ind))
  for i = 0:ind-1
    setcoeff!(g0, i, coeff(f, i) * invc)
  end
  setcoeff!(g0, ind, one(R))

  Zy, y = PolynomialRing(FlintZZ, "y", cached = false)
  f2 = lift(Zy, f)
  mod = fmpz(gcd(smod, fmpz(modulus(Rx)))) #We have the equality modulo mod
  mod = mod*mod
  R1 = ResidueRing(FlintZZ, mod, cached = false)
  R1x, x = PolynomialRing(R1, "x", cached = false)
  s = R1x(lift(inv(coeff(u0, 0))))
  t = zero(R1x)
  u = R1x(lift(Zy, u0))
  g = R1x(lift(Zy, g0))

  f1 = R1x(f2)
  u, g, s, t = _hensel(f1, u, g, s, t)
  @hassert :NfOrd 1 f1 == u*g
  i = 1
  modRx = fmpz(modulus(Rx))
  while modRx > mod

    mod = mod*mod
    if mod > modRx
      mod = modRx
    end
    R1 = ResidueRing(FlintZZ, mod, cached = false)
    R1x, x = PolynomialRing(R1, "x", cached = false)
    u = R1x(lift(Zy, u))
    g = R1x(lift(Zy, g))
    s = R1x(lift(Zy, s))
    t = R1x(lift(Zy, t))
    f1 = R1x(f2)
    i += 1

    u, g, s, t = _hensel(f1, u, g, s, t)

    @hassert :NfOrd 1 f1 == u*g

    if i>20 #in general: loop forever... one could check that the
      # contents of f-gh and s*g+t*h - 1 is nilpotent.
      # this SHOULD ensure convergence
      @show f
      @show parent(f)
      @show g, u, s, t
      @show content(f-g*u)
      @show content(g*t-s*u-1)
      error("too long")
    end
  end

  u0 = Rx(lift(Zy, u))
  g0 = Rx(lift(Zy, g))
  @hassert :NfOrd 1 g0*u0 == f
  return u0, g0
end

function _hensel(f::PolyElem{T}, g::PolyElem{T}, h::PolyElem{T}, s::PolyElem{T}, t::PolyElem{T}) where T <: RingElem #from von zur Gathen: h needs to be monic
  @assert ismonic(h)
  #@assert isnilpotent(content(f-g*h))  #this guarantees useful lifting
  #@assert isnilpotent(content(g*s+h*t-1))
  Rx = parent(f)
  aux1 = Rx()
  mul!(aux1, g, h)
  sub!(aux1, f, aux1)
  #aux1 = f-g*h
  aux2 = Rx()
  mul!(aux2, s, aux1)
  q, r = divrem(aux2, h)
  #@assert s*e == q*h+r
  mul!(aux1, aux1, t)
  add!(aux1, aux1, g)
  mul!(aux2, q, g)
  add!(aux1, aux1, aux2)
  #g = g+t*e+q*g
  g = aux1
  add!(aux2, h, r)
  h = aux2
  #h = h+r
  #@assert ismonic(h)
  aux3 = Rx()
  aux4 = Rx()
  mul!(aux3, s, aux1)
  mul!(aux4, t, aux2)
  add!(aux3, aux3, aux4)
  sub!(aux3, aux3, one(Rx))
  mul!(aux4, s, aux3)
  c, d = divrem(aux4, h)
  #@assert s*b == c*h+d
  mul!(aux3, aux3, t)
  sub!(aux3, t, aux3)
  mul!(c, c, g)
  sub!(aux3, aux3, c)
  sub!(aux4, s, d)
  return aux1, aux2, aux4, aux3
end

################################################################################
#
#  Primitive splitting
#
################################################################################

@doc Markdown.doc"""
    primsplit!(f::PolyElem{ResElem{fmpz}}) -> c, f
    primsplit!(f::PolyElem{ResElem{Integer}}) -> c, f
Computes the contents $c$ and the primitive part of $f$ destructively.
"""
function primsplit!(f::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}

  @assert !iszero(f)
  d = degree(f)
  if d == 0
    if iszero(f)
      return base_ring(parent(f))(1), f
    end
    c = canonical_unit(coeff(f, 0))
    c1 =inv(c)*coeff(f, 0)
    setcoeff!(f, 0, 1)
    return c1, f
  end

  g = coeff(f, 0)
  setcoeff!(f, 0, 1)
  for i = 1:d
    h, _, _, u, v = xxgcd(g, coeff(f, i))
    setcoeff!(f, i, v)
    if  g != h
      for j=0:i-1
        setcoeff!(f, j, u*coeff(f, j))
      end
    end
    g = h
    if isone(g)
      return g, f
    end
  end
  return g, f
end

@doc Markdown.doc"""
    primsplit(f::PolyElem{ResElem{fmpz}}}) -> c, f
    primsplit(f::PolyElem{ResElem{Integer}}}) -> c, f
Computes the contents $c$ and the primitive part of $f$.
"""
function primsplit(f::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  g = deepcopy(f)
  return primsplit!(g)
end

function quo(R::FqNmodPolyRing, f::fq_nmod_poly)
  Q = ResidueRing(R, f)
  mQ = MapFromFunc(x -> Q(x), y -> lift(y), R, Q)
  return Q, mQ
end

#= not finished
function unit_group(R::Generic.ResRing{fq_nmod_poly})
  f = modulus(R)
  lf = factor(f)
  lu = [unit_group_pp(p, k) for (p,k) = f.fac]
  return lu
end

function unit_group_pp(f::fq_nmod_poly, k::Int)
  @assert isirreducible(f)
  k, o = GF(characteristic(parent(f)), degree(base_ring(f))*degree(f))
  #o is a primitive element as per Bill's semantic...
  #however, I need it written wrt. "my" poly
end

=#
function basis(K::FqNmodFiniteField)
  b = fq_nmod[]
  for i=1:degree(K)
    x = K()
    setcoeff!(x, i-1, UInt(1))
    push!(b, x)
  end
  return b
end

function unit_group_1_part(f::fq_nmod_poly, k::Int)
  pr = [k]
  while k>1
    k = div(k+1, 2)
    push!(pr, k)
  end

  K = base_ring(f)
  x = gen(parent(f))
  p = characteristic(K)

  b = basis(K)
  gens = [1+f*x^i*c for i=0:degree(f)-1 for c = b]
  rels = identity_matrix(FlintZZ, length(gens))*p

  for i=length(pr)-1:-1:2
    p1 = pr[i]
    p2 = pr[i-1]
    f1 = f^p1
    f2 = f^p2
    # 1 -> 1+f^p1/1+f^p2 -> 1+f/1+f^p2 -> 1+f/1+f^p1 -> 1
    # by induction, we have  presentation for 1+f/1+f^p1
    # gens are polys and rels is a matrix with the relations
    #
    # new gens:
    # 1+f^p1/1+f^p2 = f^p1/f^p2 = f^(p2-p1), latter additively
    ngens = [1+(x^i)*f1*c for i=0:(degree(f)*(p2-p1)-1) for c = b]
    nr = matrix(FlintZZ, 0, length(ngens), [])
    for j=1:nrows(rels)
      g = rem(prod(powmod(gens[k], rels[j, k], f2) for k=1:ncols(rels)), f2) - 1
      q,r = divrem(g, f1)
      @assert iszero(r)
      nr = vcat(nr, matrix(FlintZZ, 1, ncols(nr), [(coeff(coeff(q, k), l)) for k = 0:degree(f)*(p2-p1)-1 for l = 0:degree(K)-1]))
    end
    rels = [rels -nr; zero_matrix(FlintZZ, length(ngens), ncols(rels)) identity_matrix(FlintZZ, length(ngens))*p]
    append!(gens, ngens)
  end
  return gens, rels
end

#=
function FlintFiniteField(f::fq_nmod_poly, s::AbstractString = "o"; cached::Bool = true, check::Bool = true)
  if check && !isirreducible(f)
    error("poly not irreducible")
  end
  k = base_ring(f)
  p = characteristic(k)
  K, o = FlintFiniteField(p, degree(k)*degree(f), s, cached = cached)
  r = roots(f, K)[1]  # not working, embeddings are missing
  fl || error("s.th. went wrong")
  return K, r
end
=#

function euler_phi(f::T) where {T <: Union{gfp_poly, fq_nmod_poly, gfp_fmpz_poly}}
  lf = factor(f)
  q = size(base_ring(f))
  return prod((q^degree(p)-1)*q^(degree(p)*k) for (p,k) = lf.fac)
end

function carmichael_lambda(f::T) where {T <: Union{gfp_poly, fq_nmod_poly, gfp_fmpz_poly}}
  lf = factor(f)
  pp = characteristic(base_ring(f))
  q = size(base_ring(f))
  #ala Auer... (Diss, DOI:10.4064/aa-95-2-97-122)
  l = reduce(lcm, [(q^degree(p)-1)*pp^ceil(Int, log(k)/log(pp)) for (p,k) = lf.fac], init = fmpz(1))
  #l = reduce(lcm, [(q^degree(p)-1)*largest_elementary_divisor(unit_group_1_part(p, k)[2]) for (p,k) = lf.fac], init = fmpz(1))
  return l
end


text:
function find_first_above_or_equal(cutoff, x::StepRangeLen)
    # this doesn't work for all cases - assert the right ones
    @assert Float64(x.step) > 0
    @assert x.offset == 1
    # basic math:
    # y = x.ref + (i-x.offset)*x.step
    # i = (y - x.ref)/x.step + x.offset
    ind_float = (cutoff - x.ref)/x.step
    ind = Int(ceil(Float64(ind_float))) + x.offset
    @assert ind <= x.len # none of the elements are above or equal the cutoff
    return max(1, ind)
end

function find_last_below_or_equal(cutoff, x::StepRangeLen)
    # this doesn't work for all cases - assert the right ones
    @assert Float64(x.step) > 0
    @assert x.offset == 1
    # basic math:
    # y = x.ref + (i-x.offset)*x.step
    # i = (y - x.ref)/x.step + x.offset
    ind_float = (cutoff - x.ref)/x.step
    ind = Int(floor(Float64(ind_float))) + x.offset
    @assert ind >= 1 # none of the elements are below or equal the cutoff
    return min(x.len, ind)
end

function find_last_below(cutoff, x::StepRangeLen)
    # this doesn't work for all cases - assert the right ones
    @assert Float64(x.step) > 0
    @assert x.offset == 1
    # basic math:
    # y = x.ref + (i-x.offset)*x.step
    # i = (y - x.ref)/x.step + x.offset
    ind_float = (cutoff - x.ref)/x.step
    ind = Int(floor(Float64(ind_float)|>prevfloat)) + x.offset
    @assert ind >= 1 # none of the elements are below or equal the cutoff
    return min(x.len, ind)
end

ext: jl


text:
export
    Circ,
    AABB,
    AABB_center_length_width,
    OBB

struct Circ{V<:AbstractVec}
    c::V # center
    r::Float64 # radius
end
Circ(x::Real, y::Real, r::Real) = Circ{VecE2}(VecE2(x, y), r)
Circ(x::Real, y::Real, z::Real, r::Real) = Circ{VecE3}(VecE3(x, y, z), r)

Base.:+(circ::Circ{VecE2}, v::VecE2) = Circ{VecE2}(circ.c + v, circ.r)
Base.:-(circ::Circ{VecE2}, v::VecE2) = Circ{VecE2}(circ.c - v, circ.r)

Base.contains(circ::Circ{V}, p::V) where {V} = normsquared(circ.c - p) ≤ circ.r*circ.r
inertial2body(circ::Circ{VecE2}, reference::VecSE2) = Circ{VecE2}(inertial2body(circ.c, reference), circ.r)

# Axis-Aligned Bounding Box
struct AABB
    center::VecE2
    len::Float64 # length along x axis
    wid::Float64 # width along y axis
end
function AABB(bot_left::VecE2, top_right::VecE2)
    center = (bot_left + top_right)/2
    Δ = top_right - bot_left
    return AABB(center, abs(Δ.x), abs(Δ.y))
end

Base.:+(box::AABB, v::VecE2) = AABB(box.center + v, box.len, box.wid)
Base.:-(box::AABB, v::VecE2) = AABB(box.center - v, box.len, box.wid)

function Base.contains(box::AABB, P::VecE2)
    -box.len/2 ≤ P.x - box.center.x ≤ box.len/2 &&
    -box.wid/2 ≤ P.y - box.center.y ≤ box.wid/2
end

# Oriented Bounding Box
struct OBB
    aabb::AABB
    θ::Float64
end
function OBB(center::VecE2, len::Float64, wid::Float64, θ::Float64)
    del = VecE2(len/2, wid/2)
    return OBB(AABB(center, len, wid), θ)
end
OBB(center::VecSE2, len::Float64, wid::Float64) = OBB(convert(VecE2, center), len, wid, center.θ)

Base.:+(box::OBB, v::VecE2) = OBB(box.aabb+v, box.θ)
Base.:-(box::OBB, v::VecE2) = OBB(box.aabb-v, box.θ)

function Base.contains(box::OBB, P::VecE2)
    C = box.aabb.center
    contains(box.aabb, rot(P-C, -box.θ)+C)
end

function inertial2body(box::OBB, reference::VecSE2)
    c = VecSE2(box.aabb.center, box.θ)
    c′ = inertial2body(c, reference)
    OBB(c′, box.aabb.len, box.aabb.wid)
end


text:
# ---
# title: 1138. Alphabet Board Path
# id: problem1138
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, String
# link: <https://leetcode.com/problems/alphabet-board-path/description/>
# hidden: true
# ---
#
# On an alphabet board, we start at position `(0, 0)`, corresponding to
# character `board[0][0]`.
#
# Here, `board = ["abcde", "fghij", "klmno", "pqrst", "uvwxy", "z"]`, as shown
# in the diagram below.
#
# ![](https://assets.leetcode.com/uploads/2019/07/28/azboard.png)
#
# We may make the following moves:
#
#   * `'U'` moves our position up one row, if the position exists on the board;
#   * `'D'` moves our position down one row, if the position exists on the board;
#   * `'L'` moves our position left one column, if the position exists on the board;
#   * `'R'` moves our position right one column, if the position exists on the board;
#   * `'!'` adds the character `board[r][c]` at our current position `(r, c)` to the answer.
#
# (Here, the only positions that exist on the board are positions with letters
# on them.)
#
# Return a sequence of moves that makes our answer equal to `target` in the
# minimum number of moves.  You may return any path that does so.
#
#
#
# **Example 1:**
#
#
#
#     Input: target = "leet"
#     Output: "DDR!UURRR!!DDD!"
#
#
# **Example 2:**
#
#
#
#     Input: target = "code"
#     Output: "RR!DDRR!UUL!R!"
#
#
#
#
# **Constraints:**
#
#   * `1 <= target.length <= 100`
#   * `target` consists only of English lowercase letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "shah" begin
    img = generate_surface(SynthSphere(), 0.5, [0,0,1], radius = 50, background = false)
    Z, p, q = retrieve_surface(Shah(), img)
    @test Z[76,76] ≈ 1.094239400597376e20
    @test Z[1,1] ≈ 0.35149454483014475
    @test Z[151,151] ≈ 31.010401226974405
    @test Z[90,90] ≈ 1.0007532507217835e20
    @test Z[101,36] ≈ 3.332266322735267e19
    @test Z[22, 136] ≈ 46.397279917579205
    @test Z[42,42] ≈ 2.4242786395625816e19
end


text:
###
# 2048

import Base: rand, <, ==

# =======================================================================
# This code implements 2048 with a state type of UInt64.
# This representations fits 16 tiles using 16 nibbles (4 bits) each.
# Based on github.com/nneonneo/2048-ai, by Robert Xiao

const Row = UInt16
const Board = UInt64 # Our 2048 state representation

const TwentyFortyEightAction = UInt8
const NONE = 0xFF

const SA2048 = Tuple{Board,TwentyFortyEightAction}
const AS2048 = Tuple{TwentyFortyEightAction, Board}

# Move tables.
# Each extracted row or column has 4 nibbles, so 2^16 = 65536 possible values.
# We can thus store lookups ahead of time.
# The lookup value is 0 if there is no move, and otherwise equals a value that can easily be
# xor'ed into the current board state to update the board.
# Rows tables map rows to rows, whereas column tables map colums to expanded columns for efficiency.
const ROW_LEFT_TABLE = Vector{Row}(undef, 2^16)
const ROW_RIGHT_TABLE = Vector{Row}(undef, 2^16)
const COL_UP_TABLE = Vector{Board}(undef, 2^16)
const COL_DOWN_TABLE = Vector{Board}(undef, 2^16)
const SCORE_TABLE = Vector{Float32}(undef, 2^16)

const ROW_MASK = 0x000000000000FFFF
const COL_MASK = 0x000F000F000F000F

# The probability that a two tile is spawned
const PROB_SPAWN_TWO = 0.80f0
const PROB_SPAWN_FOUR = 1.0f0 - PROB_SPAWN_TWO

const LEFT = 0x00
const DOWN = 0x01
const RIGHT = 0x02
const UP = 0x03
const DIRECTIONS = (LEFT, DOWN, RIGHT, UP)
const TWENTY_FORTY_EIGHT_MOVE_STRINGS = ["LEFT", "DOWN", "RIGHT", "UP"]

@with_kw struct TwentyFortyEight
    γ::Float64 = 1.0
end

function transition(mdp::TwentyFortyEight, s::Board, a::TwentyFortyEightAction)
    s′ = move(s, a)
    if s′ == s # terminal state or illegal action
        return s′
    end
    s′ = insert_tile_rand(s′, draw_tile())
    return s′
end

function reward(mdp::TwentyFortyEight, s::Board, a::TwentyFortyEightAction)
    s′ = move(s, a)
    if s′ == s # terminal state or illegal action
        return -1.0
    end
    s′ = insert_tile_rand(s′, draw_tile())
    return score_board(s′) - score_board(s)
end

function transition_and_reward(mdp::TwentyFortyEight, s::Board, a::TwentyFortyEightAction)
    s′ = move(s, a)
    if s′ == s # terminal state or illegal action
        return (s′, -1.0)
    end
    s′ = insert_tile_rand(s′, draw_tile())
    r = score_board(s′) - score_board(s)
    return (s′, r)
end

function MDP(mdp::TwentyFortyEight; γ::Float64=mdp.γ)
    return MDP(
            γ,
            nothing, # no ordered states
            DIRECTIONS,
            nothing, # no probabilistic transition function
            (s,a) -> reward(mdp, s, a),
            (s, a)-> transition_and_reward(mdp, s,a)
        )
end

"""
Print out a 2048 state.
"""
function print_board(s::Board)
    for i in 1:4
        for j in 1:4
            v = (s & 0x0F) # the tile value is 2^v, where v is the nibble value
            @printf("%6u", v == 0 ? 0 : 1 << v)
            s >>= 4
        end
        println("")
    end
    println("")
end

"""
Turn a column back into a full-sized board
"""
function unpack_col(row::Row)
    x = Board(row)
    return (x | (x << 12) | (x << 24) | (x << 36)) & COL_MASK
end

"""
Reverse a 2048 row
"""
function reverse_row(row::Row)
    return (row >> 12) | ((row >> 4) & 0x00F0)  | ((row << 4) & 0x0F00) | (row << 12)
end

"""
Transpose rows/columns in a board:
  0123       048c
  4567  -->  159d
  89ab       26ae
  cdef       37bf
This method is used to easily extract columns.
"""
function transpose_board(s::Board)
    # Set a to:
    #     0426
    #     1537
    #     8cae
    #     9dbf
    a1 = s & 0xF0F00F0FF0F00F0F
    a2 = s & 0x0000F0F00000F0F0
    a3 = s & 0x0F0F00000F0F0000
    a = a1 | (a2 << 12) | (a3 >> 12)

    b1 = a & 0xFF00FF0000FF00FF
    b2 = a & 0x00FF00FF00000000
    b3 = a & 0x00000000FF00FF00
    return b1 | (b2 >> 24) | (b3 << 24)
end

"""
Count the number of empty tiles (nibbles that are zero) in a state.
Precondition: the board cannot be fully empty.
"""
function count_empty(s::Board)
    # Any set upper two bits are copied down to lower two bits, per nibble
    s |= (s >> 2) & 0x3333333333333333
    # Any set 2nd bit is copied down to the lowest bit, per nibble
    s |= (s >> 1)
    # Only keep first bits that are not set
    s = ~s & 0x1111111111111111
    # At this point each nibble is:
    #  0 if the original nibble was non-zero
    #  1 if the original nibble was zero
    # Next sum them all
    s += s >> 32
    s += s >> 16
    s += s >>  8
    s += s >>  4 # This can overflow to the next nibble if there were 16 empty positions
    return s & 0x0F
end

"""
Populate the score and move tables, which are used for quick lookups.
Note - nneoneo's heuristic score calculation has been removed.
"""
function init_2048_tables()

    for row in 0x0000 : 0xFFFF
        line = [
            (row >>  0) & 0x0F,
            (row >>  4) & 0x0F,
            (row >>  8) & 0x0F,
            (row >> 12) & 0x0F,
        ]

        # Score
        score = 0.0f0
        for rank in line
            if rank >= 2
                # The score is the total sum of the tile and all intermediate merged tiles
                score += (rank - 1) * (1 << rank)
            end
        end
        SCORE_TABLE[row+1] = score

        # Execute a move to the left
        i = 1
        while i < 4
            j = i + 1
            while j < 5
                if line[j] != 0
                    break
                end
                j += 1
            end
            if j == 5
                # No more tiles to the right
                break
            end

            if line[i] == 0
                line[i] = line[j]
                line[j] = 0
                i -= 1 # retry this entry
            elseif line[i] == line[j]
                if line[i] != 0x0F
                    # Note: Pretend that 32768 + 32768 = 32768 (representational limit).
                    line[i] += 1
                end
                line[j] = 0
            end
            i += 1
        end

        result = (line[1] <<  0) |
                 (line[2] <<  4) |
                 (line[3] <<  8) |
                 (line[4] << 12)
        rev_result = reverse_row(result)
        rev_row = reverse_row(row)

        ROW_LEFT_TABLE[     row+1] =                row  ⊻                result
        ROW_RIGHT_TABLE[rev_row+1] =            rev_row  ⊻            rev_result
        COL_UP_TABLE[       row+1] = unpack_col(    row) ⊻ unpack_col(    result)
        COL_DOWN_TABLE[ rev_row+1] = unpack_col(rev_row) ⊻ unpack_col(rev_result)
    end
end
init_2048_tables() # Initialize our score and move tables

function move_up(s::Board)
    s′ = s
    t = transpose_board(s)
    s′ ⊻= COL_UP_TABLE[((t >>  0) & ROW_MASK)+1] <<  0
    s′ ⊻= COL_UP_TABLE[((t >> 16) & ROW_MASK)+1] <<  4
    s′ ⊻= COL_UP_TABLE[((t >> 32) & ROW_MASK)+1] <<  8
    s′ ⊻= COL_UP_TABLE[((t >> 48) & ROW_MASK)+1] << 12
    return s′
end

function move_down(s::Board)
    s′ = s
    t = transpose_board(s)
    s′ ⊻= COL_DOWN_TABLE[((t >>  0) & ROW_MASK)+1] <<  0
    s′ ⊻= COL_DOWN_TABLE[((t >> 16) & ROW_MASK)+1] <<  4
    s′ ⊻= COL_DOWN_TABLE[((t >> 32) & ROW_MASK)+1] <<  8
    s′ ⊻= COL_DOWN_TABLE[((t >> 48) & ROW_MASK)+1] << 12
    return s′
end

function move_left(s::Board)
    s′ = s
    s′ ⊻= Board(ROW_LEFT_TABLE[((s >>  0) & ROW_MASK)+1]) <<  0
    s′ ⊻= Board(ROW_LEFT_TABLE[((s >> 16) & ROW_MASK)+1]) << 16
    s′ ⊻= Board(ROW_LEFT_TABLE[((s >> 32) & ROW_MASK)+1]) << 32
    s′ ⊻= Board(ROW_LEFT_TABLE[((s >> 48) & ROW_MASK)+1]) << 48
    return s′
end

function move_right(s::Board)
    s′ = s
    s′ ⊻= Board(ROW_RIGHT_TABLE[((s >>  0) & ROW_MASK)+1]) <<  0
    s′ ⊻= Board(ROW_RIGHT_TABLE[((s >> 16) & ROW_MASK)+1]) << 16
    s′ ⊻= Board(ROW_RIGHT_TABLE[((s >> 32) & ROW_MASK)+1]) << 32
    s′ ⊻= Board(ROW_RIGHT_TABLE[((s >> 48) & ROW_MASK)+1]) << 48
    return s′
end

"""
This is our deterministic transition function
"""
function move(s::Board, a::TwentyFortyEightAction)
    if a == UP
        return move_up(s)
    elseif a == DOWN
        return move_down(s)
    elseif a == LEFT
        return move_left(s)
    elseif a == RIGHT
        return move_right(s)
    else
        return ~zero(UInt64) # Invalid move - return a filled board
    end
end

"""
True if the 2048 state has no valid actions.
"""
function isdone(s::Board)
    return move(s, UP   ) == s &&
           move(s, DOWN ) == s &&
           move(s, LEFT ) == s &&
           move(s, RIGHT) == s
end

"""
Get the maximum rank among tiles on the board.
"""
function get_max_rank(s::Board)
    maxrank = 0
    while s > 0
        maxrank = max(maxrank, Int(s & 0x0f))
        s >>= 4
    end
    return maxrank
end

"""
Get the number of distinct tile values.
"""
function count_distinct_tiles(s::Board)
    bitset = 0x00
    while s > 0
        bitset |= 1 << (s & 0x0f)
        s >>= 4
    end

    # Do not count empty tiles.
    bitset >>= 1

    n_distinct = 0
    while bitset > 0
        bitset &= bitset - 1
        n_distinct += 1
    end
    return n_distinct
end

"""
Randomly determine whether to spawn a 2 or a 4 tile.
"""
draw_tile()::Board = rand() < PROB_SPAWN_TWO ? 1 : 2

"""
Spawn a tile in a random square, which occurs after we move.
"""
function insert_tile_rand(s::Board, tile::Board)
    n_empty = count_empty(s)
    @assert n_empty > 0
    index = rand(1:n_empty)
    tmp = s
    while true
        while (tmp & 0x0F) != 0
            tmp >>= 4
            tile <<= 4
        end
        if index == 1
            break
        end
        index -= 1
        tmp >>= 4
        tile <<= 4
    end
    return s | tile
end

"""
Generate a starting board.
"""
function initial_board()
    s = draw_tile() << (4 * rand(0:15))
    return insert_tile_rand(s, draw_tile())
end

"""
A helper function for scoring boards.
Scoring is done row-wise using our score table.
"""
function score_helper(s::Board, table::Vector{Float32})
    return table[((s >>  0) & ROW_MASK)+1] +
           table[((s >> 16) & ROW_MASK)+1] +
           table[((s >> 32) & ROW_MASK)+1] +
           table[((s >> 48) & ROW_MASK)+1]
end

"""
Get the true score of a single board.
Note that we cannot tell the difference between a 4 that was spawned
and a 4 that was merged from two 2s. This scoring function assumes
all 4s were spawned from 2s.
"""
score_board(s::Board) = score_helper(s, SCORE_TABLE)

"""
Play 2048 randomly to completition from the given initial board state.
The final board score is returned.
"""
function rollout_to_end(s::Board)
    a = NONE
    n_valid = 0
    for a_candidate in DIRECTIONS
        s′ = move(s, a_candidate)
        if s′ != s
            # valid move; potentially accept
            n_valid += 1
            if rand() < 1/n_valid
                a = a_candidate
            end
        end
    end

    # No valid moves
    if a == NONE
        return score_board(s)
    end

    s′ = insert_tile_rand(move(s,a), draw_tile())
    return rollout_to_end(s′)
end

"""
Play 2048 to completion using the given policy.
The final score is returned.
Note that this core is "correct" in that we track whether 2 or 4 tiles are generated
and update the score appropriately.
"""
function play_game(π::Function)
    s = initial_board()

    # Number of moves.
    moveno = 0

    # Cumulative penalty for obtaining free 4 tiles, as
    # when computing the score of merged tiles we cannot distinguish between
    # merged 2-tiles and spawned 4 tiles.
    scorepenalty = score_board(s)

    while !isdone(s)

        moveno += 1
        println("Move #$(moveno), current score=$(score_board(s) - scorepenalty)")
        print_board(s)

        a = π(s)
        if a == NONE
            break
        end

        println("\ta = ", TWENTY_FORTY_EIGHT_MOVE_STRINGS[a+1])

        s′ = move(s, a)
        if s′ == s
            @warn "Illegal move!"
            moveno -= 1
            continue
        end

        tile = draw_tile()
        if tile == 2
            scorepenalty += 4
        end
        s = insert_tile_rand(s′, tile)
    end

    print_board(s)
    @printf("\nGame over. Your score is %.0f. The highest rank you achieved was %d.\n", score_board(s) - scorepenalty, get_max_rank(s));
end


text:
push!(empty!(LOAD_PATH), dirname(dirname(@__DIR__)))
Base.require(:Pkg3)


text:
"""
	plan_PSF(F::LinearMap, M::LinearMap [, S::LinearMap])

"""
plan_psf(F::LinearMap, M::LinearMap) = F' * M * F
plan_psf(F::LinearMap, M::LinearMap, S::LinearMap) = S' * F' * M * F * S

max_stars_repo_name: felixhorger/MRIRecon

max_issues_repo_name: felixhorger/MRIRecon

max_forks_repo_name: felixhorger/MRIRecon


text:
import Base: sum, prod, maximum, minimum, sumabs, sumabs2

cuda21 = [
("add","sum","ai+xi","xi","0"),
("mul","prod","ai*xi","xi","1"),
("max","maximum","(ai>xi?ai:xi)","xi","(-INFINITY)"),
("min","minimum","(ai<xi?ai:xi)","xi","INFINITY"),
("sumabs","sumabs","ai+xi","(xi<0?-xi:xi)","0"),
("sumabs2","sumabs2","ai+xi","(xi*xi)","0"),
]

function cuda21def(f, j=f, o...)
    J=Symbol(j)
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$(S)_21"
        @eval begin
            function $J(x::KnetArray{$T}, region)
                if length(region) == ndims(x)
                    return fill!(similar(x,ntuple(i->1,ndims(x))), $J(x))
                elseif length(region) == ndims(x)-1
                    i0 = 0
                    ysize = ntuple(ndims(x)) do i
                        if in(i,region)
                            1
                        elseif i0>0
                            error("Bad region $region")
                        else
                            i0=i
                            size(x,i)
                        end
                    end
                    y = similar(x, ysize)
                    nx = length(x); ny = length(y); sy = stride(x,i0)
                    ccall(($F,$libknet8),Void,(Cint,Ptr{$T},Cint,Cint,Ptr{$T}),nx,x,sy,ny,y)
                    return y
                else
                    error("Only scalar and vector reductions supported.")
                end
            end
        end
    end
end

#if isdefined(:libknet8)
    for f in cuda21
        isa(f,Tuple) || (f=(f,))
        cuda21def(f...)
    end
#end


text:
using MathGL
mgl = MathGL

gr = mgl.Graph(800, 500)
@mglplot gr [
    xlabel "x"
    ylabel "y"
    box
    axis
    grid
    plot 0.8*sin(linspace(-4pi, 4pi, 200))
    write "../graphs/sin_mgl.png"
]


text:
# Functions for representing graphs in GraphViz's dot format
# http://www.graphviz.org/
# http://www.graphviz.org/Documentation/dotguide.pdf
# http://www.graphviz.org/pub/scm/graphviz2/doc/info/lang.html

# based on :
# https://github.com/JuliaAttic/OldGraphs.jl/blob/master/src/dot.jl
# https://github.com/tkf/ShowGraphviz.jl
# and modified by U.Roettgermann, 12/21

# internal helper functions
to_dot(sym::String, value::Any) = "$sym=$value"
edge_op(graph::AbstractSimpleWeightedGraph) = Graphs.is_directed(graph) ? "->" : "--"
parse_attributes(attributes) = string("[", join(map(a -> to_dot(a.key, a.value), attributes), ", "))
is_all_zero(arr) = length(arr) == 0 || all(==(0), arr)

function dot(g::AbstractSimpleWeightedGraph, stream::IO, attrs::GraphvizAttributes)

    # write beginning: graph or digraph:
    directed = Graphs.is_directed(g)

    graph_type_string = directed ? "digraph" : "graph"
    write(stream, "$graph_type_string G {\n")

    set!(attrs.plot_options, "type", graph_type_string)
    set!(attrs.graph_options, "concentrate", "true")
    set!(attrs.graph_options, "layout", (directed) ? "dot" : "neato")

    # write GENERAL attributes, belongs to all elements:
    write(stream, " graph$(parse_attributes(attrs.graph_options))]\n")
    write(stream, " node$(parse_attributes(attrs.node_options))]\n")
    write(stream, " edge$(parse_attributes(attrs.edge_options))]\n")

    # write subgraphs:
    if !isempty(attrs.subgraphs)
        for subgraph in attrs.subgraphs
            write(stream, "subgraph $(subgraph.type) { \n")

            write(stream, " graph$(parse_attributes(subgraph.graph_options))]\n")
            write(stream, " node$(parse_attributes(subgraph.node_options))]\n")
            write(stream, " edge$(parse_attributes(subgraph.edge_options))]\n")

            # write node
            for n in subgraph.nodes
                write(stream, " $(n.id) $(parse_attributes(n.attributes))];\n")
            end

            # write edges
            for e in subgraph.edges
                write(stream, " $(e.from) $(edge_op(g)) $(e.to) $(parse_attributes(e.attributes))];\n")
            end

            write(stream, "}\n")
        end
    end

    # write node
    for n in attrs.nodes
        write(stream, " $(n.id) $(parse_attributes(n.attributes))];\n")
    end

    # write edges
    for e in attrs.edges
        write(stream, " $(e.from) $(edge_op(g)) $(e.to) $(parse_attributes(e.attributes))];\n")
    end

    write(stream, "}\n")
    return stream
end

function color_nodes!(attrs::GraphvizAttributes, colors)

    # hard-coded color scheme:
    set!(attrs.node_options, "colorscheme", "set19")

    # standard colorscheme (max 9)
    if (!is_all_zero(colors))  # ToDo: no hardcoded part!!!!
        if maximum(colors) > 9
            colors = _reduce_colors!(colors)
        end
    end

    for node in attrs.nodes
        if colors[node.id] > 0
            set!(node.attributes, "color", colors[node.id])
        end
    end
end

function color_path!(attrs::GraphvizAttributes, path, g::AbstractSimpleWeightedGraph; color = "red")

    for node = 1:nv(g)
        childs = Graphs.inneighbors(g, node)

        if !Base.isnothing(findfirst(isequal(node), path))
            set!(attrs.nodes, node, Property("style", "filled"))
            set!(attrs.nodes, node, Property("fillcolor", color))
            for kid in childs
                if !Base.isnothing(findfirst(isequal(kid), path)) && (kid != path[1])
                    set!(attrs.edges, node, kid, Property("color", color))
                end
            end
        end
    end
end

# helper function to reduze colors
function _reduce_colors!(components)
    n = length(components)
    cz = zeros(Int, n)
    co = deepcopy(components)
    min = 1
    while true
        for i = 1:n
            cz[i] = length(findall(x -> x == co[i], co))
        end
        idx = findall(x -> x < min, cz)
        co[idx] .= 0

        j = 1
        ma = maximum(co)
        for i = 1:ma
            idx = findall(x -> x == i, co)
            if !isempty(idx)
                co[idx] .= j
                j = j + 1
            end
        end

        # end while
        (maximum(co)<9):break:(min=min+1)
    end

    return components = co
end

max_stars_repo_name: tragisch/PlotGraphviz

max_issues_repo_name: tragisch/PlotGraphviz


text:
abstract ResourceKey

abstract AbstractResource

function isless(a::ResourceKey, b::ResourceKey)
  return (a.priority < b.priority) || (a.priority == b.priority && a.id < b.id)
end

function trigger_put(sim::Simulation, put_ev::Event, res::AbstractResource)
  queue = PriorityQueue(res.put_queue)
  while length(queue) > 0
    (put_ev, key) = peek(queue)
    proceed = do_put(sim, res, put_ev, key)
    if state(put_ev) == triggered
      dequeue!(res.put_queue, put_ev)
    end
    if proceed
      dequeue!(queue)
    else
      break
    end
  end
end

function trigger_get(sim::Simulation, get_ev::Event, res::AbstractResource)
  queue = PriorityQueue(res.get_queue)
  while length(queue) > 0
    (get_ev, key) = peek(queue)
    proceed = do_get(sim, res, get_ev, key)
    if state(get_ev) == triggered
      dequeue!(res.get_queue, get_ev)
    end
    if proceed
      dequeue!(queue)
    else
      break
    end
  end
end


text:
module Bibliography

# BibInternal
import BibInternal
import BibInternal: AbstractEntry
import BibInternal.BibTeX: Article, Book, Booklet, InBook, InCollection, InProceedings, Manual, MasterThesis, Misc, PhDThesis, Proceedings, TechReport, Unpublished

# BibParser
import BibParser, BibParser.BibTeX

# Others
import DataStructures
import DataStructures.OrderedSet

export export_bibtex, import_bibtex
export export_web, bibtex_to_web

include("bibtex.jl")
include("staticweb.jl")

function export_bibtex(target::AbstractString, bibliography::DataStructures.OrderedSet{AbstractEntry})
    data = export_bibtex(bibliography)
    if target != ""
        f = open(target, "w")
        write(f, data)
        close(f)
    end
    return data
end

function bibtex_to_web(source::AbstractString)
    export_web(import_bibtex(source))
end

end # module


text:
#This Julia file has been generated by rosetta script of MatlabCompat library from an m-file
#The code generated need further corrections by you. Execute it line by line and correct the errors.
#
# Janus is a minimalistic model image
# analysis script aimed
# to work in Julia and MATLAB without
# changing its syntax
#
# Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
# MIT license.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including without
# limitation the rights to use, copy, modify, merge, publish, distribute,
# sublicense, and/or sell copies of the Software, and to permit persons to whom
# the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using MatlabCompat
importall MatlabCompat.ImageTools
importall MatlabCompat.MathTools

tic()
img = imread("http://matlabcompat.github.io/img/example.tif"); # read the remote image
threshold = graythresh(img); # compute graysacale threshold using Otsu algorithm
imgbw = im2bw(img, threshold); # create a binary image based on the grayscale image
imshow(imgbw); # display the resulting binary image
labeledbw = bwlabel(imgbw, 4); # label each connected object in the image
numberOfCells = max(reshape(labeledbw, 1,numel(labeledbw))); # count cells
disp(strcat("number of objects:", num2str(numberOfCells)));# display the number of objects
imshow(label2rgb(labeledbw,"jet",[1 1 1],"shuffle"));
toc()


text:
x0 = 10
y0 = 10
r = 1
x = r.*[-1.1,0,0,1] .+ x0
y = r.*[0.0,1,-1,0] .+ y0
result = Circle([x0,y0],r,[x y],:kasa)
coefs = coef(result)
@test x0 ≈ coefs[1]
@test y0 ≈ coefs[2]
@test r ≈ coefs[3]
@test dof(result) == 1
@test all(coef(result) .== [result.position...,result.radius])
@test length(coefnames(result)) == length(coef(result))
@test residuals(result) ≈ [0.1,0,0,0]
@test rss(result) ≈ 0.1^2

x = r.*[-1,0,0,1] .+ x0
result = fit(Circle,x,y)
coefs = coef(result)
@test x0 ≈ coefs[1]
@test y0 ≈ coefs[2]
@test r ≈ coefs[3]
@test algorithm(result) == :kasa

result = fit(Circle,x,y,alg=:pratt)
coefs = coef(result)
@test x0 ≈ coefs[1]
@test y0 ≈ coefs[2]
@test r ≈ coefs[3]
@test algorithm(result) == :pratt

result = fit(Circle,x,y,alg=:taubin)
coefs = coef(result)
@test x0 ≈ coefs[1]
@test y0 ≈ coefs[2]
@test r ≈ coefs[3]
@test algorithm(result) == :taubin


text:
using InteractiveDynamics, GLMakie
using DynamicalSystems

i = 1
p_index = 1

systems = [(Systems.logistic(), 2.5, 4.0, "r", "logistic map"),
           (Systems.henon(), 0.8, 1.4, "a", "Hénon map"),
           (Systems.standardmap(), 0.0, 1.2, "k", "standard map")]

ds, p_min, p_max, parname, t = systems[2]
t = "orbit diagram for the "*t

oddata = interactive_orbitdiagram(ds, p_index, p_min, p_max, i;
                                  parname = parname, title = t)

ps, us = scaleod(oddata)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: service_quotas
using AWS.Compat
using AWS.UUIDs

"""
    associate_service_quota_template()
    associate_service_quota_template(params::Dict{String,<:Any})

Associates your quota request template with your organization. When a new account is
created in your organization, the quota increase requests in the template are automatically
applied to the account. You can add a quota increase request for any adjustable quota to
your template.

"""
associate_service_quota_template(; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("AssociateServiceQuotaTemplate"; aws_config=aws_config)
associate_service_quota_template(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("AssociateServiceQuotaTemplate", params; aws_config=aws_config)

"""
    delete_service_quota_increase_request_from_template(aws_region, quota_code, service_code)
    delete_service_quota_increase_request_from_template(aws_region, quota_code, service_code, params::Dict{String,<:Any})

Deletes the quota increase request for the specified quota from your quota request template.

# Arguments
- `aws_region`: The AWS Region.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
delete_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("DeleteServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
delete_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("DeleteServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    disassociate_service_quota_template()
    disassociate_service_quota_template(params::Dict{String,<:Any})

Disables your quota request template. After a template is disabled, the quota increase
requests in the template are not applied to new accounts in your organization. Disabling a
quota request template does not apply its quota increase requests.

"""
disassociate_service_quota_template(; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("DisassociateServiceQuotaTemplate"; aws_config=aws_config)
disassociate_service_quota_template(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("DisassociateServiceQuotaTemplate", params; aws_config=aws_config)

"""
    get_association_for_service_quota_template()
    get_association_for_service_quota_template(params::Dict{String,<:Any})

Retrieves the status of the association for the quota request template.

"""
get_association_for_service_quota_template(; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetAssociationForServiceQuotaTemplate"; aws_config=aws_config)
get_association_for_service_quota_template(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetAssociationForServiceQuotaTemplate", params; aws_config=aws_config)

"""
    get_awsdefault_service_quota(quota_code, service_code)
    get_awsdefault_service_quota(quota_code, service_code, params::Dict{String,<:Any})

Retrieves the default value for the specified quota. The default value does not reflect any
quota increases.

# Arguments
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
get_awsdefault_service_quota(QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetAWSDefaultServiceQuota", Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_awsdefault_service_quota(QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetAWSDefaultServiceQuota", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    get_requested_service_quota_change(request_id)
    get_requested_service_quota_change(request_id, params::Dict{String,<:Any})

Retrieves information about the specified quota increase request.

# Arguments
- `request_id`: The ID of the quota increase request.

"""
get_requested_service_quota_change(RequestId; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetRequestedServiceQuotaChange", Dict{String, Any}("RequestId"=>RequestId); aws_config=aws_config)
get_requested_service_quota_change(RequestId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetRequestedServiceQuotaChange", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RequestId"=>RequestId), params)); aws_config=aws_config)

"""
    get_service_quota(quota_code, service_code)
    get_service_quota(quota_code, service_code, params::Dict{String,<:Any})

Retrieves the applied quota value for the specified quota. For some quotas, only the
default values are available. If the applied quota value is not available for a quota, the
quota is not retrieved.

# Arguments
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
get_service_quota(QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetServiceQuota", Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_service_quota(QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetServiceQuota", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    get_service_quota_increase_request_from_template(aws_region, quota_code, service_code)
    get_service_quota_increase_request_from_template(aws_region, quota_code, service_code, params::Dict{String,<:Any})

Retrieves information about the specified quota increase request in your quota request
template.

# Arguments
- `aws_region`: The AWS Region.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
get_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
get_service_quota_increase_request_from_template(AwsRegion, QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("GetServiceQuotaIncreaseRequestFromTemplate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AwsRegion"=>AwsRegion, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    list_awsdefault_service_quotas(service_code)
    list_awsdefault_service_quotas(service_code, params::Dict{String,<:Any})

Lists the default values for the quotas for the specified AWS service. A default value does
not reflect any quota increases.

# Arguments
- `service_code`: The service identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
"""
list_awsdefault_service_quotas(ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListAWSDefaultServiceQuotas", Dict{String, Any}("ServiceCode"=>ServiceCode); aws_config=aws_config)
list_awsdefault_service_quotas(ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListAWSDefaultServiceQuotas", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    list_requested_service_quota_change_history()
    list_requested_service_quota_change_history(params::Dict{String,<:Any})

Retrieves the quota increase requests for the specified service.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
- `"ServiceCode"`: The service identifier.
- `"Status"`: The status of the quota increase request.
"""
list_requested_service_quota_change_history(; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistory"; aws_config=aws_config)
list_requested_service_quota_change_history(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistory", params; aws_config=aws_config)

"""
    list_requested_service_quota_change_history_by_quota(quota_code, service_code)
    list_requested_service_quota_change_history_by_quota(quota_code, service_code, params::Dict{String,<:Any})

Retrieves the quota increase requests for the specified quota.

# Arguments
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
- `"Status"`: The status value of the quota increase request.
"""
list_requested_service_quota_change_history_by_quota(QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistoryByQuota", Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
list_requested_service_quota_change_history_by_quota(QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListRequestedServiceQuotaChangeHistoryByQuota", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    list_service_quota_increase_requests_in_template()
    list_service_quota_increase_requests_in_template(params::Dict{String,<:Any})

Lists the quota increase requests in the specified quota request template.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AwsRegion"`: The AWS Region.
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
- `"ServiceCode"`: The service identifier.
"""
list_service_quota_increase_requests_in_template(; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListServiceQuotaIncreaseRequestsInTemplate"; aws_config=aws_config)
list_service_quota_increase_requests_in_template(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListServiceQuotaIncreaseRequestsInTemplate", params; aws_config=aws_config)

"""
    list_service_quotas(service_code)
    list_service_quotas(service_code, params::Dict{String,<:Any})

Lists the applied quota values for the specified AWS service. For some quotas, only the
default values are available. If the applied quota value is not available for a quota, the
quota is not retrieved.

# Arguments
- `service_code`: The service identifier.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
"""
list_service_quotas(ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListServiceQuotas", Dict{String, Any}("ServiceCode"=>ServiceCode); aws_config=aws_config)
list_service_quotas(ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListServiceQuotas", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    list_services()
    list_services(params::Dict{String,<:Any})

Lists the names and codes for the services integrated with Service Quotas.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of results to return with a single call. To retrieve
  the remaining results, if any, make another call with the token returned from this call.
- `"NextToken"`: The token for the next page of results.
"""
list_services(; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListServices"; aws_config=aws_config)
list_services(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListServices", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Returns a list of the tags assigned to the specified applied quota.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the applied quota for which you want
  to list tags. You can get this information by using the Service Quotas console, or by
  listing the quotas using the list-service-quotas AWS CLI command or the ListServiceQuotas
  AWS API operation.

"""
list_tags_for_resource(ResourceARN; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListTagsForResource", Dict{String, Any}("ResourceARN"=>ResourceARN); aws_config=aws_config)
list_tags_for_resource(ResourceARN, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("ListTagsForResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARN"=>ResourceARN), params)); aws_config=aws_config)

"""
    put_service_quota_increase_request_into_template(aws_region, desired_value, quota_code, service_code)
    put_service_quota_increase_request_into_template(aws_region, desired_value, quota_code, service_code, params::Dict{String,<:Any})

Adds a quota increase request to your quota request template.

# Arguments
- `aws_region`: The AWS Region.
- `desired_value`: The new, increased value for the quota.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
put_service_quota_increase_request_into_template(AwsRegion, DesiredValue, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("PutServiceQuotaIncreaseRequestIntoTemplate", Dict{String, Any}("AwsRegion"=>AwsRegion, "DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
put_service_quota_increase_request_into_template(AwsRegion, DesiredValue, QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("PutServiceQuotaIncreaseRequestIntoTemplate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AwsRegion"=>AwsRegion, "DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    request_service_quota_increase(desired_value, quota_code, service_code)
    request_service_quota_increase(desired_value, quota_code, service_code, params::Dict{String,<:Any})

Submits a quota increase request for the specified quota.

# Arguments
- `desired_value`: The new, increased value for the quota.
- `quota_code`: The quota identifier.
- `service_code`: The service identifier.

"""
request_service_quota_increase(DesiredValue, QuotaCode, ServiceCode; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("RequestServiceQuotaIncrease", Dict{String, Any}("DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode); aws_config=aws_config)
request_service_quota_increase(DesiredValue, QuotaCode, ServiceCode, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("RequestServiceQuotaIncrease", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DesiredValue"=>DesiredValue, "QuotaCode"=>QuotaCode, "ServiceCode"=>ServiceCode), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds tags to the specified applied quota. You can include one or more tags to add to the
quota.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the applied quota. You can get this
  information by using the Service Quotas console, or by listing the quotas using the
  list-service-quotas AWS CLI command or the ListServiceQuotas AWS API operation.
- `tags`: The tags that you want to add to the resource.

"""
tag_resource(ResourceARN, Tags; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("TagResource", Dict{String, Any}("ResourceARN"=>ResourceARN, "Tags"=>Tags); aws_config=aws_config)
tag_resource(ResourceARN, Tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("TagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARN"=>ResourceARN, "Tags"=>Tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes tags from the specified applied quota. You can specify one or more tags to remove.

# Arguments
- `resource_arn`: The Amazon Resource Name (ARN) for the applied quota that you want to
  untag. You can get this information by using the Service Quotas console, or by listing the
  quotas using the list-service-quotas AWS CLI command or the ListServiceQuotas AWS API
  operation.
- `tag_keys`: The keys of the tags that you want to remove from the resource.

"""
untag_resource(ResourceARN, TagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("UntagResource", Dict{String, Any}("ResourceARN"=>ResourceARN, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resource(ResourceARN, TagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = service_quotas("UntagResource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARN"=>ResourceARN, "TagKeys"=>TagKeys), params)); aws_config=aws_config)


text:
add_brackets(arr::AbstractArray, vars) = [add_brackets(element, vars) for element in arr]

add_brackets(s::Any, vars) = return s

function add_brackets(ex::Expr, vars)
    ex = postwalk(x -> x in vars ? "\\left[ $(convertSubscript(x)) \\right]" : x, ex)
    return ex
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

module AbaqusReader

using Nullables

include("parse_mesh.jl")
include("keyword_register.jl")
include("parse_model.jl")
include("create_surface_elements.jl")
include("abaqus_download.jl")

export abaqus_read_mesh, abaqus_read_model, create_surface_elements
export abaqus_download

end


text:
"""

SAVE AND READ SAMPLES

"""

"""
    writesamples(data::Matrix{Int},
                 [model::Union{MPS, MPO, LPDO, Nothing},]
                 output_path::String)

Save data and model on file:

# Arguments:
  - `data`: array of measurement data
  - `model`: (optional) MPS, MPO, or Choi
  - `output_path`: path to file
"""
function writesamples(
  data::Matrix{Int}, model::Union{MPS,MPO,LPDO,Nothing}, output_path::String
)
  # Make the path the file will sit in, if it doesn't exist
  mkpath(dirname(output_path))
  h5rewrite(output_path) do fout
    write(fout, "outcomes", data)
    if isnothing(model)
      write(fout, "state", "nothing")
    else
      write(fout, "state", model)
    end
  end
end

function writesamples(data::Matrix{Int}, output_path::String)
  # Make the path the file will sit in, if it doesn't exist
  mkpath(dirname(output_path))
  h5rewrite(output_path) do fout
    write(fout, "outcomes", data)
  end
end

function writesamples(
  data::Matrix{Pair{String,Int}}, model::Union{MPS,MPO,LPDO,Nothing}, output_path::String
)
  # Make the path the file will sit in, if it doesn't exist
  mkpath(dirname(output_path))
  h5rewrite(output_path) do fout
    write(fout, "bases", first.(data))
    write(fout, "outcomes", last.(data))
    if isnothing(model)
      write(fout, "state", "nothing")
    else
      write(fout, "state", model)
    end
  end
end

function writesamples(data::Matrix{Pair{String,Int}}, output_path::String)
  # Make the path the file will sit in, if it doesn't exist
  mkpath(dirname(output_path))
  h5rewrite(output_path) do fout
    write(fout, "bases", first.(data))
    write(fout, "outcomes", last.(data))
  end
end

function writesamples(
  data::Matrix{Pair{String,Pair{String,Int}}},
  model::Union{MPS,MPO,LPDO,Nothing},
  output_path::String,
)
  # Make the path the file will sit in, if it doesn't exist
  mkpath(dirname(output_path))
  h5rewrite(output_path) do fout
    write(fout, "inputs", first.(data))
    write(fout, "bases", first.(last.(data)))
    write(fout, "outcomes", last.(last.(data)))
    if isnothing(model)
      write(fout, "state", "nothing")
    else
      write(fout, "state", model)
    end
  end
end

function writesamples(data::Matrix{Pair{String,Pair{String,Int}}}, output_path::String)
  # Make the path the file will sit in, if it doesn't exist
  mkpath(dirname(output_path))
  h5rewrite(output_path) do fout
    write(fout, "inputs", first.(data))
    write(fout, "bases", first.(last.(data)))
    write(fout, "outcomes", last.(last.(data)))
  end
end

"""
    readsamples(input_path::String)

Load data and model from file:

# Arguments:
  - `input_path`: path to file
"""
function readsamples(input_path::String)
  fin = h5open(input_path, "r")
  # Check if the data is for state tomography or process tomography
  # Process tomography
  if haskey(fin, "inputs")
    inputs = read(fin, "inputs")
    bases = read(fin, "bases")
    outcomes = read(fin, "outcomes")
    data = inputs .=> (bases .=> outcomes)
    # Measurements in bases
  elseif haskey(fin, "bases")
    bases = read(fin, "bases")
    outcomes = read(fin, "outcomes")
    data = bases .=> outcomes
    # Measurements in Z basis
  elseif haskey(fin, "outcomes")
    data = read(fin, "outcomes")
  else
    close(fin)
    error(
      "File must contain either \"data\" for quantum state tomography data or \"data_first\" and \"data_second\" for quantum process tomography.",
    )
  end

  # Check if a model is saved, if so read it and return it
  if haskey(fin, "state")
    g = fin["state"]
    if haskey(attributes(g), "type")
      typestring = read(attributes(g)["type"])
      modeltype = eval(Meta.parse(typestring))
      model = read(fin, "state", modeltype)
    else
      model = read(fin, "state")
      if model == "nothing"
        model = nothing
      else
        error("model must be MPS, LPDO, or Nothing")
      end
    end
    close(fin)
    return data, model
  end

  close(fin)
  return data
end

"""
Various printing functionalities
"""

printmetric(name::String, metric::Int) = @printf("%s = %d  ", name, metric)
printmetric(name::String, metric::Float64) = @printf("%s = %-4.4f  ", name, metric)
printmetric(name::String, metric::AbstractArray) = @printf("%s = [...]  ", name)

function printmetric(name::String, metric::Complex)
  if imag(metric) < 1e-8
    @printf("%s = %-4.4f  ", name, real(metric))
  else
    @printf("%s = %.4f±i%-4.4f  ", name, real(metric), imag(metric))
  end
end

function printobserver(observer::Observer, print_metrics::Union{String,AbstractArray}
)
  if !isempty(print_metrics)
    if print_metrics isa String
      printmetric(print_metrics, results(observer, print_metrics)[end])
    else
      for metric in print_metrics
        printmetric(metric, results(observer, metric)[end])
      end
    end
  end
  return
end


text:
# This file is a part of KCenters.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export size_criterion, sqrt_criterion, change_criterion, fun_criterion, log2_criterion, epsilon_criterion, salesman_criterion

"""
    epsilon_criterion(e)

Creates a function that evaluates the stop criterion when the distance between far items achieves the given `e`
"""
function epsilon_criterion(e)
   (dmaxlist, database) -> dmaxlist[end] < e
end

"""
    size_criterion(maxsize)

Creates a function that stops when the number of far items are equal or larger than the given `maxsize`
"""
function size_criterion(maxsize)
   (dmaxlist, database) -> length(dmaxlist) >= maxsize
end

"""
    fun_criterion(fun::Function)

Creates a stop-criterion function that stops whenever the number of far items reaches ``\\lceil fun(|database|)\\rceil``.
Already defined examples:
```julia
    sqrt_criterion() = fun_criterion(sqrt)
    log2_criterion() = fun_criterion(log2)
```
"""
function fun_criterion(fun)
   (dmaxlist, database) -> length(dmaxlist) >= ceil(Int, length(database) |> fun |> round)
end

"""
    sqrt_criterion()

Stops after \$sqrt(n)\$ centers
"""
sqrt_criterion() = fun_criterion(sqrt)

"""
    log2_criterion()

Stops after \$log_2(n)\$ centers
"""
log2_criterion() = fun_criterion(log2)

"""
    change_criterion(tol=0.001, window=3)

Creates a fuction that stops the process whenever the maximum distance converges (averaging `window` far items).
The `tol` parameter defines the tolerance range.
"""
function change_criterion(tol=0.001, window=3)
    mlist = Float64[]
    count = 0.0
    function stop(dmaxlist, database)
        count += dmaxlist[end]

        if length(dmaxlist) % window != 1
            return false
        end
        push!(mlist, count)
        count = 0.0
        if length(dmaxlist) < 2
            return false
        end

        s = abs(mlist[end] - mlist[end-1])
        return s <= tol
    end

    return stop
end

"""
    salesman_criterion()

It creates a function that explores the entire dataset making a full farthest first traversal approximation
"""
function salesman_criterion()
    function stop(dmaxlist, dataset)
        return false
    end
end


text:
# module Bukdu

export routes

export get, post, delete, patch, put
import Base: get
using .System: HaltedError, NotApplicableError, InternalError, SystemController, internal_error, halted_error, not_applicable

const server_info = string("Bukdu/", BUKDU_VERSION)
const routing_verbs = [:get, :post, :delete, :patch, :put, :options]

# HTTP.URIs: queryparams
# HTTP.jl - URIs.jl
function queryparams(q::AbstractString)
    Dict(Deps.HTTP.URIs.unescapeuri(k) => Deps.HTTP.URIs.unescapeuri(v)
        for (k,v) in ([split(e, "=")..., ""][1:2]
            for e in split(q, "&", keepempty=false)))
end

function fetch_query_params(req::Deps.Request)::Vector{Pair{String,Any}}
    params = queryparams(Deps.HTTP.URIs.URI(req.target).query)
    collect(params)
end

function parsed_path_params(route::Route)::Vector{Pair{String,Any}}
    if isempty(route.param_types)
        route.path_params
    else
        map(route.path_params) do (name, val)
            typ = get(route.param_types, Symbol(name), String)
            Pair{String,Any}(name, typ === String ? val : parse(typ, val))
        end
    end
end

function _apply_action(route::Route, conn::Conn)::Tuple{Route, Union{Nothing,Render}}
    if conn.halted
        err = HaltedError("halted on pipelines")
        rou = Route(SystemController, halted_error, route.param_types, route.path_params, route.pipelines)
        obj = halted_error(SystemController(conn, err))
        (rou, obj)
    else
        controller = route.C(conn)
        if applicable(route.action, controller)
            conn.request.response.status = 200
            ret = route.action(controller)
            if ret isa Render
                obj = ret
            else
                obj = render(Julia, ret)
            end
            (route, obj)
        else
            err = NotApplicableError(string(route.action, "(::", route.C, ")"))
            rou = Route(SystemController, not_applicable, route.param_types, route.path_params, route.pipelines)
            obj = not_applicable(SystemController(conn, err))
            (rou, obj)
        end
    end
end

function _catch_internal_error(block, route, conn::Conn)::Tuple{Route, Union{Nothing,Render}}
    try
        block()
    catch ex
        stackframes = stacktrace(catch_backtrace())
        err = InternalError(ex, stackframes)
        rou = Route(SystemController, internal_error, route.param_types, route.path_params, route.pipelines)
        obj = internal_error(SystemController(conn, err))
        (rou, obj)
    end
end

function _proc_request(route::Route, conn::Conn)::Tuple{Route, Union{Nothing,Render}}
    System.catch_request(route, conn) # System
    _catch_internal_error(route, conn) do
    # begin
        _apply_action(route, conn)
    end
end

function _proc_response(route::Route, conn::Conn)
    Plug.Loggers.info_response(conn, (controller=route.C, action=route.action))
    System.catch_response(route, conn) # System
end

function put_response_header_and_body(req::Deps.Request, obj::Union{Nothing, Render})
    push!(req.response.headers, Pair("Server", server_info))
    if obj isa Render
        push!(req.response.headers, Pair("Content-Type", obj.content_type))
        if "HEAD" == req.method
            body = Vector{UInt8}()
        else
            body = Vector{UInt8}(obj.writer(obj.data))
        end
        push!(req.response.headers, Pair("Content-Length", string(length(body))))
        req.response.body = body
    end
end

function request_handler(route::Route, conn::Conn)::NamedTuple{(:got, :resp, :route)}
    (rou, obj) = _proc_request(route, conn)
    put_response_header_and_body(conn.request, obj)
    _proc_response(rou, conn)
    (got=obj isa Render ? obj.data : obj, resp=conn.request.response, route=rou)
end

"""
    routes(block::Function)
"""
function routes(block::Function)
    block()
end

"""
    routes(block::Function, pipe::Symbol)
"""
function routes(block::Function, pipe::Symbol)
    Routing.store[:pipe] = pipe
    block()
    Routing.store[:pipe] = nothing
end

"""
    get(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...)
"""
function get
end

"""
    post(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...)
"""
function post
end

"""
    delete(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...)
"""
function delete
end

"""
    patch(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...)
"""
function patch
end

"""
    put(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...)
"""
function put
end

function head
end

"""
    options(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...)
"""
function options
end

function trace
end

for verb in routing_verbs
    @eval ($verb)(url::String, C::Type{<:ApplicationController}, action, param_types::Pair{Symbol,DataType}...) = Routing.add_route($verb, url, C, action, Dict{Symbol,DataType}(param_types...))
end

for (verb, action) in [(:get, :index), (:post, :create)]
    @eval function ($verb)(f::Function, url::String, param_types::Pair{Symbol,DataType}...)
        $action(c::System.AnonymousController) = f(c.conn)
        Routing.add_route($verb, url, System.AnonymousController, $action, Dict{Symbol,DataType}(param_types...))
    end
end

# module Bukdu


text:
function gradcheck(net;showgrad=false,returngrad=false)
    #ONLY WORKS ON THE CPU AT THE MOMENT
    # (c) David Barber, University College London 2015
    println("gradcheck...")
    value=net.value
    node=net.node
    net2=deepcopy(net) ## CAREFUL -- DEEPCOPY IS BUGGY FOR CUDAARRAY

    ADforward!(net)
    ADbackward!(net)
    g=net.gradient

    epsilon=1e-7;
    tol=1e-6
    gemp=deepcopy(g)
    for par=1:length(net.node) # loop over the nodes in the network
        if node[par].returnderivative
            println("------------------------------------------------------------------------------")
            println("node $par:")
            for linearidx in 1:prod(size(value[par]))
                net2.value[par][linearidx]=value[par][linearidx]+epsilon
                fplus=ADforward!(net2,returnf=true)
                net2.value[par][linearidx]=value[par][linearidx]-epsilon
                fminus=ADforward!(net2,returnf=true)
                net2.value[par][linearidx]=value[par][linearidx]
                gemp[par][linearidx]=0.5*(fplus-fminus)/epsilon
            end
            diff=mean(abs(g[par]-gemp[par]))
            reldiff=mean(abs((g[par]-gemp[par])./(realmin()+g[par])))
            println("absolute difference between analytic and empirical gradient=$diff")
            println("relative difference between analytic and empirical gradient=$reldiff")

                #if showgrad
                #println("analytic gradient:"); println(g[par])
                #println("empirical gradient:"); println(gemp[par])
                #end

            if diff>tol
                if showgrad
                println("analytic gradient:"); println(g[par])
                println("empirical gradient:"); println(gemp[par])
                end
                print_with_color(:red,"failed: analytic and empiricial gradient mismatch more than $tol\n")
            else
                print_with_color(:green,"passed\n")
            end
        end
    end
    if returngrad
        return gemp
    end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff