text:
################################################################################
#
#  Integer functions
#
################################################################################

function rem(a::fmpz, b::UInt)
  return ccall((:fmpz_fdiv_ui, libflint), UInt, (Ref{fmpz}, UInt), a, b)
end

function isless(a::BigFloat, b::Nemo.fmpz)
  if _fmpz_is_small(b)
    c = ccall((:mpfr_cmp_si, :libmpfr), Int32, (Ref{BigFloat}, Int), a, b.d)
  else
    c = ccall((:mpfr_cmp_z, :libmpfr), Int32, (Ref{BigFloat}, UInt), a, unsigned(b.d) << 2)
  end
  return c < 0
end

function mulmod(a::UInt, b::UInt, n::UInt, ni::UInt)
  ccall((:n_mulmod2_preinv, libflint), UInt, (UInt, UInt, UInt, UInt), a, b, n, ni)
end

# TODO (CF):
# should be Bernstein'ed: this is slow for large valuations
# returns the maximal v s.th. z mod p^v == 0 and z div p^v
#   also useful if p is not prime....
#
# TODO: what happens to z = 0???

function remove(z::T, p::T) where T <: Integer
  z == 0 && return (0, z)
  v = 0
  @assert p > 1
  while mod(z, p) == 0
    z = div(z, p)
    v += 1
  end
  return (v, z)
end

function remove(z::Rational{T}, p::T) where T <: Integer
  z == 0 && return(0, z)
  v, d = remove(denominator(z), p)
  w, n = remove(numerator(z), p)
  return w-v, n//d
end

function valuation(z::T, p::T) where T <: Integer
  z == 0 && return 0
  v = 0
  @assert p > 1
  while mod(z, p) == 0
    z = div(z, p)
    v += 1
  end
  return v
end

function valuation(z::Rational{T}, p::T) where T <: Integer
  z == 0 && error("Not yet implemented")
  v = valuation(denominator(z), p)
  w = valuation(numerator(z), p)
  return w-v
end

function remove!(a::fmpz, b::fmpz)
  v = ccall((:fmpz_remove, libflint), Clong, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), a, a, b)
  return v, a
end

function remove!(a::fmpq, b::fmpz)
  nr = ccall((:fmpq_numerator_ptr, libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  #fmpq's are simplified: either num OR den will be non-trivial
  if vn != 0
    return vn, a
  end
  nr = ccall((:fmpq_denominator_ptr, libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  return -vn, a
end

function valuation!(a::fmpq, b::fmpz)
  nr = ccall((:fmpq_numerator_ptr, libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  #fmpq's are simplified: either num OR den will be non-trivial
  if vn != 0
    return vn
  end
  nr = ccall((:fmpq_denominator_ptr, libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  return -vn
end

function *(a::fmpz, b::BigFloat)
  return BigInt(a)*b
end

function BigFloat(a::fmpq)
  r = BigFloat(0)
  ccall((:fmpq_get_mpfr, libflint), Nothing, (Ref{BigFloat}, Ref{fmpq}, Int32), r, a, __get_rounding_mode())
  return r
end

function isless(a::Float64, b::fmpq) return a<b*1.0; end
function isless(a::fmpq, b::Float64) return a*1.0<b; end

function isless(a::Float64, b::fmpz) return a<b*1.0; end
function isless(a::fmpz, b::Float64) return a*1.0<b; end

function (a::FlintIntegerRing)(b::fmpq)
  !isone(denominator(b)) && error("Denominator not 1")
  return deepcopy(numerator(b))
end

iscommutative(::FlintIntegerRing) = true

#function ^(a::fmpz, k::fmpz)
#  if a == 0
#    if k == 0
#      return fmpz(1)
#    end
#    return fmpz(0)
#  end
#
#  if a == 1
#    return fmpz(1)
#  end
#  if a == -1
#    if isodd(k)
#      return fmpz(-1)
#    else
#      return fmpz(1)
#    end
#  end
#  return a^Int(k)
#end

function ^(a::fmpq, k::fmpz)
  if a == 0
    if k == 0
      return fmpq(1)
    end
    return fmpq(0)
  end

  if a == 1
    return fmpq(1)
  end
  if a == -1
    if isodd(k)
      return fmpq(-1)
    else
      return fmpq(1)
    end
  end
  return a^Int(k)
end

function //(a::fmpz, b::fmpq)
  return fmpq(a)//b
end

function *(a::fmpz, b::Float64)
  return BigInt(a)*b
end

function *(b::Float64, a::fmpz)
  return BigInt(a)*b
end

function *(a::fmpq, b::Float64)
  return Rational(a)*b
end

function *(b::Float64, a::fmpq)
  return Rational(a)*b
end

function Float64(a::fmpq)
  b = a*fmpz(2)^53
  Float64(div(numerator(b), denominator(b)))/(Float64(2)^53) #CF 2^53 is bad in 32bit
end

function convert(R::Type{Rational{Base.GMP.BigInt}}, a::Nemo.fmpz)
  return R(BigInt(a))
end

log(a::fmpz) = log(BigInt(a))
log(a::fmpq) = log(numerator(a)) - log(denominator(a))

function one(::Type{Nemo.fmpz})
  return fmpz(1)
end

function divisible(x::Integer, y::Integer)
  return iszero(rem(x, y))
end

@doc Markdown.doc"""
    modord(a::fmpz, m::fmpz) -> Int
    modord(a::Integer, m::Integer)

  The multiplicative order of a modulo $m$ (not a good algorithm).
"""
function modord(a::fmpz, m::fmpz)
  gcd(a,m)!=1 && throw("1st agrument not a unit")
  i = 1
  b = a % m
  while b != 1
    i += 1
    b = b*a % m
  end
  return i
end

function modord(a::Integer, m::Integer)
  gcd(a,m)!=1 && throw("1st agrument not a unit")
  i = 1
  b = a % m
  while b != 1
    i += 1
    b = b*a % m
  end
  return i
end

if Nemo.version() <= v"0.15.1"
  function isodd(a::fmpz)
    ccall((:fmpz_is_odd, libflint), Int, (Ref{fmpz},), a) == 1
  end

  function iseven(a::fmpz)
    ccall((:fmpz_is_even, libflint), Int, (Ref{fmpz},), a) == 1
  end
end

function neg!(a::fmpz)
  ccall((:fmpz_neg, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}), a, a)
  return a
end

##
## to support rand(fmpz:fmpz)....
##
# note, we cannot get a UnitRange as this is only legal for subtypes of Real

function colon(a::Int, b::fmpz)
  return fmpz(a):b
end

function one(::fmpz)
  return fmpz(1)
end

function zero(::fmpz)
  return fmpz(0)
end

#the build-in show fails due to Integer(a::fmpz) not defined
# I don't think it would efficient to provide this for printing
#display() seems to NOT call my show function, why, I don't know
function Integer(a::fmpz)
  return BigInt(a)
end

function show(io::IO, a::StepRange{fmpz, fmpz})
  println(io, "2-element ", typeof(a), ":\n ", a.start, ",", a.stop)
end

#TODO
# need to be mapped onto proper Flint primitives
# flints needs a proper interface to randomness - I think
# currently one simply cannot use it at all
#
# should be tied(?) to the Julia rng stuff?
# similarly, all the derived rand functions should probably also do this
#
# inspired by/ copied from the Base/random.jl
#

function rand(a::StepRange{fmpz, fmpz})
  return rand(Random.GLOBAL_RNG, a)
end

function rand(rng::AbstractRNG, a::StepRange{fmpz, fmpz})
  m = Base.last(a) - Base.first(a)
  m < 0 && throw("range empty")
  nd = ndigits(m, 2)
  nl, high = divrem(nd, 8*sizeof(Base.GMP.Limb))
  if high>0
    mask = m>>(nl*8*sizeof(Base.GMP.Limb))
  end
  s = fmpz(0)
  while true
    s = fmpz(0)
    for i=1:nl
      s = s << (8*sizeof(Base.GMP.Limb))
      s += rand(rng, Base.GMP.Limb)
    end
    if high >0
      s = s << high
      s += rand(rng, 0:Base.GMP.Limb(mask))
    end
    if s <= m break; end
  end
  return s + first(a)
end

function length(a::StepRange{fmpz, fmpz})
  return a.stop - a.start +1
end

struct RangeGeneratorfmpz# <: Base.Random.RangeGenerator
  a::StepRange{fmpz, fmpz}
end

function Random.RangeGenerator(r::StepRange{fmpz,fmpz})
    m = last(r) - first(r)
    m < 0 && throw(ArgumentError("range must be non-empty"))
    return RangeGeneratorfmpz(r)
end

function rand(rng::AbstractRNG, g::RangeGeneratorfmpz)
  return rand(rng, g.a)
end

function Base.getindex(a::StepRange{fmpz,fmpz}, i::fmpz)
  a.start+(i-1)*Base.step(a)
end

function Base.divrem(a::fmpz, b::Int)
  return (div(a, b), rem(a, b))
end

################################################################################
#
#  Should go to Nemo?
#
################################################################################

one(::Type{fmpq}) = fmpq(1)

############################################################
# more unsafe function that Bill does not want to have....
############################################################

function mod!(z::fmpz, x::fmpz, y::fmpz)
  ccall((:fmpz_mod, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
  return z
end

function divexact!(z::fmpz, x::fmpz, y::fmpz)
  iszero(y) && throw(DivideError())
  ccall((:fmpz_divexact, libflint), Nothing,
        (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
  return z
end

function lcm!(z::fmpz, x::fmpz, y::fmpz)
   ccall((:fmpz_lcm, libflint), Nothing,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function gcd!(z::fmpz, x::fmpz, y::fmpz)
   ccall((:fmpz_gcd, libflint), Nothing,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function mul!(z::fmpz, x::fmpz, y::Int)
  ccall((:fmpz_mul_si, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Int), z, x, y)
  return z
end

function mul!(z::fmpz, x::fmpz, y::UInt)
  ccall((:fmpz_mul_ui, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, UInt), z, x, y)
  return z
end

function mul!(z::fmpz, x::fmpz, y::Integer)
  mul!(z, x, fmpz(y))
  return z
end

function one!(a::fmpz)
  ccall((:fmpz_one, libflint), Nothing, (Ref{fmpz}, ), a)
  return a
end

################################################################################
#
#  power detection
#
################################################################################
#compare to Oscar/examples/PerfectPowers.jl which is, for large input,
#far superiour over gmp/ fmpz_is_perfect_power

@doc Markdown.doc"""
    ispower(a::fmpz) -> Int, fmpz
    ispower(a::Integer) -> Int, Integer

Returns $e$, $r$ such that $a = r^e$ with $e$ maximal. Note: $1 = 1^0$.
"""
function ispower(a::fmpz)
  if iszero(a)
    error("must not be zero")
  end
  if isone(a)
    return 0, a
  end
  if a < 0
    e, r = ispower(-a)
    if isone(e)
      return 1, a
    end
    v, s = remove(e, 2)
    return s, -r^(2^v)
  end
  rt = fmpz()
  e = 1
  while true
    ex = ccall((:fmpz_is_perfect_power, libflint), Int, (Ref{fmpz}, Ref{fmpz}), rt, a)
    if ex == 1 || ex == 0
      return e, a
    end
    e *= ex
    a = rt
  end
end

function ispower(a::Integer)
  e,r = ispower(fmpz(a))
  return e, typeof(a)(r)
end

@doc Markdown.doc"""
    ispower(a::fmpq) -> Int, fmpq
    ispower(a::Rational) -> Int, Rational

Writes $a = r^e$ with $e$ maximal. Note: $1 = 1^0$.
"""
function ispower(a::fmpq)
  e, r = ispower(numerator(a))
  if e==1
    return e, a
  end
  f, s = ispower(denominator(a))
  g = gcd(e, f)
  return g, r^div(e, g)//s^div(f, g)
end

function ispower(a::Rational)
  T = typeof(denominator(a))
  e, r = ispower(fmpq(a))
  return e, T(numerator(r))//T(denominator(r))
end

@doc Markdown.doc"""
    ispower(a::fmpz, n::Int) -> Bool, fmpz
    ispower(a::fmpq, n::Int) -> Bool, fmpq
    ispower(a::Integer, n::Int) -> Bool, Integer

Tests if $a$ is an $n$-th power. Return `true` and the root if successful.
"""
function ispower(a::fmpz, n::Int)
   if a < 0 && iseven(n)
    return false, a
  end
  b = root(a, n)
  return b^n==a, b
end

function ispower(a::fmpq, n::Int)
  fl, nu = ispower(numerator(a), n)
  if !fl
    return fl, a
  end
  fl, de = ispower(denominator(a), n)
  return fl, fmpq(nu, de)
end

function issquare(x::fmpq)
  return x > 0 && issquare(numerator(x)) && issquare(denominator(x))
end

function issquare_with_sqrt(x::fmpq)
  if x < 0
    return false, x
  else
    fl, n = issquare_with_sqrt(numerator(x))
    !fl && return false, x
    fl, d = issquare_with_sqrt(denominator(x))
    !fl && return false, x
    return true, fmpq(n, d)
  end
end

################################################################################
#
#  Chinese remaindering modulo UInts to fmpz
#
################################################################################

mutable struct fmpz_comb
  primes::Ptr{UInt}
  num_primes::Int
  n::Int
  comb::Ptr{Ptr{fmpz}}
  res::Ptr{Ptr{fmpz}}
  mod_n::UInt
  mod_ninv::UInt
  mod_norm::UInt

  function fmpz_comb(primes::Vector{UInt})
    z = new()
    ccall((:fmpz_comb_init, libflint), Nothing, (Ref{fmpz_comb}, Ptr{UInt}, Int),
            z, primes, length(primes))
    finalizer(z, _fmpz_comb_clear_fn)
    return z
  end
end

function _fmpz_comb_clear_fn(z::fmpz_comb)
  ccall((:fmpz_comb_clear, libflint), Nothing, (Ref{fmpz_comb}, ), z)
end

mutable struct fmpz_comb_temp
  n::Int
  comb_temp::Ptr{Ptr{fmpz}}
  temp::Ptr{fmpz}
  temp2::Ptr{fmpz}

  function fmpz_comb_temp(comb::fmpz_comb)
    z = new()
    ccall((:fmpz_comb_temp_init, libflint), Nothing,
            (Ref{fmpz_comb_temp}, Ref{fmpz_comb}), z, comb)
    finalizer(z, _fmpz_comb_temp_clear_fn)
    return z
  end
end

function _fmpz_comb_temp_clear_fn(z::fmpz_comb_temp)
  ccall((:fmpz_comb_temp_clear, libflint), Nothing, (Ref{fmpz_comb_temp}, ), z)
end

function fmpz_multi_crt_ui!(z::fmpz, a::Vector{UInt}, b::fmpz_comb, c::fmpz_comb_temp)
  ccall((:fmpz_multi_CRT_ui, libflint), Nothing,
          (Ref{fmpz}, Ptr{UInt}, Ref{fmpz_comb}, Ref{fmpz_comb_temp}, Cint),
          z, a, b, c, 1)
  return z
end

function _fmpz_preinvn_struct_clear_fn(z::fmpz_preinvn_struct)
  ccall((:fmpz_preinvn_clear, libflint), Nothing, (Ref{fmpz_preinvn_struct}, ), z)
end

function fdiv_qr_with_preinvn!(q::fmpz, r::fmpz, g::fmpz, h::fmpz, hinv::fmpz_preinvn_struct)
  ccall((:fmpz_fdiv_qr_preinvn, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz_preinvn_struct}), q, r, g, h, hinv)
end

function submul!(z::fmpz, x::fmpz, y::fmpz)
  ccall((:fmpz_submul, libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
end

################################################################################
#
#  Number of bits
#
################################################################################

function nbits(a::BigInt)
  return ndigits(a, 2)
end

@doc Markdown.doc"""
    nbits(a::Int) -> Int
    nbits(a::UInt) -> Int
    nbits(a::BigInt) -> Int

  Returns the number of bits necessary to represent $a$.
"""
function nbits(a::Int)
  a==0 && return 0
  return floor(Int, (log(abs(a))/log(2)))+1
end

function nbits(a::UInt)
  a==0 && return 0
  return floor(Int, (log(a)/log(2)))+1
end

################################################################################
#
#  Modular reduction with symmetric residue system
#
################################################################################

function mod_sym(a::fmpz, b::fmpz)
  c = mod(a,b)
  @assert c>=0
  if b > 0 && 2*c > b
    return c-b
  elseif b < 0 && 2*c > -b
    return c+b
  else
    return c
  end
end

@doc Markdown.doc"""
    isinteger(a::fmpq) -> Bool

Returns `true` iff the denominator of $a$ is one.
"""
function isinteger(a::fmpq)
  return isone(denominator(a))
end

################################################################################
#
#  sunit group
#
################################################################################

mutable struct MapSUnitGrpZFacElem <: Map{GrpAbFinGen, FacElemMon{FlintRationalField}, HeckeMap, MapSUnitGrpZFacElem}
  header::MapHeader{GrpAbFinGen, FacElemMon{FlintRationalField}}
  idl::Vector{fmpz}

  function MapSUnitGrpZFacElem()
    return new()
  end
end

function show(io::IO, mC::MapSUnitGrpZFacElem)
  println(io, "SUnits (in factored form) map of $(codomain(mC)) for $(mC.idl)")
end

mutable struct MapSUnitGrpZ <: Map{GrpAbFinGen, FlintRationalField, HeckeMap, MapSUnitGrpZ}
  header::MapHeader{GrpAbFinGen, FlintRationalField}
  idl::Vector{fmpz}

  function MapSUnitGrpZ()
    return new()
  end
end

function show(io::IO, mC::MapSUnitGrpZ)
  println(io, "SUnits map of $(codomain(mC)) for $(mC.idl)")
end

@doc Markdown.doc"""
    sunit_group_fac_elem(S::Vector{fmpz}) -> GrpAbFinGen, Map
    sunit_group_fac_elem(S::Vector{Integer}) -> GrpAbFinGen, Map

The $S$-unit group of $Z$ supported at $S$: the group of
rational numbers divisible only by primes in $S$.
The second return value is the map mapping group elements to rationals
in factored form or rationals back to group elements.
"""
function sunit_group_fac_elem(S::Vector{T}) where T <: Integer
  return sunit_group_fac_elem(fmpz[x for x=S])
end

function sunit_group_fac_elem(S::Vector{fmpz})
  S = coprime_base(S)  #TODO: for S-units use factor???
  G = abelian_group(vcat(fmpz[2], fmpz[0 for i=S]))
  S = vcat(fmpz[-1], S)

  mp = MapSUnitGrpZFacElem()
  mp.idl = S

  Sq = fmpq[x for x=S]

  function dexp(a::GrpAbFinGenElem)
    return FacElem(Sq, fmpz[a.coeff[1,i] for i=1:length(S)])
  end

  mp.header = MapHeader(G, FacElemMon(FlintQQ), dexp)

  return G, mp
end

function preimage(f::MapSUnitGrpZFacElem, a::fmpz)
  g = Int[a>=0 ? 0 : 1]
  S = f.idl
  g = vcat(g, Int[valuation(a, x) for x=S[2:end]])
  return domain(f)(g)
end

function preimage(f::MapSUnitGrpZFacElem, a::Integer)
  return preimage(f, fmpz(a))
end

function preimage(f::MapSUnitGrpZFacElem, a::Rational)
  return preimage(f, fmpq(a))
end

function preimage(f::MapSUnitGrpZFacElem, a::fmpq)
  return preimage(f, numerator(a)) - preimage(f, denominator(a))
end

function preimage(f::MapSUnitGrpZFacElem, a::FacElem)
  return sum(GrpAbFinGenElem[e*preimage(f, k) for (k,e) = a.fac])
end

@doc Markdown.doc"""
    sunit_group(S::Vector{fmpz}) -> GrpAbFinGen, Map
    sunit_group(S::Vector{Integer}) -> GrpAbFinGen, Map

The $S$-unit group of $Z$ supported at $S$: the group of
rational numbers divisible only by primes in $S$.
The second return value is the map mapping group elements to rationals
or rationals back to group elements.
"""
function sunit_group(S::Vector{T}) where T <: Integer
  return sunit_group(fmpz[x for x=S])
end

function sunit_group(S::Vector{fmpz})
  u, mu = sunit_group_fac_elem(S)

  mp = MapSUnitGrpZ()
  mp.idl = S

  function dexp(a::GrpAbFinGenElem)
    return evaluate(image(mu, a))
  end

  mp.header = MapHeader(u, FlintQQ, dexp, mu.header.preimage)

  return u, mp
end

@doc Markdown.doc"""
    isprime_power(n::fmpz) -> Bool
    isprime_power(n::Integer) -> Bool

Tests if $n$ is the exact power of a prime number.
"""
function isprime_power(n::fmpz)
  e, p = ispower(n)
  return isprime(p)
end

function isprime_power(n::Integer)
  return isprime_power(fmpz(n))
end

################################################################################
# random and factor
################################################################################

factor(a...; b...) = Nemo.factor(a...; b...)

factor(a::Integer) = factor(fmpz(a))

mutable struct flint_rand_ctx_t
  a::Ptr{Nothing}
  function flint_rand_ctx_t()
    return new()
  end
end

function show(io::IO, A::flint_rand_ctx_t)
  println(io, "Flint random state")
end

function flint_rand_state()
  A = flint_rand_ctx_t()
  A.a = ccall((:flint_rand_alloc, libflint), Ptr{Nothing}, (Int, ), 1)
  ccall((:flint_randinit, libflint), Nothing, (Ptr{Nothing}, ), A.a)

  function clean_rand_state(A::flint_rand_ctx_t)
    ccall((:flint_randclear, libflint), Nothing, (Ptr{Nothing}, ), A.a)
    ccall((:flint_rand_free, libflint), Nothing, (Ptr{Nothing}, ), A.a)
    nothing
  end
  finalizer(clean_rand_state, A)
  return A
end

global flint_rand_ctx

function ecm(a::fmpz, B1::UInt, B2::UInt, ncrv::UInt, rnd = flint_rand_ctx)
  f = fmpz()
  r = ccall((:fmpz_factor_ecm, libflint), Int32, (Ref{fmpz}, UInt, UInt, UInt, Ptr{Nothing}, Ref{fmpz}), f, ncrv, B1, B2, rnd.a, a)
  return r, f
end

function ecm(a::fmpz, B1::Int, B2::Int, ncrv::Int, rnd = flint_rand_ctx)
  return ecm(a, UInt(B1), UInt(B2), UInt(ncrv), rnd)
end

#data from http://www.mersennewiki.org/index.php/Elliptic_Curve_Method
B1 = [2, 11, 50, 250, 1000, 3000, 11000, 43000, 110000, 260000, 850000, 2900000];
nC = [25, 90, 300, 700, 1800, 5100, 10600, 19300, 49000, 124000, 210000, 340000];

function ecm(a::fmpz, max_digits::Int = div(ndigits(a), 3), rnd = flint_rand_ctx)
  n = ndigits(a, 10)
  B1s = 15

  i = 1
  s = max(div(max_digits - 10, 5), 1)
  #i = s = max(i, s)
  while i <= s
    e, f = ecm(a, B1[i]*1000, B1[i]*1000*100, nC[i], rnd)
    if e != 0
      return (e,f)
    end
    i += 1
    if i > length(B1)
      return (e, f)
    end
  end
  return (Int32(0), a)
end

function factor_trial_range(N::fmpz, start::Int=0, np::Int=10^5)
   F = Nemo.fmpz_factor()
   ccall((:fmpz_factor_trial_range, libflint), Nothing, (Ref{Nemo.fmpz_factor}, Ref{fmpz}, UInt, UInt), F, N, start, np)
   res = Dict{fmpz, Int}()
   for i in 1:F.num
     z = fmpz()
     ccall((:fmpz_factor_get_fmpz, libflint), Nothing,
           (Ref{fmpz}, Ref{Nemo.fmpz_factor}, Int), z, F, i - 1)
     res[z] = unsafe_load(F.exp, i)
   end
   return res, canonical_unit(N)
end

const big_primes = fmpz[]

function factor(N::fmpz)
  if iszero(N)
    throw(ArgumentError("Argument is not non-zero"))
  end
  N_in = N
  global big_primes
  r, c = factor_trial_range(N)
  for (p, v) = r
    N = divexact(N, p^v)
  end
  if isunit(N)
    @assert N == c
    return Nemo.Fac(c, r)
  end
  N *= c
  @assert N > 0

  for p = big_primes
    v, N = remove(N, p)
    if v > 0
      @assert !haskey(r, p)
      r[p] = v
    end
  end
  factor_insert!(r, N)
  for p = keys(r)
    if nbits(p) > 60 && !(p in big_primes)
      push!(big_primes, p)
    end
  end
  return Nemo.Fac(c, r)
end

function factor_insert!(r::Dict{fmpz, Int}, N::fmpz, scale::Int = 1)
  #assumes N to be positive
  #        no small divisors
  #        no big_primes
  if isone(N)
    return r
  end
  fac, N = ispower(N)
  if fac > 1
    return factor_insert!(r, N, fac)
  end
  if isprime(N)
    @assert !haskey(r, N)
    r[N] = scale
    return r
  end
  if ndigits(N) < 60
    s = Nemo.factor(N) #MPQS
    for (p, k) in s
      if haskey(r, p)
        r[p] += k*scale
      else
        r[p] = k*scale
      end
    end
    return r
  end

  e, f = ecm(N)
  if e == 0
    s = Nemo.factor(N)
    for (p, k) in s
      if haskey(r, p)
        r[p] += k*scale
      else
        r[p] = k*scale
      end
    end
    return r
  end
  cp = coprime_base([N, f])
  for i = cp
    factor_insert!(r, i, scale*valuation(N, i))
  end
  return r
end

#TODO: problem(s)
# Nemo.factor = mpqs is hopeless if > n digits, but asymptotically and practically
# faster than ecm.
# ecm is much better if there are "small" factors.
# p-1 and p+1 methods are missing
# so probably
# if n is small enough -> Nemo
# if n is too large: ecm
# otherwise
#  need ecm to find small factors
# then recurse...

function _factors_trial_division(n::fmpz, np::Int = 10^5)
  res, u = factor_trial_range(n, 0, np)
  factors = fmpz[]
  for (p, v) in res
    push!(factors, p)
    n = divexact(n, p^v)
  end
  return factors, n

end

function ceil(::Type{fmpz}, a::BigFloat)
  return fmpz(ceil(BigInt, a))
end

function ceil(::Type{Int}, a::fmpq)
  return Int(ceil(fmpz, a))
end

function floor(::Type{fmpz}, a::BigFloat)
  return fmpz(floor(BigInt, a))
end

function floor(::Type{Int}, a::fmpq)
  return Int(floor(fmpz, a))
end

function round(::Type{fmpz}, a::BigFloat)
  return fmpz(round(BigInt, a))
end

function round(::Type{Int}, a::BigFloat)
  return Int(round(fmpz, a))
end

/(a::BigFloat, b::fmpz) = a/BigInt(b)

function rand!(A::Vector{fmpz}, v::StepRange{fmpz, fmpz})
  for i in 1:length(A)
    A[i] = rand(v)
  end
  return A
end

Base.isless(a::Int, b::fmpz) = a < b

Base.isless(a::fmpz, b::Int) = a < b

function (::Type{Base.Rational{BigInt}})(x::fmpq)
  return Rational{BigInt}(BigInt(numerator(x)), BigInt(denominator(x)))
end

export euler_phi_inv, Divisors, carmichael_lambda

@doc Markdown.doc"""
    Divisors{T}

An iterator for the divisors of a given object.
Create using
    `Divisors(A, power::Int = 1)`
where $A$ is either a FacElem or a direct element. Power can be used
to restrict to objects $B$ s.th. $B$^power still divides $A$, e.g.
    `Divisors(12, powers = 2)`
will produce square divisors.

For rings where this makes sense, i.e. where the unit group is finite,
```units = true``` can be passed in to also take into account
the units.
"""
mutable struct Divisors{T}
  n::T
  lf::MSet{T}
  s#::Iterator
  f::Function
  U::GrpAbFinGen
  function Divisors(a::T; units::Bool = false, power::Int = 1) where {T}
    r = new{T}()
    r.n = a
    r.lf = MSet{T}()
    for (p, k) = factor(a).fac
      k = div(k, power)
      if k > 0
        push!(r.lf, p, k)
      end
    end
    r.f = x -> length(x) == 0 ? one(parent(a)) : prod(x)
    r.s = subsets(r.lf)
    if units
      U, mU = unit_group(parent(a))
      r.U = U
      g = r.f
      r.f = x->mU(x[1]) * g(x[2])
      r.s = Base.Iterators.ProductIterator((U, r.s))
    end
   return r
  end
  function Divisors(a::NfOrdIdl; units::Bool = false, power::Int = 1)
    r = new{NfOrdIdl}()
    r.n = a
    r.lf = MSet{NfOrdIdl}()
    for (p, k) = factor(a)
      k = div(k, power)
      if k > 0
        push!(r.lf, p, k)
      end
    end
    r.f = x -> length(x) == 0 ? one(parent(a)) : prod(x)
    r.s = subsets(r.lf)
    return r
  end
  function Divisors(a::FacElem{NfOrdIdl}; units::Bool = false, power::Int = 1)
    r = new{NfOrdIdl}()
    r.n = evaluate(a)
    r.lf = MSet{NfOrdIdl}()
    for (p, k) = factor(a)
      k = div(k, power)
      if k > 0
        push!(r.lf, p, k)
      end
    end
    r.f = x -> length(x) == 0 ? one(parent(a)) : prod(x)
    r.s = subsets(r.lf)
    return r
  end

  function Divisors(a::FacElem{fmpz, FlintIntegerRing}; units::Bool = false, power::Int = 1)
    r = new{fmpz}()
    r.n = evaluate(a)
    r.lf = MSet{fmpz}()
    for (p, k) = factor(a).fac
      k = div(k, power)
      if k > 0
        push!(r.lf, p, k)
      end
    end
    r.f = x -> length(x) == 0 ? one(parent(r.n)) : prod(x)
    r.s = subsets(r.lf)
    if units
      U, mU = unit_group(parent(a))
      r.U = U
      g = r.f
      r.f = x->mU(x[1]) * g(x[2])
      r.s = Base.Iterators.ProductIterator((U, r.s))
    end
   return r
  end
  function Divisors(a::Fac{fmpz}; units::Bool = false, power::Int = 1)
    return Divisors(FacElem(a), units = units, power = power)
  end
end
Base.IteratorSize(::Divisors) = Base.HasLength()
Base.length(D::Divisors) = length(D.s)

function Base.iterate(D::Divisors)
  x = iterate(D.s)
  if x == nothing
    return x
  end
  return D.f(x[1]), x[2]
end

function Base.iterate(D::Divisors, i)
  x = iterate(D.s, i)
  if x == nothing
    return x
  end
  return D.f(x[1]), x[2]
end

function Base.show(io::IO, D::Divisors)
  print(io, "Divisors of $(D.n) = $(D.lf)")
  if isdefined(D, :U)
    print(io, " times $(D.U)")
  end
  print(io, "\n")
end

@doc Markdown.doc"""
    unit_group(::FlintIntegerRing) -> GrpAbFinGen, Map

The unit group of $\mathbb{Z}$, i.e. $C_2$ and the map translating between the group and $\mathbb{Z}$.
"""
function unit_group(::FlintIntegerRing)
  G = abelian_group([2])
  exp = function(z::GrpAbFinGenElem)
    return isodd(z[1]) ? fmpz(-1) : fmpz(1)
  end
  log = function(z::fmpz)
    return z == -1 ? G[1] : G[0]
  end
  return G, MapFromFunc(exp, log, G, FlintZZ)
end

@doc Markdown.doc"""
    unit_group(::Integers{T}) -> GrpAbFinGen, Map

The unit group of , i.e. $C_2$ and the map translating between the group and $\mathbb{Z}$.
"""
function unit_group(R::AbstractAlgebra.Integers{T}) where {T}
  G = abelian_group([2])
  exp = function(z::GrpAbFinGenElem)
    return isodd(z[1]) ? T(-1) : T(1)
  end
  log = function(z::T)
    return z == -1 ? G[1] : G[0]
  end
  return G, MapFromFunc(exp, log, G, R)
end

#Nemo.GaloisField = nmod?
# PolyRing

#basically from
#http://people.math.gatech.edu/~ecroot/shparlinski_final.pdf
#Contini, Croot, Shparlinski: Complexity of inverting the Euler function
@doc Markdown.doc"""
    euler_phi_inv_fac_elem(n::fmpz)

The inverse of the Euler totient functions: find all $x$ s.th. $phi(x) = n$
holds. The elements are returned in factored form.
"""
function euler_phi_inv_fac_elem(n::fmpz)
  lp = fmpz[]
  for d = Divisors(n)
    if isprime(d+1)
      push!(lp, d+1)
    end
  end
#  println("possible primes: ", lp)

  E = Tuple{fmpz, Vector{Tuple{fmpz, Int}}}[]
  res = FacElem{fmpz, FlintIntegerRing}[]
  for p = lp
    v = valuation(n, p)
    for i=0:v
      push!(E, ((p-1)*p^i, [(p, i+1)]))
      if E[end][1] == n
        push!(res, FacElem(Dict(E[end][2])))
      end
    end
  end
  while true
    F = []
    for e = E
      nn = divexact(n, e[1])
      x = e[2]
      pm = x[end][1]
      for p = lp
        if p <= pm
          continue
        end
        if nn % (p-1) == 0
          v = valuation(nn, p)
          for i = 0:v
            push!(F, (e[1]*(p-1)*p^i, vcat(e[2], [(p, i+1)])))
            if F[end][1] == n
              push!(res, FacElem(Dict(F[end][2])))
            end
          end
        end
      end
    end
    if length(F) == 0
      return res
    end
    E = F
  end
end

function euler_phi(x::Fac{fmpz})
  return prod((p-1)*p^(v-1) for (p,v) = x.fac)
end

function euler_phi(x::FacElem{fmpz, FlintIntegerRing})
  x = factor(x)
  return prod((p-1)*p^(v-1) for (p,v) = x.fac)
end

function carmichael_lambda(x::Fac{fmpz})
  if haskey(x.fac, fmpz(2))
    y = deepcopy(x.fac)
    v = y[fmpz(2)]
    delete!(y, fmpz(2))
    if v > 2
      c = fmpz(2)^(v-2)
    else
      c = fmpz(1)
    end
  else
    c = fmpz(1)
    y = x.fac
  end
  if length(y) == 0
    return c
  end
  return c * reduce(lcm, (p-1)*p^(v-1) for (p,v) = y)
end

function carmichael_lambda(x::fmpz)
  v, x = remove(x, fmpz(2))
  if isone(x)
    c = x
  else
    x = factor(x)
    c = reduce(lcm, (p-1)*p^(v-1) for (p,v) = x.fac)
  end
  if v < 2
    return c
  else
    return fmpz(2)^(v-2)*c
  end
end

function carmichael_lambda(x::FacElem{fmpz, FlintIntegerRing})
  x = factor(x)
  return carmichael_lambda(x)
end

function carmichael_lambda(n::T) where {T <: Integer}
  return T(carmichael_lambda(fmpz(n)))
end

@doc Markdown.doc"""
    euler_phi_inv(n::Integer) -> Vector{fmpz}

The inverse of the Euler totient functions: find all $x$ s.th. $phi(x) = n$
holds.
"""
function euler_phi_inv(n::Integer)
  return euler_phi_inv(fmpz(n))
end

@doc Markdown.doc"""
    euler_phi_inv(n::fmpz) -> Vector{fmpz}

The inverse of the Euler totient functions: find all $x$ s.th. $phi(x) = n$
holds.
"""
function euler_phi_inv(n::fmpz)
  return [ evaluate(x) for x = euler_phi_inv_fac_elem(n)]
end

function factor(a::FacElem{fmpz, FlintIntegerRing})
  b = simplify(a)
  c = Dict{fmpz, Int}()
  s = fmpz(1)
  for (p,k) = b.fac
    lp = factor(p)
    s *= lp.unit
    for (q,w) = lp.fac
      c[q] = w*k
    end
  end
  l = Fac{fmpz}()
  l.fac = c
  l.unit = s
  return l
end

function FacElem(a::Fac{fmpz})
  f = FacElem(a.fac)
  if a.unit == -1
    return a.unit * f
  end
  return f
end

#= for torsion units:

   [maximum([maximum(vcat([fmpz(-1)], euler_phi_inv(x))) for x = Divisors(fmpz(n))]) for n = 1:250]

=#

radical(a::fmpz) = prod(keys(factor(a).fac))
function radical(a::T) where {T <: Integer}
  return T(radical(fmpz(a)))
end

function quo(::FlintIntegerRing, a::fmpz)
  R = ResidueRing(FlintZZ, a)
  f = MapFromFunc(x -> R(x), y->lift(y), FlintZZ, R)
  return R, f
end

function quo(::FlintIntegerRing, a::Integer)
  R = ResidueRing(FlintZZ, a)
  f = MapFromFunc(x -> R(x), y->lift(y), FlintZZ, R)
  return R, f
end

function (::FlintIntegerRing)(x::Rational{Int})
  @assert denominator(x) == 1
  return fmpz(numerator(x))
end

module BitsMod

using Hecke
import Nemo
import Base: ^, show, getindex, iterate, length
import Hecke.bits
export bits, Limbs

const hb = UInt(1)<<63

#= not used - lacks length
struct BitsUInt
  a::UInt
end

function bits(a::UInt)
  l = nbits(a)
  return BitsUInt(a<<(sizeof(a)*8-l))
end

function Base.iterate(x::BitsUInt)
  return iterate(x, x.a)
end

@inline function Base.iterate(x::BitsUInt, u::UInt)
  iszero(u) && return nothing
  return (u&hb) != 0, u<<1
end
=#

struct Limbs
  a::fmpz
  len::Int
  b::Ptr{UInt}
  function Limbs(a::fmpz; MSW::Bool = true)
    if Nemo._fmpz_is_small(a)
      return new(a, 0, convert(Ptr{UInt}, 0))
    end
    z = convert(Ptr{Cint}, unsigned(a.d)<<2)
    len = unsafe_load(z, 2)
    d = convert(Ptr{Ptr{UInt}}, unsigned(a.d) << 2) + 2*sizeof(Cint)
    p = unsafe_load(d)
    if !MSW
      new(a, -len, p)
    else
      new(a, len, p)
    end
  end
end

function show(io::IO, L::Limbs)
  print(io, "limb-access for: ", L.a)
end

@inline function getindex(L::Limbs, i::Int)
  if L.len == 0
    return UInt(abs(L.a.d)) #error???
  end
  @boundscheck @assert i <= abs(L.len)
  return unsafe_load(L.b, i)
end

function iterate(L::Limbs)
  L.len < 0 && return L[1], 1

  return L[L.len], L.len
end

function iterate(L::Limbs, i::Int)
  if L.len < 0
    i > -L.len && return nothing
    return L[i+1], i+1
  end
  i == 0 && return nothing
  return L[i-1], i-1
end

length(L::Limbs) = L.len+1

#=
#from https://github.com/JuliaLang/julia/issues/11592
#compiles directly down to the ror/rol in assembly
for T in Base.BitInteger_types
  mask = UInt8(sizeof(T) << 3 - 1)
  @eval begin
    ror(x::$T, k::Integer) = (x >>> ($mask & k)) | (x <<  ($mask & -k))
    rol(x::$T, k::Integer) = (x <<  ($mask & k)) | (x >>> ($mask & -k))
  end
end
=#

struct BitsFmpz
  L::Limbs

  function BitsFmpz(b::fmpz)
    return new(Limbs(b))
  end
end

function iterate(B::BitsFmpz)
  L = B.L
  a = L[L.len]
  b = UInt(1) << (nbits(a)-1)
  return true, (b, L.len)
end

@inline function iterate(B::BitsFmpz, s::Tuple{UInt, Int})
  b = s[1] >> 1
  if b == 0
    l = s[2] - 1
    if l < 1
      return nothing
    end
    b = hb
    a = B.L[l]
    return a&b != 0, (b, l)
  end
  return B.L[s[2]]&b != 0, (b, s[2])
end

function show(io::IO, B::BitsFmpz)
  print(io, "bit iterator for:", B.L.a)
end

length(B::BitsFmpz) = nbits(B.L.a)

bits(a::fmpz) = BitsFmpz(a)
#= wrong order, thus disabled

function getindex(B::BitsFmpz, i::Int)
  return ccall((:fmpz_tstbit, libflint), Int, (Ref{fmpz}, Int), B.L.a, i) != 0
end
=#

end

using .BitsMod
export bits, Limbs

^(a::T, n::Union{Integer, fmpz}) where {T <: RingElem} = _generic_power(a, n)

^(a::NfAbsOrdIdl, n::Union{Integer, fmpz})  = _generic_power(a, n)

#^(a::NfRelOrdIdl, n::Union{Integer, fmpz})  = _generic_power(a, n)

function _generic_power(a, n::Union{Integer, fmpz})
  fits(Int, n) && return a^Int(n)
  if isnegative(n)
    a = inv(a)
    n = -n
  end
  r = one(parent(a))
  for b = bits(n)
    r = mul!(r, r, r)
    if b
      r = mul!(r, r, a)
    end
  end
  return r
end

#square-and-multiply algorithm to compute f^e mod g
function powermod(f::T, e::fmpz, g::T) where {T}
    #small exponent -> use powermod
    if nbits(e) <= 63
        return powermod(f, Int(e), g)
    else
        #go through binary representation of exponent and multiply with res
        #or (res and f)
        res = parent(f)(1)
        for b=bits(e)
            res = mod(res^2, g)
            if b
                res = mod(res*f, g)
            end
        end
        return res
    end
end

################################################################################
#
#  Prime numbers up to
#
################################################################################

@doc Markdown.doc"""
    primes_up_to(n::Int) -> Vector{Int}

Returns a vector containing all the prime numbers up to $n$.
"""
function primes_up_to(n::Int)
  list = trues(n)
  list[1] = false
  i = 2
  s = 4
  while s <= n
    list[s] = false
    s += 2
  end
  i = 3
  b = sqrt(n)
  while i <= b
    if list[i]
      j = 3*i
      s = 2*i
      while j <= n
        list[j] = false
        j += s
      end
    end
    i += 1
  end
  return findall(list)
end

################################################################################
#
#  Squarefree numbers up to
#
################################################################################

@doc Markdown.doc"""
    squarefree_up_to(n::Int) -> Vector{Int}

Returns a vector containing all the squarefree numbers up to $n$.
"""
function squarefree_up_to(n::Int; coprime_to::Vector{fmpz} = fmpz[], prime_base::Vector{fmpz} = fmpz[])

  @assert isempty(coprime_to) || isempty(prime_base)
  if !isempty(prime_base)
    listi = Int[1]
    for x in prime_base
      if x > n
        continue
      end
      el = Int(x)
      to_add = Int[]
      for i = 1:length(listi)
        eln = el * listi[i]
        if eln <= n
          push!(to_add, eln)
        else
          break
        end
      end
      append!(listi, to_add)
      sort!(listi)
    end
    return listi
  end
  list = trues(n)
  for x in coprime_to
    if x > n
      continue
    end
    t = Int(x)
    while t <= n
      list[t] = false
      t += Int(x)
    end
  end
  i = 2
  b = isqrt(n)
  lp = primes_up_to(b)
  for i = 1:length(lp)
    p2 = lp[i]^2
    ind = p2
    while ind <= n
      list[ind] = false
      ind += p2
    end
  end
  return findall(list)
end

################################################################################
#
#  issquarefree
#
################################################################################

#TODO (Hard): Implement this properly.
@doc Markdown.doc"""
    issquarefree(n::Union{Int, fmpz}) -> Bool

Returns true if $n$ is squarefree, false otherwise.
"""
function issquarefree(n::Union{Int,fmpz})
  if iszero(n)
    error("Argument must be non-zero")
  end
  if isone(abs(n))
    return true
  end
  e, b = ispower(n)
  if e > 1
    return false
  end
  return isone(maximum(values(factor(n).fac)))
end

################################################################################
#
#  Rounding and friends
#
################################################################################

Base.floor(::Type{fmpz}, x::fmpz) = x

Base.ceil(::Type{fmpz}, x::fmpz) = x

Base.floor(::Type{fmpz}, x::Int) = fmpz(x)

Base.ceil(::Type{fmpz}, x::Int) = fmpz(x)

Base.floor(::Type{fmpz}, x::fmpq) = fdiv(numerator(x), denominator(x))

Base.ceil(::Type{fmpz}, x::fmpq) = cdiv(numerator(x), denominator(x))

Base.round(x::fmpq, ::RoundingMode{:Up}) = ceil(x)

Base.round(::Type{fmpz}, x::fmpq, ::RoundingMode{:Up}) = ceil(fmpz, x)

Base.round(x::fmpq, ::RoundingMode{:Down}) = floor(x)

Base.round(::Type{fmpz}, x::fmpq, ::RoundingMode{:Down}) = floor(fmpz, x)

function Base.round(x::fmpq, ::RoundingMode{:Nearest})
  d = denominator(x)
  n = numerator(x)
  if d == 2
    if mod(n, 4) == 1
      if n > 0
        return div(n, d)
      else
        return div(n, d) - 1
      end
    else
      if n > 0
        return div(n, d) + 1
      else
        return div(n, d)
      end
    end
  end

  return floor(x + 1//2)
end

Base.round(x::fmpq, ::RoundingMode{:NearestTiesAway}) = sign(x) * floor(abs(x) + 1//2)

Base.round(::Type{fmpz}, x::fmpq, ::RoundingMode{:NearestTiesAway}) = sign(x) == 1 ? floor(fmpz, abs(x) + 1//2) : -floor(fmpz, abs(x) + 1//2)

function Base.round(::Type{fmpz}, a::fmpq)
  return round(fmpz, a, RoundNearestTiesAway)
end

function Base.round(a::fmpq)
  return round(fmpz, a)
end

################################################################################
#
#  Squarefree part
#
################################################################################

@doc Markdown.doc"""
    squarefree_part(a::fmpz) -> fmpz

Returns the squarefee part $b$ of $a$, which is the smallest (absolute value)
integer $b$ such that $a/b$ is positive and squarefree.
"""
function squarefree_part(a::fmpz)
  f = factor(a)
  s = sign(a)
  for (p, e) in f
    if isodd(e)
      s = s * p
    end
  end
  return s
end

################################################################################
#
#  Factorization of a rational
#
################################################################################

@doc Markdown.doc"""
    factor(a::fmpq, ::FlintIntegerRing) -> Fac{fmpz}

Factor the rational number $a$ into prime numbers.
"""
function factor(a::fmpq, ::FlintIntegerRing)
  fn = factor(numerator(a))
  fd = factor(denominator(a))
  for (p, e) = fd
    fn.fac[p] = -e
  end
  return fn
end

#missing in Nemo...
Hecke.clog(a::Int, b::Int) = clog(fmpz(a), b)

################################################################################
#
#   Support
#
################################################################################

function support(d::fmpz)
  return collect(keys(factor(d).fac))
end

function support(a::fmpq)
  d = denominator(a)
  n = numerator(a)
  res = fmpz[]
  for (p, _) in factor(d)
    push!(res, p)
  end
  for (p, _) in factor(n)
    push!(res, p)
  end
  return res
end


text:
"""
    base_L2!(result, img, depth=length(result))

Compute L2 over vector `img` and *add* it to `result`.

Simple and fast L2 algorithm,
`O(n), n = length(img)`.

Meaningfull values for elements of `img`:
* -1 -- not a part of the image
* 0 -- void
* 1 -- solid, inside original boundaries
* 2 -- solid, outside original boundaries
"""
function base_L2!(result, img, depth=length(result))
    len = 0
    len_in_original = 0

    for (i, x) in enumerate(img)
        if x >= 1
            len += 1
            if x == 1
                len_in_original += 1
            end
        end

        if x < 1 || i == length(img)
            for r in 1:min(depth, len)
                result[r] += min(len - r + 1, len_in_original)
            end
            len = 0
            len_in_original = 0
        end
    end
end

"""
    x_L2!(result, img, depth)

Compute L2 over first dim of `img` and *add* it to `result`.

Simple and fast single-CPU algorithm: `O(N)`,
where `N = length(img)`.
"""
function x_L2!(result, img, depth)
    result .= 0

    indxs = CartesianIndices(size_slice(img, 1))
    for indx in indxs
        data_x = view(img, :, indx)
        result_x = view(result, :, indx)
        base_L2!(result_x, data_x, depth)
    end
    return
end

"""
    x_L2_CUDA!(result, img, depth)

Compute L2 over first dim of `img` and *add* it to `result`
using CUDA.
"""
function x_L2_CUDA!(result, img, depth)
    indxs = CartesianIndices(size_slice(img, 1))
    nindxs = length(indxs)

    cu_index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    cu_stride = blockDim().x * gridDim().x

    for i in cu_index:cu_stride:nindxs
        indx = indxs[i]
        vres = view(result, :, indx)
        vimg = view(img, :, indx)
        base_L2!(vres, vimg, depth)
    end
    return
end

"""
when applied to CuArrays, performs computations on GPU
"""
function x_L2!(result::CuArray, img::CuArray, depth)
    result .= 0

    NTHREADS = 256
    N = reduce(*, size_slice(img, 1); init=1)

    NBLOCKS = div(N - 1, NTHREADS) + 1
    @cuda threads = NTHREADS blocks = NBLOCKS x_L2_CUDA!(result, img, depth)
end

"""
    align!(aimg, img, side, ray; periodic=true)

"""
function align!(aimg, img, side, ray; periodic=true)
# TODO well documented align
n = size(img, side)
    for i in 1:n
        aslice = selectdim(aimg, 1, i)
        slice = selectdim(img, side, i)

        j = round.(Int, (i - 1) .* ray ./ (n - 1))

        if periodic
            circshift!(aslice, slice, .-j)

            aslice2 = selectdim(aimg, 1, i + n)
            k = round.(Int, (n + i - 1) .* ray ./ (n - 1))
            circshift!(aslice2, slice, .-k)
            aslice2 .*= 2
        else
            ix = CartesianIndex(size(slice)) .+ CartesianIndices(slice) .- CartesianIndex(j)
            aslice .= -1
            aslice[ix] .= slice
        end
    end
    aimg
end

"""
    align!(aimg, img, side, ray; periodic=true)

Simple align for 1D
"""
function align!(aimg::AbstractVector, img::AbstractVector, side, ray; periodic=true)
    if periodic
        n = length(img)
        aimg[1:n] .= img
        aimg[n + 1:end] .= 2 .* img
    else
        aimg .= img
    end
    aimg
end

function L2_side!(
    side_result,
    img,
    depth,
    side,
    aligned_result,
    aligned_img;
    periodic::Bool
)
    ray_ixs = CartesianIndices(size_slice(img, side))

    for ray_ix in ray_ixs
        ray_result = view(side_result, :, ray_ix)
        ray_projection = Tuple(ray_ix) .- 1

        align!(aligned_img, img, side, ray_projection; periodic)
        x_L2!(aligned_result, aligned_img, depth)
        sum!(ray_result, aligned_result)
    end
end

function L2_positive_sides!(
    side_results,
    img::AbstractArray{<:Integer,N},
    side_depths,
    side_align_results,
    side_align_imgs;
    periodic
) where N
    for side in 1:N
        result = side_results[side]
        depth = side_depths[side]
        aligned_result = side_align_results[side]
        aligned_img = side_align_imgs[side]
        L2_side!(result, img, depth, side, aligned_result, aligned_img; periodic)
    end
end

function map_ix!(original_ixs, ray_ixs, img)
    N = length(original_ixs)
    for i in 1:N
        original_ixs[i] = []
        ray_ixs[i] = []
    end

    box = size(img)

    for orig_ix in CartesianIndices(img)
        coords = Tuple(orig_ix) .- 1
        side = argmax(coords ./ box)

        n = box[side] - 1
        q = coords[side]

        if q == 0
            rayend = box
        else
            rayend = round.(Int, coords ./ q .* n) .+ 1
        end

        mask = [1:side - 1; side + 1:N]
        rayend = rayend[mask]

        ray_ix = CartesianIndex(q + 1, rayend...)
        push!(original_ixs[side], orig_ix)
        push!(ray_ixs[side], ray_ix)
    end
end

function restore!(result, side_results, original_ixs, ray_ixs)
    for (side_result, orig_ix, ray_ix) in zip(side_results, original_ixs, ray_ixs)
        result[orig_ix] .= side_result[ray_ix]
    end
end

struct Params_L2{Total,Result,AlignImg,N}
    # boundary conditions
    periodic::Bool
    # normalization
    total::Total

    # algorithm-specific

    side_results::Vector{Result}
    side_depths::Tuple{Vararg{Int,N}}
    side_align_imgs::Vector{AlignImg}
    side_align_results::Vector{Result}
    original_ixs::Vector{Vector{CartesianIndex{N}}}
    ray_ixs::Vector{Vector{CartesianIndex{N}}}
end

function Params_L2(img::AbstractArray{<:Integer,N};
            periodic::Bool=true,
            depth::Int=img |> size |> maximum) where N

    side_sizes = [(size(img, d), size_slice(img, d)...) for d in 1:N]
    side_results = map(s -> similar(img, Int64, s), side_sizes)
    side_depths = map(s -> min(s, depth), size(img))

    if periodic
        side_align_img_sizes = map(s -> (2s[1], s[2:end]...), side_sizes)

        side_align_result_sizes = side_sizes

    else
        side_align_img_sizes = map(s -> (s[1], 2 .* s[2:end]...), side_sizes)

        side_align_result_sizes = side_align_img_sizes
    end
    side_align_imgs = map(s -> similar(img, Int8, s), side_align_img_sizes)
    side_align_results = map(s -> similar(img, Int, s...), side_align_result_sizes)

    original_ixs = Vector{Vector{CartesianIndex{N}}}(undef, N)
    ray_ixs = Vector{Vector{CartesianIndex{N}}}(undef, N)
    map_ix!(original_ixs, ray_ixs, img)

    p = Params_L2(
        periodic,
        cnt_total(img, periodic),
        side_results,
        side_depths,
        side_align_imgs,
        side_align_results,
        original_ixs,
        ray_ixs
    )
    cf_type = :central_symmetry
    p, cf_type
end

"""
    correllation_function!(res, img, params::Params_L2)

Compute L₂ correlation function in positive directions
"""
function correllation_function!(res, img, params::Params_L2)
    L2_positive_sides!(
        params.side_results,
        img,
        params.side_depths,
        params.side_align_results,
        params.side_align_imgs;
        params.periodic
    )
    restore!(res, params.side_results, params.original_ixs, params.ray_ixs)
    res ./= params.total
end

"""
    l2(image; periodic = false)

Calculate `L₂` (lineal path) correlation function map
for the 1,2,3-dimensional image and return a `CFMap` object.

The `image` contains the probability of the voxel being in the correct phase.

*For L2, only probabilities 0 and 1 are currently implemented.*

# Examples
```jldoctest
julia> l2([1 0; 0 1]; periodic=true).result
3×2 Matrix{Float32}:
 0.0  0.5
 0.5  0.0
 0.0  0.5
```
"""
function l2(image; periodic=false)
    if !(eltype(image) <: Integer)
        error("For L2, only arrays of integers and Boolean are supported")
    elseif !(eltype(image) == Bool) && any(s -> s ∉ [0, 1], image)
        @warn ("Correct work is possible only for 0 and 1")
    end

    corr_function_map(image, Params_L2; periodic)
end

function l2(image, phase::Int; periodic=false)
    one_phase_img = image .== phase
    corr_function_map(one_phase_img, Params_L2; periodic)
end

ext: jl


text:
"""
    makebuffer(tfm, item)

Create a buffer `buf` that can be used in a call to `apply!(buf, tfm, item)`.
Default to `buffer = apply(tfm, item)`.

You only need to implement this if the default `apply(tfm, item)` isn't
enough. See `apply(tfm::Sequence, item)` for an example of this.
"""
makebuffer(tfm::Transform, items) = apply(tfm, items)

"""
    apply!(buffer::I, tfm, item::I)

Applies `tfm` to `item`, mutating the preallocated `buffer`.

`buffer` can be obtained with `buffer = makebuffer(tfm, item)`

    apply!(buffer, tfm::Transform, item::I; randstate) = apply(tfm, item; randstate)

Default to `apply(tfm, item)` (non-mutating version).
"""
apply!(buf, tfm::Transform, items; randstate = getrandstate(tfm)) = apply(tfm, items, randstate = randstate)
apply!(buf, tfm::Transform, item::Item; randstate = getrandstate(tfm)) = apply(tfm, item, randstate = randstate)
function apply!(bufs::Tuple, tfm::Transform, items::Tuple; randstate = getrandstate(tfm))
    return map((item, buf) -> apply!(buf, tfm, item; randstate = randstate), items, bufs)
end

# ## Buffered transforms

mutable struct Buffered{T<:Transform} <: Transform
    tfm::T
    buffer
    Buffered(tfm::T, buffer = nothing) where T = new{T}(tfm, buffer)
end

getrandstate(buffered::Buffered) = getrandstate(buffered.tfm)

function apply(buffered::Buffered, items::Tuple; randstate = getrandstate(buffered))
    if isnothing(buffered.buffer)
        buffered.buffer = makebuffer(buffered.tfm, items)
    end
    titems = apply!(buffered.buffer, buffered.tfm, items, randstate = randstate)
    return titems
end

apply(buffered::Buffered, item::Item; randstate = getrandstate(buffered)) =
    apply(buffered, (item,); randstate = randstate) |> only

function apply!(buf::Tuple, buffered::Buffered, items::Tuple; randstate = getrandstate(buffered))
    if isnothing(buffered.buffer)
        buffered.buffer = makebuffer(buffered.tfm, items)
    end
    buffered.buffer = apply!(buffered.buffer, buffered.tfm, items; randstate = randstate)
    copyitemdata!(buf, buffered.buffer)
    return buf
end

function apply!(buf::I, buffered::Buffered, item::I; randstate = getrandstate(buffered)) where {I<:Item}
    bufs, items = (buf,), (item,)
    titems = apply!(bufs, buffered, items; randstate = randstate)
    return only(titems)
end

struct BufferedThreadsafe
    buffereds::Vector{Buffered}
    function BufferedThreadsafe(tfm; n = Threads.nthreads())
        @assert n >= 1
        return new([Buffered(tfm) for _ in 1:n])
    end
end

function apply(bufferedt::BufferedThreadsafe, items; kwargs...)
    bufferedthread = bufferedt.buffereds[Threads.threadid()]
    return apply(bufferedthread, items; kwargs...)
end

function apply!(buf, bufferedt::BufferedThreadsafe, items; kwargs...)
    bufferedthread = bufferedt.buffereds[Threads.threadid()]
    return apply!(buf, bufferedthread, items; kwargs...)
end

# Utils

copyitemdata!(buf::I, item::I) where I<:Item = copy!(itemdata(buf), itemdata(item))
copyitemdata!(bufs::T, items::T) where T<:Tuple = (copyitemdata!.(bufs, items); bufs)
copyitemdata!(bufs::T, items::T) where T<:AbstractVector = (copyitemdata!.(bufs, items); bufs)


text:
# Period testing
@test -Dates.Year(1) == Dates.Year(-1)
@test Dates.Year(1) > Dates.Year(0)
@test (Dates.Year(1) < Dates.Year(0)) == false
@test Dates.Year(1) == Dates.Year(1)
@test Dates.Year(1) + Dates.Year(1) == Dates.Year(2)
@test Dates.Year(1) - Dates.Year(1) == Dates.Year(0)
@test Dates.Year(1) * Dates.Year(1) == Dates.Year(1)
@test Dates.Year(10) % Dates.Year(4) == Dates.Year(2)
@test div(Dates.Year(10),Dates.Year(3)) == Dates.Year(3)
@test div(Dates.Year(10),Dates.Year(4)) == Dates.Year(2)
t = Dates.Year(1)
t2 = Dates.Year(2)
@test ([t,t,t,t,t] .+ Dates.Year(1)) == ([t2,t2,t2,t2,t2])
@test (Dates.Year(1) .+ [t,t,t,t,t]) == ([t2,t2,t2,t2,t2])
@test ([t2,t2,t2,t2,t2] .- Dates.Year(1)) == ([t,t,t,t,t])
@test ([t,t,t,t,t] .* Dates.Year(1)) == ([t,t,t,t,t])
@test ([t,t,t,t,t] .% t2) == ([t,t,t,t,t])
@test div([t,t,t,t,t],Dates.Year(1)) == ([t,t,t,t,t])

#Period arithmetic
y = Dates.Year(1)
m = Dates.Month(1)
w = Dates.Week(1)
d = Dates.Day(1)
h = Dates.Hour(1)
mi = Dates.Minute(1)
s = Dates.Second(1)
ms = Dates.Millisecond(1)
@test Dates.Year(y) == y
@test Dates.Month(m) == m
@test Dates.Week(w) == w
@test Dates.Day(d) == d
@test Dates.Hour(h) == h
@test Dates.Minute(mi) == mi
@test Dates.Second(s) == s
@test Dates.Millisecond(ms) == ms
@test typeof(int8(y)) <: Int8
@test typeof(uint8(y)) <: Uint8
@test typeof(int16(y)) <: Int16
@test typeof(uint16(y)) <: Uint16
@test typeof(int32(y)) <: Int32
@test typeof(uint32(y)) <: Uint32
@test typeof(int64(y)) <: Int64
@test typeof(uint64(y)) <: Uint64
@test typeof(int128(y)) <: Int128
@test typeof(uint128(y)) <: Uint128
@test typeof(convert(BigInt,y)) <: BigInt
@test typeof(convert(BigFloat,y)) <: BigFloat
@test typeof(convert(Complex,y)) <: Complex
@test typeof(convert(Rational,y)) <: Rational
@test typeof(float16(y)) <: Float16
@test typeof(float32(y)) <: Float32
@test typeof(float64(y)) <: Float64
@test Dates.Year(convert(Int8,1)) == y
@test Dates.Year(convert(Uint8,1)) == y
@test Dates.Year(convert(Int16,1)) == y
@test Dates.Year(convert(Uint16,1)) == y
@test Dates.Year(convert(Int32,1)) == y
@test Dates.Year(convert(Uint32,1)) == y
@test Dates.Year(convert(Int64,1)) == y
@test Dates.Year(convert(Uint64,1)) == y
@test Dates.Year(convert(Int128,1)) == y
@test Dates.Year(convert(Uint128,1)) == y
@test Dates.Year(convert(BigInt,1)) == y
@test Dates.Year(convert(BigFloat,1)) == y
@test Dates.Year(convert(Complex,1)) == y
@test Dates.Year(convert(Rational,1)) == y
@test Dates.Year(convert(Float16,1)) == y
@test Dates.Year(convert(Float32,1)) == y
@test Dates.Year(convert(Float64,1)) == y
@test y == y
@test m == m
@test w == w
@test d == d
@test h == h
@test mi == mi
@test s == s
@test ms == ms
@test_throws ArgumentError y != m
@test_throws ArgumentError m != w
@test_throws ArgumentError w != d
@test_throws ArgumentError d != h
@test_throws ArgumentError h != mi
@test_throws ArgumentError mi != s
@test_throws ArgumentError s != ms
@test_throws ArgumentError ms != y
y2 = Dates.Year(2)
@test y < y2
@test y2 > y
@test y != y2

@test Dates.Year(int8(1)) == y
@test Dates.Year(uint8(1)) == y
@test Dates.Year(int16(1)) == y
@test Dates.Year(uint16(1)) == y
@test Dates.Year(int(1)) == y
@test Dates.Year(uint(1)) == y
@test Dates.Year(int64(1)) == y
@test Dates.Year(uint64(1)) == y
@test Dates.Year(int128(1)) == y
@test Dates.Year(uint128(1)) == y
@test Dates.Year(big(1)) == y
@test Dates.Year(BigFloat(1)) == y
@test Dates.Year(float(1)) == y
@test Dates.Year(float32(1)) == y
@test Dates.Year(Rational(1)) == y
@test Dates.Year(complex(1)) == y
@test_throws InexactError Dates.Year(BigFloat(1.2)) == y
@test_throws InexactError Dates.Year(1.2) == y
@test_throws InexactError Dates.Year(float32(1.2)) == y
@test_throws InexactError Dates.Year(3//4) == y
@test_throws InexactError Dates.Year(complex(1.2)) == y
@test_throws InexactError Dates.Year(float16(1.2)) == y
@test Dates.Year(true) == y
@test Dates.Year(false) != y
@test Dates.Year('\x01') == y
@test_throws MethodError Dates.Year(:hey) == y
@test Dates.Year(real(1)) == y
@test_throws ArgumentError Dates.Year(m) == y
@test_throws ArgumentError Dates.Year(w) == y
@test_throws ArgumentError Dates.Year(d) == y
@test_throws ArgumentError Dates.Year(h) == y
@test_throws ArgumentError Dates.Year(mi) == y
@test_throws ArgumentError Dates.Year(s) == y
@test_throws ArgumentError Dates.Year(ms) == y
@test Dates.Year(Dates.Date(2013,1,1)) == Dates.Year(2013)
@test Dates.Year(Dates.DateTime(2013,1,1)) == Dates.Year(2013)
@test typeof(y+m) <: Dates.CompoundPeriod
@test typeof(m+y) <: Dates.CompoundPeriod
@test typeof(y+w) <: Dates.CompoundPeriod
@test typeof(y+d) <: Dates.CompoundPeriod
@test typeof(y+h) <: Dates.CompoundPeriod
@test typeof(y+mi) <: Dates.CompoundPeriod
@test typeof(y+s) <: Dates.CompoundPeriod
@test typeof(y+ms) <: Dates.CompoundPeriod
@test_throws ArgumentError y > m
@test_throws ArgumentError d < w
@test typemax(Dates.Year) == Dates.Year(typemax(Int64))
@test typemax(Dates.Year) + y == Dates.Year(-9223372036854775808)
@test typemin(Dates.Year) == Dates.Year(-9223372036854775808)
#Period-Real arithmetic
@test y + 1 == Dates.Year(2)
@test 1 + y == Dates.Year(2)
@test y + true == Dates.Year(2)
@test true + y == Dates.Year(2)
@test y + '\x01' == Dates.Year(2)
@test '\x01' + y == Dates.Year(2)
@test y + 1.0 == Dates.Year(2)
@test_throws InexactError y + 1.2
@test y + 1f0 == Dates.Year(2)
@test_throws InexactError y + 1.2f0
@test y + BigFloat(1) == Dates.Year(2)
@test_throws InexactError y + BigFloat(1.2)
@test y + 1.0 == Dates.Year(2)
@test_throws InexactError y + 1.2
@test y * 4 == Dates.Year(4)
@test y * 4f0 == Dates.Year(4)
@test_throws InexactError y * 3//4 == Dates.Year(1)
@test div(y,2) == Dates.Year(0)
@test div(2,y) == Dates.Year(2)
@test div(y,y) == Dates.Year(1)
@test y*10 % 5 == Dates.Year(0)
@test 5 % y*10 == Dates.Year(0)
@test (y > 3) == false
@test (4 < y) == false
@test 1 == y
t = [y,y,y,y,y]
@test t .+ Dates.Year(2) == [Dates.Year(3),Dates.Year(3),Dates.Year(3),Dates.Year(3),Dates.Year(3)]
dt = Dates.DateTime(2012,12,21)
# Associativity
test = ((((((((dt + y) - m) + w) - d) + h) - mi) + s) - ms)
@test test == dt + y - m + w - d + h - mi + s - ms
@test test == y - m + w - d + dt + h - mi + s - ms
@test test == dt - m + y - d + w - mi + h - ms + s
@test test == dt + (y - m + w - d + h - mi + s - ms)
@test test == dt + y - m + w - d + (h - mi + s - ms)
@test (dt + Dates.Year(4)) + Dates.Day(1) == dt + (Dates.Year(4) + Dates.Day(1))
@test Dates.Date(2014,1,29) + Dates.Month(1) + Dates.Day(1) + Dates.Month(1) + Dates.Day(1) ==
    Dates.Date(2014,1,29) + Dates.Day(1) + Dates.Month(1) + Dates.Month(1) + Dates.Day(1)
@test Dates.Date(2014,1,29) + Dates.Month(1) + Dates.Day(1) == Dates.Date(2014,1,29) + Dates.Day(1) + Dates.Month(1)
# traits
@test Dates._units(Dates.Year(0)) == " years"
@test Dates._units(Dates.Year(1)) == " year"
@test Dates._units(Dates.Year(-1)) == " year"
@test Dates._units(Dates.Year(2)) == " years"
@test Dates.string(Dates.Year(0)) == "0 years"
@test Dates.string(Dates.Year(1)) == "1 year"
@test Dates.string(Dates.Year(-1)) == "-1 year"
@test Dates.string(Dates.Year(2)) == "2 years"
@test zero(Dates.Year) == Dates.Year(0)
@test zero(Dates.Year(10)) == Dates.Year(0)
@test zero(Dates.Month) == Dates.Month(0)
@test zero(Dates.Month(10)) == Dates.Month(0)
@test zero(Dates.Day) == Dates.Day(0)
@test zero(Dates.Day(10)) == Dates.Day(0)
@test zero(Dates.Hour) == Dates.Hour(0)
@test zero(Dates.Hour(10)) == Dates.Hour(0)
@test zero(Dates.Minute) == Dates.Minute(0)
@test zero(Dates.Minute(10)) == Dates.Minute(0)
@test zero(Dates.Second) == Dates.Second(0)
@test zero(Dates.Second(10)) == Dates.Second(0)
@test zero(Dates.Millisecond) == Dates.Millisecond(0)
@test zero(Dates.Millisecond(10)) == Dates.Millisecond(0)
@test one(Dates.Year) == Dates.Year(1)
@test one(Dates.Year(10)) == Dates.Year(1)
@test one(Dates.Month) == Dates.Month(1)
@test one(Dates.Month(10)) == Dates.Month(1)
@test one(Dates.Day) == Dates.Day(1)
@test one(Dates.Day(10)) == Dates.Day(1)
@test one(Dates.Hour) == Dates.Hour(1)
@test one(Dates.Hour(10)) == Dates.Hour(1)
@test one(Dates.Minute) == Dates.Minute(1)
@test one(Dates.Minute(10)) == Dates.Minute(1)
@test one(Dates.Second) == Dates.Second(1)
@test one(Dates.Second(10)) == Dates.Second(1)
@test one(Dates.Millisecond) == Dates.Millisecond(1)
@test one(Dates.Millisecond(10)) == Dates.Millisecond(1)
@test Dates.Year(-1) < Dates.Year(1)
@test !(Dates.Year(-1) > Dates.Year(1))
@test Dates.Year(1) == Dates.Year(1)
@test Dates.Year(1) == 1
@test 1 == Dates.Year(1)
@test (Dates.Year(1) < 1) == false
@test (1 < Dates.Year(1)) == false
@test Dates.Month(-1) < Dates.Month(1)
@test !(Dates.Month(-1) > Dates.Month(1))
@test Dates.Month(1) == Dates.Month(1)
@test Dates.Month(1) == 1
@test 1 == Dates.Month(1)
@test (Dates.Month(1) < 1) == false
@test (1 < Dates.Month(1)) == false
@test Dates.Day(-1) < Dates.Day(1)
@test !(Dates.Day(-1) > Dates.Day(1))
@test Dates.Day(1) == Dates.Day(1)
@test Dates.Day(1) == 1
@test 1 == Dates.Day(1)
@test (Dates.Day(1) < 1) == false
@test (1 < Dates.Day(1)) == false
@test Dates.Hour(-1) < Dates.Hour(1)
@test !(Dates.Hour(-1) > Dates.Hour(1))
@test Dates.Hour(1) == Dates.Hour(1)
@test Dates.Hour(1) == 1
@test 1 == Dates.Hour(1)
@test (Dates.Hour(1) < 1) == false
@test (1 < Dates.Hour(1)) == false
@test Dates.Minute(-1) < Dates.Minute(1)
@test !(Dates.Minute(-1) > Dates.Minute(1))
@test Dates.Minute(1) == Dates.Minute(1)
@test Dates.Minute(1) == 1
@test 1 == Dates.Minute(1)
@test (Dates.Minute(1) < 1) == false
@test (1 < Dates.Minute(1)) == false
@test Dates.Second(-1) < Dates.Second(1)
@test !(Dates.Second(-1) > Dates.Second(1))
@test Dates.Second(1) == Dates.Second(1)
@test Dates.Second(1) == 1
@test 1 == Dates.Second(1)
@test (Dates.Second(1) < 1) == false
@test (1 < Dates.Second(1)) == false
@test Dates.Millisecond(-1) < Dates.Millisecond(1)
@test !(Dates.Millisecond(-1) > Dates.Millisecond(1))
@test Dates.Millisecond(1) == Dates.Millisecond(1)
@test Dates.Millisecond(1) == 1
@test 1 == Dates.Millisecond(1)
@test (Dates.Millisecond(1) < 1) == false
@test (1 < Dates.Millisecond(1)) == false
@test_throws ArgumentError Dates.Year(1) < Dates.Millisecond(1)
@test_throws ArgumentError Dates.Year(1) == Dates.Millisecond(1)

@test Dates.Year("1") == y
@test Dates.Month("1") == m
@test Dates.Week("1") == w
@test Dates.Day("1") == d
@test Dates.Hour("1") == h
@test Dates.Minute("1") == mi
@test Dates.Second("1") == s
@test Dates.Millisecond("1") == ms
@test_throws ErrorException Dates.Year("1.0")
@test Dates.Year(parsefloat("1.0")) == y

dt = Dates.DateTime(2014)
@test typeof(Dates.Year(dt)) <: Dates.Year
@test typeof(Dates.Month(dt)) <: Dates.Month
@test typeof(Dates.Week(dt)) <: Dates.Week
@test typeof(Dates.Day(dt)) <: Dates.Day
@test typeof(Dates.Hour(dt)) <: Dates.Hour
@test typeof(Dates.Minute(dt)) <: Dates.Minute
@test typeof(Dates.Second(dt)) <: Dates.Second
@test typeof(Dates.Millisecond(dt)) <: Dates.Millisecond

# Default values
@test Dates.default(Dates.Year) == y
@test Dates.default(Dates.Month) == m
@test Dates.default(Dates.Week) == w
@test Dates.default(Dates.Day) == d
@test Dates.default(Dates.Hour) == zero(Dates.Hour)
@test Dates.default(Dates.Minute) == zero(Dates.Minute)
@test Dates.default(Dates.Second) == zero(Dates.Second)
@test Dates.default(Dates.Millisecond) == zero(Dates.Millisecond)

# Conversions
@test Dates.toms(ms) == Dates.value(Dates.Millisecond(ms)) == 1
@test Dates.toms(s)  == Dates.value(Dates.Millisecond(s)) == 1000
@test Dates.toms(mi) == Dates.value(Dates.Millisecond(mi)) == 60000
@test Dates.toms(h)  == Dates.value(Dates.Millisecond(h)) == 3600000
@test Dates.toms(d)  == Dates.value(Dates.Millisecond(d)) == 86400000
@test Dates.toms(w)  == Dates.value(Dates.Millisecond(w)) == 604800000

@test Dates.days(ms) == Dates.days(s) == Dates.days(mi) == Dates.days(h) == 0
@test Dates.days(Dates.Millisecond(86400000)) == 1
@test Dates.days(Dates.Second(86400)) == 1
@test Dates.days(Dates.Minute(1440)) == 1
@test Dates.days(Dates.Hour(24)) == 1
@test Dates.days(d) == 1
@test Dates.days(w) == 7


text:
include("test_stencil.jl")
include("test_core.jl")


text:
struct NewCodeInfo
    src::CodeInfo
    code::Vector{Any}
    nvariables::Int
    codelocs::Vector{Int32}
    newslots::Dict{Int,Symbol}
    slotnames::Vector{Symbol}
    changemap::Vector{Int}
    slotmap::Vector{Int}

    function NewCodeInfo(ci::CodeInfo, nargs::Int)
        code = []
        codelocs = Int32[]
        newslots = Dict{Int,Symbol}()
        slotnames = copy(ci.slotnames)
        changemap = fill(0, length(ci.code))
        slotmap = fill(0, length(ci.slotnames))
        new(ci, code, nargs + 1, codelocs, newslots, slotnames, changemap, slotmap)
    end
end

source_slot(ci::NewCodeInfo, i::Int) = Core.SlotNumber(i + ci.slotmap[i])

function slot(ci::NewCodeInfo, name::Symbol)
    return Core.SlotNumber(findfirst(isequal(name), ci.slotnames))
end

function unpack_closure!(ci::NewCodeInfo, closure::Int)
    spec = Core.SlotNumber(closure)
    codeloc = ci.src.codelocs[1]
    # unpack closure
    # %1 = get variables
    push!(ci.code, Expr(:call, GlobalRef(Base, :getfield), spec, QuoteNode(:variables)))
    push!(ci.codelocs, codeloc)
    ci.changemap[1] += 1

    # %2 = get parent
    push!(
        ci.code,
        Expr(
            :(=),
            source_slot(ci, 2),
            Expr(:call, GlobalRef(Base, :getfield), spec, QuoteNode(:parent)),
        ),
    )
    push!(ci.codelocs, codeloc)
    # unpack variables
    for i in 2:ci.nvariables
        push!(
            ci.code,
            Expr(
                :(=),
                source_slot(ci, i + 1),
                Expr(:call, GlobalRef(Base, :getindex), NewSSAValue(1), i - 1),
            ),
        )
        push!(ci.codelocs, codeloc)
    end
    ci.changemap[1] += ci.nvariables
    return ci
end

function insert_slot!(ci::NewCodeInfo, v::Int, slot::Symbol)
    ci.newslots[v] = slot
    insert!(ci.slotnames, v, slot)
    prev = length(filter(x -> x < v, keys(ci.newslots)))
    for k in v-prev:length(ci.slotmap)
        ci.slotmap[k] += 1
    end
    return ci
end

function push_stmt!(ci::NewCodeInfo, stmt, codeloc::Int32 = Int32(1))
    push!(ci.code, stmt)
    push!(ci.codelocs, codeloc)
    return ci
end

function insert_stmt!(ci::NewCodeInfo, v::Int, stmt)
    push_stmt!(ci, stmt, ci.src.codelocs[v])
    ci.changemap[v] += 1
    return NewSSAValue(length(ci.code))
end

function update_slots(e, slotmap)
    if e isa Core.SlotNumber
        return Core.SlotNumber(e.id + slotmap[e.id])
    elseif e isa Expr
        return Expr(e.head, map(x -> update_slots(x, slotmap), e.args)...)
    elseif e isa Core.NewvarNode
        return Core.NewvarNode(Core.SlotNumber(e.slot.id + slotmap[e.slot.id]))
    else
        return e
    end
end

function finish(ci::NewCodeInfo)
    Core.Compiler.renumber_ir_elements!(ci.code, ci.changemap)
    replace_new_ssavalue!(ci.code)
    new_ci = copy(ci.src)
    new_ci.code = ci.code
    new_ci.codelocs = ci.codelocs
    new_ci.slotnames = ci.slotnames
    new_ci.slotflags = [0x00 for _ in new_ci.slotnames]
    new_ci.inferred = false
    new_ci.inlineable = true
    new_ci.ssavaluetypes = length(ci.code)
    return new_ci
end

function is_quantum_statement(@nospecialize(e))
    e isa Function && parentmodule(e) == Semantic && return true
    e isa Expr || return false
    e.head === :quantum && return true
    if e.head === :call
        f = e.args[1]
        f isa Function && parentmodule(f) === Semantic && return true
        f isa GlobalRef && f.mod === Semantic && return true
    elseif e.head === :invoke
        f = e.args[2]
        f isa Function && parentmodule(f) === Semantic && return true
        f isa GlobalRef && f.mod === Semantic && return true
    elseif e.head === :(=)
        return is_quantum_statement(e.args[2])
    end
    return false
end

function quantum_stmt_type(e::Expr)
    if e.head === :call
        e.args[1] isa Function && return nameof(e.args[1])
        return e.args[1].name
    elseif e.head === :invoke
        e.args[2] isa Function && return nameof(e.args[2])
        return e.args[2].name
    elseif e.head === :quantum
        return e.args[1]
    else
        return quantum_stmt_type(e.args[2])
    end
end

function _replace_new_ssavalue(e)
    if e isa NewSSAValue
        return SSAValue(e.id)
    elseif e isa Expr
        return Expr(e.head, map(_replace_new_ssavalue, e.args)...)
    elseif e isa Core.GotoIfNot
        cond = e.cond
        if cond isa NewSSAValue
            cond = SSAValue(cond.id)
        end
        return Core.GotoIfNot(cond, e.dest)
    elseif e isa Core.ReturnNode && isdefined(e, :val) && isa(e.val, NewSSAValue)
        return Core.ReturnNode(SSAValue(e.val.id))
    else
        return e
    end
end

function replace_new_ssavalue!(code::Vector)
    for idx in 1:length(code)
        code[idx] = _replace_new_ssavalue(code[idx])
    end
    return code
end

function obtain_codeinfo(::Type{RoutineSpec{P,Sigs}}) where {P,Sigs}
    nargs = length(Sigs.parameters)
    tt = Tuple{P,Sigs.parameters...}
    ms = methods(routine_stub, tt)
    @assert length(ms) == 1
    method = first(ms)
    method_args = Tuple{RoutineStub,tt.parameters...}
    mi = Core.Compiler.specialize_method(method, method_args, Core.svec())
    ci = Core.Compiler.retrieve_code_info(mi)

    name = routine_name(P)
    linetable = Any[]
    for lineinfo in ci.linetable
        push!(
            linetable,
            Core.LineInfoNode(
                lineinfo.module,
                name,
                lineinfo.file,
                lineinfo.line,
                lineinfo.inlined_at,
            ),
        )
    end
    ci.linetable = linetable
    ci.edges = Core.MethodInstance[mi]
    return ci, nargs
end

function create_codeinfo(::typeof(Semantic.main), S::Type{<:RoutineSpec})
    ci, nargs = obtain_codeinfo(S)
    new = NewCodeInfo(ci, nargs)
    insert_slot!(new, 2, :spec)
    unpack_closure!(new, 2)

    for (v, stmt) in enumerate(ci.code)
        push_stmt!(new, update_slots(stmt, new.slotmap), ci.codelocs[v])
    end
    return finish(new)
end

function _extract_measure(e)
    if e.head === :(=)
        return e.args[1], e.args[2]
    else
        return nothing, e
    end
end

function create_codeinfo(::typeof(Semantic.gate), S::Type{<:RoutineSpec})
    ci, nargs = obtain_codeinfo(S)
    new = NewCodeInfo(ci, nargs)
    insert_slot!(new, 2, :spec)
    insert_slot!(new, 3, :locations)
    unpack_closure!(new, 2)
    locations = slot(new, :locations)

    for (v, stmt) in enumerate(new.src.code)
        codeloc = new.src.codelocs[v]
        stmt = update_slots(stmt, new.slotmap)
        e = nothing
        if is_quantum_statement(stmt)
            type = quantum_stmt_type(stmt)
            if type === :gate
                local_location = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[3]),
                )
                e = Expr(:call, GlobalRef(Semantic, :gate), stmt.args[2], local_location)
            elseif type === :ctrl
                local_location = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[3]),
                )
                local_ctrl = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[4]),
                )
                e = Expr(:call, GlobalRef(Semantic, :ctrl), stmt.args[2], local_location, local_ctrl)
            elseif type === :measure
                cvar, measure = _extract_measure(stmt)
                # no location specified
                if length(measure.args) == 2
                    measure_locs = locations
                else
                    measure_locs = insert_stmt!(
                        new,
                        v,
                        Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[3]),
                    )
                end

                # TODO: handle measure operator
                measure_ex = Expr(:call, GlobalRef(Semantic, :measure), measure_locs)
                if isnothing(cvar)
                    e = measure_ex
                else
                    e = Expr(:(=), cvar, measure_ex)
                end
            elseif type === :barrier
                local_location = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[2]),
                )
                e = Expr(:call, GlobalRef(Semantic, :barrier), local_location)
            end
        end

        if isnothing(e)
            push_stmt!(new, stmt, codeloc)
        else
            push_stmt!(new, e, codeloc)
        end
    end
    return finish(new)
end

function create_codeinfo(::typeof(Semantic.ctrl), S::Type{<:RoutineSpec})
    ci, nargs = obtain_codeinfo(S)
    new = NewCodeInfo(ci, nargs)
    insert_slot!(new, 2, :spec)
    insert_slot!(new, 3, :locations)
    insert_slot!(new, 4, :ctrl)
    unpack_closure!(new, 2)
    locations = slot(new, :locations)
    ctrl = slot(new, :ctrl)

    for (v, stmt) in enumerate(new.src.code)
        codeloc = new.src.codelocs[v]
        stmt = update_slots(stmt, new.slotmap)
        e = nothing
        if is_quantum_statement(stmt)
            type = quantum_stmt_type(stmt)
            if type === :gate
                local_location = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[3]),
                )
                e = Expr(:call, GlobalRef(Semantic, :ctrl), stmt.args[2], local_location, ctrl)
            elseif type === :ctrl
                local_location = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[3]),
                )
                local_ctrl = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[4]),
                )
                real_ctrl = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(YaoCompiler, :merge_locations), local_ctrl, ctrl),
                )
                e = Expr(:call, GlobalRef(Semantic, :ctrl), stmt.args[2], local_location, real_ctrl)
            elseif type === :measure
                error("cannot use measure under a quantum control context")
            elseif type === :barrier
                local_location = insert_stmt!(
                    new,
                    v,
                    Expr(:call, GlobalRef(Base, :getindex), locations, stmt.args[2]),
                )
                e = Expr(:call, GlobalRef(Semantic, :barrier), local_location)
            end
        end

        if isnothing(e)
            push_stmt!(new, stmt, codeloc)
        else
            push_stmt!(new, e, codeloc)
        end
    end
    return finish(new)
end

# NOTE: these two functions are mainly compile time stubs
# so we can attach this piece of CodeInfo to certain method
@generated function Semantic.main(spec::RoutineSpec)
    return create_codeinfo(Semantic.main, spec)
end

@generated function Semantic.gate(spec::RoutineSpec, ::Locations)
    return create_codeinfo(Semantic.gate, spec)
end

@generated function Semantic.ctrl(spec::RoutineSpec, ::Locations, ::CtrlLocations)
    return create_codeinfo(Semantic.ctrl, spec)
end

function _prepare_frame(f, spec, args...)
    method = methods(f, Tuple{spec,args...}) |> first
    atypes = Tuple{typeof(f),spec,args...}
    mi = Core.Compiler.specialize_method(method, atypes, Core.svec())
    result = Core.Compiler.InferenceResult(mi, Any[Core.Const(f), spec, args...])
    world = Core.Compiler.get_world_counter()
    interp = YaoCompiler.YaoInterpreter(;)
    frame = Core.Compiler.InferenceState(result, true, interp) #=cached=#
    return interp, frame
end


text:
using Test
using BigArrays
using BigArrays.S3Dicts
using OffsetArrays

d = S3Dict( "s3://seunglab/jpwu/test/image/" )
ba = BigArray(d)

@show ndims(ba)
@show size(ba)
@show eltype(ba)

a = rand(UInt8, 256,256,16)

@testset "test 3D image reading and saving" begin
    # ba = BigArray(d, UInt8, (128,128,8))
    @time ba[257:512, 257:512, 17:32] = a
    # BigArrays.mysetindex!(ba, a, (201:400, 201:400, 161:116))
    @time b = ba[257:512, 257:512, 17:32] |> parent
    @show a[:][end-10:end]
    @show b[:][end-10:end]
    @test all(a.==b)
end

@testset "test 3D UInt32 segmentation reading and saving" begin
    d = S3Dict( "s3://seunglab/jpwu/test/segmentation/" )
    # ba = BigArray( d, configDict )
    ba = BigArray(d)

    a = rand(UInt32, 256,256,16)
    @time ba[257:512, 257:512, 17:32] = a
    @time b = ba[257:512, 257:512, 17:32] |> parent
    @test all(a.==b)
end

@testset "test UInt64 segmenation with uint64" begin
    d = S3Dict( "s3://seunglab/jpwu/test/segmentation-uint64/" )
    ba = BigArray(d)

    a = rand(UInt64, 256,256,16)
    @time ba[257:512, 257:512, 17:32] = a
    @time b = ba[257:512, 257:512, 17:32] |> parent
    @test all(a.==b)
end

@testset "test affinity map" begin
    d = S3Dict( "s3://seunglab/jpwu/test/affinitymap/" )
    ba = BigArray(d)

    a = rand(Float32, 256,256,16,3)
    @time ba[257:512, 257:512, 17:32, 1:3] = a
    @time b = ba[257:512, 257:512, 17:32, 1:3] |> parent

    @show size(a)
    @show size(b)
    @test all(a.==b)
end

@testset "test semantic map" begin
    d = S3Dict( "s3://seunglab/jpwu/test/semanticmap/" )
    a = rand(Float32, 256,256,16,4)
    ba = BigArray(d)

    @time ba[257:512, 257:512, 17:32, 1:4] = a
    @time b = ba[257:512, 257:512, 17:32, 1:4] |> parent

    @show size(a)
    @show size(b)
    @test all(a.==b)
end


text:
# =========================================================================== #
# Compliant julia 1.x

# Using the following packages
using JuMP, GLPK
using LinearAlgebra

include("loadSPP.jl")
include("setSPP.jl")
include("getfname.jl")
#include("DM1_1.jl")
#include("DM1_2.jl")
#include("DM2_1.jl")
include("DM2_2.jl")
#include("brouillon.jl")

# =========================================================================== #

# Setting the data
#fname = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data/didactic.dat"  # path for a standard config on macOS
#cost, matrix, n, m = loadSPP(fname)

# Proceeding to the optimization
#solverSelected = GLPK.Optimizer
#ip, ip_x = setSPP(solverSelected, cost, matrix)
#println("Solving..."); optimize!(ip)

# Displaying the results
#println("z  = ", objective_value(ip))
#sprint("x  = "); println(value.(ip_x))

# =========================================================================== #

# Collecting the names of instances to solve
target = "B:/Cours/Nantes/Metaheuristique/DM1_Metaheuristique/solveSPP-master/Data"            # path for a standard config on macOS
fnames = getfname(target)

#===================================================#
for f in fnames
    println("=================")
    println(f)
    cost, matrix, n, m = loadSPP(f)

    #for i in 1:5
    #    @time (SOL,z) =Glouton(cost, matrix, n, m)
    #    @time (SOL,z) =GRASP(cost, matrix, n, m, 0.9)
    #    @time (SOL,z) =Glouton(cost, matrix, n, m)
        @time proba = ReactiveGRASP(matrix,cost, n, m, 15, 5, 1)
    #end
end
cd("../")

# ============================================================================ #
#n = 9
#m = 7
#a  = zeros(Float64, n)
#cost = [10, 5, 8, 6, 9, 13, 11, 4, 6]
#matrix = [1 1 1 0 1 0 1 1 0; 0 1 1 0 0 0 0 1 0;0 1 0 0 1 1 0 1 1;0 0 0 1 0 0 0 0 0;1 0 1 0 1 1 0 0 1;0 1 1 0 0 0 1 0 1; 1 0 0 1 1 0 0 1 1]
#a = Glouton(cost, matrix, n, m)
#println("Solution: ",a)

max_stars_repo_path: Archives/before the merge/main.jl

max_issues_repo_path: Archives/before the merge/main.jl

max_forks_repo_path: Archives/before the merge/main.jl


text:
get_connection_type(s) = getmetadata(unwrap(s), VariableConnectType, Equality)

function with_connector_type(expr)
    @assert expr isa Expr && (expr.head == :function || (expr.head == :(=) &&
                                       expr.args[1] isa Expr &&
                                       expr.args[1].head == :call))

    sig = expr.args[1]
    body = expr.args[2]

    fname = sig.args[1]
    args = sig.args[2:end]

    quote
        function $fname($(args...))
            function f()
                $body
            end
            res = f()
            $isdefined(res, :connector_type) && $getfield(res, :connector_type) === nothing ? $Setfield.@set!(res.connector_type = $connector_type(res)) : res
        end
    end
end

macro connector(expr)
    esc(with_connector_type(expr))
end

abstract type AbstractConnectorType end
struct StreamConnector <: AbstractConnectorType end
struct RegularConnector <: AbstractConnectorType end

function connector_type(sys::AbstractSystem)
    sts = get_states(sys)
    n_stream = 0
    n_flow = 0
    n_regular = 0 # state that is not input, output, stream, or flow.
    for s in sts
        vtype = get_connection_type(s)
        if vtype === Stream
            isarray(s) && error("Array stream variables are not supported. Got $s.")
            n_stream += 1
        elseif vtype === Flow
            n_flow += 1
        elseif !(isinput(s) || isoutput(s))
            n_regular += 1
        end
    end
    (n_stream > 0 && n_flow > 1) && error("There are multiple flow variables in $(nameof(sys))!")
    if n_flow != n_regular
        @warn "$(nameof(sys)) contains $n_flow variables, yet $n_regular regular " *
        "(non-flow, non-stream, non-input, non-output) variables." *
        "This could lead to imbalanced model that are difficult to debug." *
        "Consider marking some of the regular variables as input/output variables."
    end
    n_stream > 0 ? StreamConnector() : RegularConnector()
end

get_systems(c::Connection) = c.systems

instream(a) = term(instream, unwrap(a), type=symtype(a))
SymbolicUtils.promote_symtype(::typeof(instream), _) = Real

isconnector(s::AbstractSystem) = has_connector_type(s) && get_connector_type(s) !== nothing

function flowvar(sys::AbstractSystem)
    sts = get_states(sys)
    for s in sts
        vtype = get_connection_type(s)
        vtype === Flow && return s
    end
    error("There in no flow variable in $(nameof(sys))")
end

collect_instream!(set, eq::Equation) = collect_instream!(set, eq.lhs) | collect_instream!(set, eq.rhs)

function collect_instream!(set, expr, occurs=false)
    istree(expr) || return occurs
    op = operation(expr)
    op === instream && (push!(set, expr); occurs = true)
    for a in SymbolicUtils.unsorted_arguments(expr)
        occurs |= collect_instream!(set, a, occurs)
    end
    return occurs
end

#positivemax(m, ::Any; tol=nothing)= max(m, something(tol, 1e-8))
#_positivemax(m, tol) = ifelse((-tol <= m) & (m <= tol), ((3 * tol - m) * (tol + m)^3)/(16 * tol^3) + tol, max(m, tol))
function _positivemax(m, si)
    T = typeof(m)
    relativeTolerance = 1e-4
    nominal = one(T)
    eps = relativeTolerance * nominal
    alpha = if si > eps
        one(T)
    else
        if si > 0
            (si/eps)^2*(3-2* si/eps)
        else
            zero(T)
        end
    end
    alpha * max(m, 0) + (1-alpha)*eps
end
@register _positivemax(m, tol)
positivemax(m, ::Any; tol=nothing) = _positivemax(m, tol)
mydiv(num, den) = if den == 0
    error()
else
    num / den
end
@register mydiv(n, d)

function generate_isouter(sys::AbstractSystem)
    outer_connectors = Symbol[]
    for s in get_systems(sys)
        n = nameof(s)
        isconnector(s) && push!(outer_connectors, n)
    end
    let outer_connectors=outer_connectors
        function isouter(sys)::Bool
            s = string(nameof(sys))
            isconnector(sys) || error("$s is not a connector!")
            idx = findfirst(isequal('₊'), s)
            parent_name = Symbol(idx === nothing ? s : s[1:prevind(s, idx)])
            parent_name in outer_connectors
        end
    end
end

struct LazyNamespace
    namespace::Union{Nothing,Symbol}
    sys
end

Base.copy(l::LazyNamespace) = renamespace(l.namespace, l.sys)
Base.nameof(l::LazyNamespace) = renamespace(l.namespace, nameof(l.sys))

struct ConnectionElement
    sys::LazyNamespace
    v
    isouter::Bool
end
Base.hash(l::ConnectionElement, salt::UInt) = hash(nameof(l.sys)) ⊻ hash(l.v) ⊻ hash(l.isouter) ⊻ salt
Base.isequal(l1::ConnectionElement, l2::ConnectionElement) = l1 == l2
Base.:(==)(l1::ConnectionElement, l2::ConnectionElement) = nameof(l1.sys) == nameof(l2.sys) && isequal(l1.v, l2.v) && l1.isouter == l2.isouter
namespaced_var(l::ConnectionElement) = states(l, l.v)
states(l::ConnectionElement, v) = states(copy(l.sys), v)

struct ConnectionSet
    set::Vector{ConnectionElement} # namespace.sys, var, isouter
end

function Base.show(io::IO, c::ConnectionSet)
    print(io, "<")
    for i in 1:length(c.set)-1
        @unpack sys, v, isouter = c.set[i]
        print(io, nameof(sys), ".", v, "::", isouter ? "outer" : "inner", ", ")
    end
    @unpack sys, v, isouter = last(c.set)
    print(io, nameof(sys), ".", v, "::", isouter ? "outer" : "inner", ">")
end

@noinline connection_error(ss) = error("Different types of connectors are in one conenction statement: <$(map(nameof, ss))>")

function connection2set!(connectionsets, namespace, ss, isouter)
    nn = map(nameof, ss)
    sts1 = Set(states(first(ss)))
    T = ConnectionElement
    csets = [T[] for _ in 1:length(sts1)]
    for (i, s) in enumerate(ss)
        sts = states(s)
        i != 1 && ((length(sts1) == length(sts) && all(Base.Fix2(in, sts1), sts)) || connection_error(ss))
        io = isouter(s)
        for (j, v) in enumerate(sts)
            push!(csets[j], T(LazyNamespace(namespace, s), v, io))
        end
    end
    for cset in csets
        vtype = get_connection_type(first(cset).v)
        for k in 2:length(cset)
            vtype === get_connection_type(cset[k].v) || connection_error(ss)
        end
        push!(connectionsets, ConnectionSet(cset))
    end
end

function generate_connection_set(sys::AbstractSystem)
    connectionsets = ConnectionSet[]
    sys = generate_connection_set!(connectionsets, sys)
    sys, merge(connectionsets)
end

function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace=nothing)
    subsys = get_systems(sys)

    isouter = generate_isouter(sys)
    eqs′ = get_eqs(sys)
    eqs = Equation[]

    cts = [] # connections
    for eq in eqs′
        if eq.lhs isa Connection
            push!(cts, get_systems(eq.rhs))
        else
            push!(eqs, eq) # split connections and equations
        end
    end

    # all connectors are eventually inside connectors.
    T = ConnectionElement
    for s in subsys
        isconnector(s) || continue
        for v in states(s)
            Flow === get_connection_type(v) || continue
            push!(connectionsets, ConnectionSet([T(LazyNamespace(namespace, s), v, false)]))
        end
    end

    for ct in cts
        connection2set!(connectionsets, namespace, ct, isouter)
    end

    # pre order traversal
    @set! sys.systems = map(s->generate_connection_set!(connectionsets, s, renamespace(namespace, nameof(s))), subsys)
    @set! sys.eqs = eqs
end

function Base.merge(csets::AbstractVector{<:ConnectionSet})
    mcsets = ConnectionSet[]
    ele2idx = Dict{ConnectionElement,Int}()
    cacheset = Set{ConnectionElement}()
    for cset in csets
        idx = nothing
        for e in cset.set
            idx = get(ele2idx, e, nothing)
            idx !== nothing && break
        end
        if idx === nothing
            push!(mcsets, cset)
            for e in cset.set
                ele2idx[e] = length(mcsets)
            end
        else
            for e in mcsets[idx].set
                push!(cacheset, e)
            end
            for e in cset.set
                push!(cacheset, e)
            end
            empty!(mcsets[idx].set)
            for e in cacheset
                ele2idx[e] = idx
                push!(mcsets[idx].set, e)
            end
            empty!(cacheset)
        end
    end
    mcsets
end

function generate_connection_equations_and_stream_connections(csets::AbstractVector{<:ConnectionSet})
    eqs = Equation[]
    stream_connections = ConnectionSet[]

    for cset in csets
        v = cset.set[1].v
        if hasmetadata(v, Symbolics.GetindexParent)
            v = getparent(v)
        end
        vtype = get_connection_type(v)
        if vtype === Stream
            push!(stream_connections, cset)
            continue
        elseif vtype === Flow
            rhs = 0
            for ele in cset.set
                v = namespaced_var(ele)
                rhs += ele.isouter ? -v : v
            end
            push!(eqs, 0 ~ rhs)
        else # Equality
            base = namespaced_var(cset.set[1])
            for i in 2:length(cset.set)
                v = namespaced_var(cset.set[i])
                push!(eqs, base ~ v)
            end
        end
    end
    eqs, stream_connections
end

function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10)
    sys, csets = generate_connection_set(sys)
    ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets)
    _sys = expand_instream(instream_csets, sys; debug=debug, tol=tol)
    sys = flatten(sys, true)
    @set! sys.eqs = [equations(_sys); ceqs]
end

function unnamespace(root, namespace)
    root === nothing && return namespace
    root = string(root)
    namespace = string(namespace)
    if length(namespace) > length(root)
        @assert root == namespace[1:length(root)]
        Symbol(namespace[nextind(namespace, length(root)):end])
    else
        @assert root == namespace
        nothing
    end
end

function expand_instream(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSystem, namespace=nothing, prevnamespace=nothing; debug=false, tol=1e-8)
    subsys = get_systems(sys)
    # post order traversal
    @set! sys.systems = map(s->expand_instream(csets, s, renamespace(namespace, nameof(s)), namespace; debug, tol), subsys)
    subsys = get_systems(sys)

    if debug
        @info "Expanding" namespace
    end

    sub = Dict()
    eqs = Equation[]
    instream_eqs = Equation[]
    instream_exprs = Set()
    for s in subsys
        for eq in get_eqs(s)
            eq = namespace_equation(eq, s)
            if collect_instream!(instream_exprs, eq)
                push!(instream_eqs, eq)
            else
                push!(eqs, eq)
            end
        end

    end

    for ex in instream_exprs
        cset, idx_in_set, sv = get_cset_sv(namespace, ex, csets)

        n_inners = n_outers = 0
        for (i, e) in enumerate(cset)
            if e.isouter
                n_outers += 1
            else
                n_inners += 1
            end
        end
        if debug
            @info "Expanding at [$idx_in_set]" ex ConnectionSet(cset)
            @show n_inners, n_outers
        end
        if n_inners == 1 && n_outers == 0
            sub[ex] = sv
        elseif n_inners == 2 && n_outers == 0
            other = idx_in_set == 1 ? 2 : 1
            sub[ex] = get_current_var(namespace, cset[other], sv)
        elseif n_inners == 1 && n_outers == 1
            if !cset[idx_in_set].isouter
                other = idx_in_set == 1 ? 2 : 1
                outerstream = get_current_var(namespace, cset[other], sv)
                sub[ex] = instream(outerstream)
            end
        else
            if !cset[idx_in_set].isouter
                fv = flowvar(first(cset).sys.sys)
                # mj.c.m_flow
                innerfvs = [get_current_var(namespace, s, fv) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter]
                innersvs = [get_current_var(namespace, s, sv) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter]
                # ck.m_flow
                outerfvs = [get_current_var(namespace, s, fv) for s in cset if s.isouter]
                outersvs = [get_current_var(namespace, s, sv) for s in cset if s.isouter]

                sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...)
            end
        end
    end

    # additional equations
    additional_eqs = Equation[]
    csets = filter(cset->any(e->e.sys.namespace === namespace, cset.set), csets)
    for cset′ in csets
        cset = cset′.set
        connectors = Vector{Any}(undef, length(cset))
        n_inners = n_outers = 0
        for (i, e) in enumerate(cset)
            connectors[i] = e.sys.sys
            if e.isouter
                n_outers += 1
            else
                n_inners += 1
            end
        end
        iszero(n_outers) && continue
        connector_representative = first(cset).sys.sys
        fv = flowvar(connector_representative)
        sv = first(cset).v
        vtype = get_connection_type(sv)
        vtype === Stream || continue
        if n_inners == 1 && n_outers == 1
            push!(additional_eqs, states(cset[1].sys.sys, sv) ~ states(cset[2].sys.sys, sv))
        elseif n_inners == 0 && n_outers == 2
            # we don't expand `instream` in this case.
            v1 = states(cset[1].sys.sys, sv)
            v2 = states(cset[2].sys.sys, sv)
            push!(additional_eqs, v1 ~ instream(v2))
            push!(additional_eqs, v2 ~ instream(v1))
        else
            sq = 0
            s_inners = (s for s in cset if !s.isouter)
            s_outers = (s for s in cset if s.isouter)
            for (q, oscq) in enumerate(s_outers)
                sq += sum(s->max(-states(s, fv), 0), s_inners)
                for (k, s) in enumerate(s_outers); k == q && continue
                    f = states(s.sys.sys, fv)
                    sq += max(f, 0)
                end

                num = 0
                den = 0
                for s in s_inners
                    f = states(s.sys.sys, fv)
                    tmp = positivemax(-f, sq; tol=tol)
                    den += tmp
                    num += tmp * states(s.sys.sys, sv)
                end
                for (k, s) in enumerate(s_outers); k == q && continue
                    f = states(s.sys.sys, fv)
                    tmp = positivemax(f, sq; tol=tol)
                    den += tmp
                    num += tmp * instream(states(s.sys.sys, sv))
                end
                push!(additional_eqs, states(oscq.sys.sys, sv) ~ num / den)
            end
        end
    end

    subed_eqs = substitute(instream_eqs, sub)
    if debug && !(isempty(csets) && isempty(additional_eqs) && isempty(instream_eqs))
        println("======================================")
        @info "Additional equations" csets
        display(additional_eqs)
        println("======================================")
        println("Substitutions")
        display(sub)
        println("======================================")
        println("Substituted equations")
        foreach(i->println(instream_eqs[i] => subed_eqs[i]), eachindex(subed_eqs))
        println("======================================")
    end

    @set! sys.systems = []
    @set! sys.eqs = [get_eqs(sys); eqs; subed_eqs; additional_eqs]
    sys
end

function get_current_var(namespace, cele, sv)
    states(renamespace(unnamespace(namespace, cele.sys.namespace), cele.sys.sys), sv)
end

function get_cset_sv(namespace, ex, csets)
    ns_sv = only(arguments(ex))
    full_name_sv = renamespace(namespace, ns_sv)

    cidx = -1
    idx_in_set = -1
    sv = ns_sv
    for (i, c) in enumerate(csets)
        crep = first(c.set)
        current = namespace == crep.sys.namespace
        for (j, v) in enumerate(c.set)
            if isequal(namespaced_var(v), full_name_sv) && (current || !v.isouter)
                cidx = i
                idx_in_set = j
                sv = v.v
            end
        end
    end
    cidx < 0 && error("$ns_sv is not a variable inside stream connectors")
    cset = csets[cidx].set
    #if namespace != first(cset).sys.namespace
    #    cset = map(c->@set(c.isouter = false), cset)
    #end
    cset, idx_in_set, sv
end

# instream runtime
@generated function _instream_split(::Val{inner_n}, ::Val{outer_n}, vars::NTuple{N,Any}) where {inner_n, outer_n, N}
    #instream_rt(innerfvs..., innersvs..., outerfvs..., outersvs...)
    ret = Expr(:tuple)
    # mj.c.m_flow
    inner_f = :(Base.@ntuple $inner_n i -> vars[i])
    offset = inner_n
    inner_s = :(Base.@ntuple $inner_n i -> vars[$offset+i])
    offset += inner_n
    # ck.m_flow
    outer_f = :(Base.@ntuple $outer_n i -> vars[$offset+i])
    offset += outer_n
    outer_s = :(Base.@ntuple $outer_n i -> vars[$offset+i])
    Expr(:tuple, inner_f, inner_s, outer_f, outer_s)
end

function instream_rt(ins::Val{inner_n}, outs::Val{outer_n}, vars::Vararg{Any,N}) where {inner_n, outer_n, N}
    @assert N == 2*(inner_n + outer_n)

    # inner: mj.c.m_flow
    # outer: ck.m_flow
    inner_f, inner_s, outer_f, outer_s = _instream_split(ins, outs, vars)

    T = float(first(inner_f))
    si = zero(T)
    num = den = zero(T)
    for f in inner_f
        si += max(-f, 0)
    end
    for f in outer_f
        si += max(f, 0)
    end
    #for (f, s) in zip(inner_f, inner_s)
    for j in 1:inner_n
        @inbounds f = inner_f[j]
        @inbounds s = inner_s[j]
        num += _positivemax(-f, si) * s
        den += _positivemax(-f, si)
    end
    #for (f, s) in zip(outer_f, outer_s)
    for j in 1:outer_n
        @inbounds f = outer_f[j]
        @inbounds s = outer_s[j]
        num += _positivemax(-f, si) * s
        den += _positivemax(-f, si)
    end
    return num / den
    #=
    si = sum(max(-mj.c.m_flow,0) for j in cat(1,1:i-1, i+1:N)) +
            sum(max(ck.m_flow ,0) for k  in 1:M)

    inStream(mi.c.h_outflow) =
       (sum(positiveMax(-mj.c.m_flow,si)*mj.c.h_outflow)
      +  sum(positiveMax(ck.m_flow,s_i)*inStream(ck.h_outflow)))/
     (sum(positiveMax(-mj.c.m_flow,s_i))
        +  sum(positiveMax(ck.m_flow,s_i)))
                  for j in 1:N and i <> j and mj.c.m_flow.min < 0,
                  for k in 1:M and ck.m_flow.max > 0
    =#
end
SymbolicUtils.promote_symtype(::typeof(instream_rt), ::Vararg) = Real


text:
import Pkg
Pkg.pkg"activate ."

Settings = []

#=
    Multi-processing case.
    Requires Distributed to be installed.
    - Uncomment below, specifying number of desired processes (default is number of system cores)
    - Comment `import Dance` below
=#
#=using Distributed
addprocs(length(Sys.cpu_info()))
@everywhere import Dance=#
import Dance

#=
    STEP 1
    If apart from `static_dir` as defined in `routes.jl` there are dirs to ignore adding to LOAD_PATH
        => add to `ignore_dirs`
    By default all other sub-dirs of project structure will be added
=#
ignore_dirs = [""]
if Dance.populate_settings(abspath(@__DIR__))
    global Settings
    Settings = Dance.Configuration.Settings
    for item in Dance.populate_load_path(abspath(@__DIR__); ignore_dirs=ignore_dirs)
        push!(LOAD_PATH, item)
    end
end

#=
    STEP 2
    Import modules here below, not above
    (STEP 1 populated `LOAD_PATH`)
=#
function main(args)
    if Dance.populate_router(abspath(@__DIR__))
        include(Settings[:routes_filename]*".jl")

        ## STEP 3: Add custom scripts here that need be run before launching Dance ##
        ##############################################

        start_server::Bool = true
        if length(args)>0
            if args[1]=="repl"
                start_server = false
            end
        end
        Dance.launch(start_server)
    end
end

main(ARGS)