text:
################################################################################
#
#  fq_poly.jl: Flint fq_poly (Polynomials over FqFiniteField)
#
################################################################################

export fq_poly, FqPolyRing

################################################################################
#
#  Type and parent object methods
#
################################################################################

parent_type(::Type{fq_poly}) = FqPolyRing

elem_type(::Type{FqPolyRing}) = fq_poly

base_ring(a::FqPolyRing) = a.base_ring

parent(a::fq_poly) = a.parent

var(a::FqPolyRing) = a.S

function check_parent(a::fq_poly, b::fq_poly)
   a.parent != b.parent &&
         error("Operations on distinct polynomial rings not supported")
end

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

length(x::fq_poly) = ccall((:fq_poly_length, :libflint), Int,
                                (Ref{fq_poly},), x)

function coeff(x::fq_poly, n::Int)
   n < 0 && throw(DomainError())
   F = (x.parent).base_ring
   temp = F(1)
   ccall((:fq_poly_get_coeff, :libflint), Void,
         (Ref{fq}, Ref{fq_poly}, Int, Ref{FqFiniteField}),
         temp, x, n, F)
   return temp
end

set_length!(x::fq_poly, n::Int) = ccall((:_fq_poly_set_length, :libflint), Void,
                              (Ref{fq_poly}, Int), x, n)

zero(a::FqPolyRing) = a(zero(base_ring(a)))

one(a::FqPolyRing) = a(one(base_ring(a)))

gen(a::FqPolyRing) = a([zero(base_ring(a)), one(base_ring(a))])

isgen(x::fq_poly) = ccall((:fq_poly_is_gen, :libflint), Bool,
                              (Ref{fq_poly}, Ref{FqFiniteField}),
                              x, base_ring(x.parent))

iszero(x::fq_poly) = ccall((:fq_poly_is_zero, :libflint), Bool,
                              (Ref{fq_poly}, Ref{FqFiniteField}),
                              x, base_ring(x.parent))

isone(x::fq_poly) = ccall((:fq_poly_is_one, :libflint), Bool,
                              (Ref{fq_poly}, Ref{FqFiniteField}),
                              x, base_ring(x.parent))

degree(f::fq_poly) = f.length - 1

function deepcopy_internal(a::fq_poly, dict::ObjectIdDict)
   z = fq_poly(a)
   z.parent = a.parent
   return z
end

################################################################################
#
#   Canonicalisation
#
################################################################################

canonical_unit(a::fq_poly) = canonical_unit(lead(a))

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::fq_poly)
   if length(x) == 0
      print(io, "0")
   else
      cstr = ccall((:fq_poly_get_str_pretty, :libflint), Ref{UInt8},
                  (Ref{fq_poly}, Ref{UInt8}, Ref{FqFiniteField}),
                  x, string(var(parent(x))),
                  ((x.parent).base_ring))
      print(io, unsafe_string(cstr))
      ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
   end
end

function show(io::IO, R::FqPolyRing)
   print(io, "Univariate Polynomial Ring in ")
   print(io, string(var(R)))
   print(io, " over ")
   show(io, base_ring(R))
end

show_minus_one(::Type{fq_poly}) = show_minus_one(fq)

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::fq_poly)
   z = parent(x)()
   ccall((:fq_poly_neg, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}),
         z, x, base_ring(parent(x)))
   return z
end

################################################################################
#
#  Binary operations
#
################################################################################

function +(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_add, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly},
         Ref{fq_poly}, Ref{FqFiniteField}),
         z, x, y, base_ring(parent(x)))
   return z
end

function -(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_sub, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly},
         Ref{fq_poly}, Ref{FqFiniteField}),
         z, x, y, base_ring(parent(x)))
   return z
end

function *(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_mul, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly},
         Ref{fq_poly}, Ref{FqFiniteField}),
         z, x, y, base_ring(parent(x)))
   return z
end

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

function *(x::fq, y::fq_poly)
   parent(x) != base_ring(parent(y)) &&
         error("Coefficient rings must be equal")
   z = parent(y)()
   ccall((:fq_poly_scalar_mul_fq, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly},
         Ref{fq}, Ref{FqFiniteField}),
         z, y, x, parent(x))
  return z
end

*(x::fq_poly, y::fq) = y*x

*(x::fmpz, y::fq_poly) = base_ring(parent(y))(x) * y

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

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

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

+(x::fq, y::fq_poly) = parent(y)(x) + y

+(x::fq_poly, y::fq) = y + x

+(x::fmpz, y::fq_poly) = base_ring(parent(y))(x) + y

+(x::fq_poly, y::fmpz) = y + x

+(x::fq_poly, y::Integer) = x + fmpz(y)

+(x::Integer, y::fq_poly) = y + x

-(x::fq, y::fq_poly) = parent(y)(x) - y

-(x::fq_poly, y::fq) = x - parent(x)(y)

-(x::fmpz, y::fq_poly) = base_ring(parent(y))(x) - y

-(x::fq_poly, y::fmpz) = x - base_ring(parent(x))(y)

-(x::fq_poly, y::Integer) = x - fmpz(y)

-(x::Integer, y::fq_poly) = fmpz(x) - y

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

function ^(x::fq_poly, y::Int)
   y < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fq_poly_pow, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Int, Ref{FqFiniteField}),
         z, x, y, base_ring(parent(x)))
   return z
end

################################################################################
#
#   Comparisons
#
################################################################################

function ==(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   r = ccall((:fq_poly_equal, :libflint), Cint,
             (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}),
             x, y, base_ring(parent(x)))
   return Bool(r)
end

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

function ==(x::fq_poly, y::fq)
   base_ring(parent(x)) != parent(y) && return false
   if length(x) > 1
      return false
   elseif length(x) == 1
      r = ccall((:fq_poly_equal_fq, :libflint), Cint,
                (Ref{fq_poly}, Ref{fq}, Ref{FqFiniteField}),
                x, y, base_ring(parent(x)))
      return Bool(r)
   else
      return iszero(y)
  end
end

==(x::fq, y::fq_poly) = y == x

==(x::fq_poly, y::fmpz) = x == base_ring(parent(x))(y)

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

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

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

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

function truncate(x::fq_poly, n::Int)
   n < 0 && throw(DomainError())
   if length(x) <= n
      return x
   end
   z = parent(x)()
   ccall((:fq_poly_set_trunc, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Int, Ref{FqFiniteField}),
         z, x, n, base_ring(parent(x)))
   return z
end

function mullow(x::fq_poly, y::fq_poly, n::Int)
   check_parent(x,y)
   n < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fq_poly_mullow, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Int, Ref{FqFiniteField}),
         z, x, y, n, base_ring(parent(x)))
   return z
end

################################################################################
#
#   Reversal
#
################################################################################

function reverse(x::fq_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fq_poly_reverse, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Int, Ref{FqFiniteField}),
         z, x, len, base_ring(parent(x)))
   return z
end

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

function shift_left(x::fq_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fq_poly_shift_left, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Int, Ref{FqFiniteField}),
         z, x, len, base_ring(parent(x)))
   return z
end

function shift_right(x::fq_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fq_poly_shift_right, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Int, Ref{FqFiniteField}),
         z, x, len, base_ring(parent(x)))
   return z
end

################################################################################
#
#   Euclidean division
#
################################################################################

function div(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_div_basecase, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
  return z
end

function rem(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_rem, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
  return z
end

mod(x::fq_poly, y::fq_poly) = rem(x, y)

function divrem(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   r = parent(x)()
   ccall((:fq_poly_divrem, :libflint), Void, (Ref{fq_poly},
         Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, r, x, y, base_ring(parent(x)))
   return z, r
end

################################################################################
#
#   Remove and valuation
#
################################################################################

doc"""
    remove(z::fq_poly, p::fq_poly)
> Computes the valuation of $z$ at $p$, that is, the largest $k$ such that
> $p^k$ divides $z$. Additionally, $z/p^k$ is returned as well.
>
> See also `valuation`, which only returns the valuation.
"""
function remove(z::fq_poly, p::fq_poly)
   check_parent(z,p)
   iszero(z) && error("Not yet implemented")
   z = deepcopy(z)
   v = ccall((:fq_poly_remove, :libflint), Int,
            (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}),
             z,  p, base_ring(parent(z)))
   return v, z
end

function divides(z::fq_poly, x::fq_poly)
   check_parent(z, x)
   q = parent(z)()
   v = Bool(ccall((:fq_poly_divides, :libflint), Cint,
            (Ref{fq_poly}, Ref{fq_poly},
             Ref{fq_poly}, Ref{FqFiniteField}),
             q, z, x, base_ring(parent(z))))
   return v, q
end

################################################################################
#
#   Modular arithmetic
#
################################################################################

function powmod(x::fq_poly, n::Int, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()

   if n < 0
      g, x = gcdinv(x, y)
      if g != 1
         error("Element not invertible")
      end
      n = -n
   end

   ccall((:fq_poly_powmod_ui_binexp, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Int, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, n, y, base_ring(parent(x)))
  return z
end

################################################################################
#
#   GCD
#
################################################################################

function gcd(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_gcd, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
   return z
end

function gcdinv(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   s = parent(x)()
   t = parent(x)()
   ccall((:fq_poly_xgcd, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, s, t, x, y, base_ring(parent(x)))
   return z, s
end

function gcdx(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   s = parent(x)()
   t = parent(x)()
   ccall((:fq_poly_xgcd, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, s, t, x, y, base_ring(parent(x)))
   return z, s, t
end

################################################################################
#
#   Evaluation
#
################################################################################

function evaluate(x::fq_poly, y::fq)
   base_ring(parent(x)) != parent(y) && error("Incompatible coefficient rings")
   z = parent(y)()
   ccall((:fq_poly_evaluate_fq, :libflint), Void,
         (Ref{fq}, Ref{fq_poly}, Ref{fq},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
   return z
end

################################################################################
#
#   Composition
#
################################################################################

function compose(x::fq_poly, y::fq_poly)
   check_parent(x,y)
   z = parent(x)()
   ccall((:fq_poly_compose, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
   return z
end

################################################################################
#
#   Derivative
#
################################################################################

function derivative(x::fq_poly)
   z = parent(x)()
   ccall((:fq_poly_derivative, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{FqFiniteField}),
         z, x, base_ring(parent(x)))
   return z
end

################################################################################
#
#  Inflation and deflation
#
################################################################################

function inflate(x::fq_poly, n::Int)
   z = parent(x)()
   ccall((:fq_poly_inflate, :libflint), Void, (Ref{fq_poly},
         Ref{fq_poly}, Culong, Ref{FqFiniteField}),
         z, x, UInt(n), base_ring(parent(x)))
   return z
end

function deflate(x::fq_poly, n::Int)
   z = parent(x)()
   ccall((:fq_poly_deflate, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Culong, Ref{FqFiniteField}),
         z, x, UInt(n), base_ring(parent(x)))
  return z
end

################################################################################
#
#  Irreducibility
#
################################################################################

doc"""
    isirreducible(x::fq_poly)
> Return `true` if $x$ is irreducible, otherwise return `false`.
"""
function isirreducible(x::fq_poly)
  return Bool(ccall((:fq_poly_is_irreducible, :libflint), Int32,
                    (Ref{fq_poly}, Ref{FqFiniteField} ),
                    x, base_ring(parent(x))))
end

################################################################################
#
#  Squarefree testing
#
################################################################################

doc"""
    issquarefree(x::fq_poly)
> Return `true` if $x$ is squarefree, otherwise return `false`.
"""
function issquarefree(x::fq_poly)
   return Bool(ccall((:fq_poly_is_squarefree, :libflint), Int32,
       (Ref{fq_poly}, Ref{FqFiniteField}), x, base_ring(parent(x))))
end

################################################################################
#
#  Factorization
#
################################################################################

doc"""
    factor(x::fq_poly)
> Return the factorisation of $x$.
"""
function factor(x::fq_poly)
   fac, z = _factor(x)
   return Fac(parent(x)(z), fac)
end

function _factor(x::fq_poly)
   R = parent(x)
   F = base_ring(R)
   a = F()
   fac = fq_poly_factor(F)
   ccall((:fq_poly_factor, :libflint), Void, (Ref{fq_poly_factor},
         Ref{fq}, Ref{fq_poly}, Ref{FqFiniteField}),
         fac, a, x, F)
   res = Dict{fq_poly,Int}()
   for i in 1:fac.num
      f = R()
      ccall((:fq_poly_factor_get_poly, :libflint), Void,
            (Ref{fq_poly}, Ref{fq_poly_factor}, Int,
            Ref{FqFiniteField}), f, fac, i-1, F)
      e = unsafe_load(fac.exp,i)
      res[f] = e
   end
   return res, a
end

doc"""
    factor_squarefree(x::fq_poly)
> Return the squarefree factorisation of $x$.
"""
function factor_squarefree(x::fq_poly)
  # _factor_squareefree does weird things if the polynomial is not monic
  return Fac(parent(x)(lead(x)), _factor_squarefree(divexact(x, lead(x))))
end

function _factor_squarefree(x::fq_poly)
  F = base_ring(parent(x))
  fac = fq_poly_factor(F)
  ccall((:fq_poly_factor_squarefree, :libflint), UInt,
        (Ref{fq_poly_factor}, Ref{fq_poly}, Ref{FqFiniteField}), fac, x, F)
  res = Dict{fq_poly,Int}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:fq_poly_factor_get_poly, :libflint), Void,
          (Ref{fq_poly}, Ref{fq_poly_factor}, Int,
          Ref{FqFiniteField}), f, fac, i-1, F)
    e = unsafe_load(fac.exp, i)
    res[f] = e
  end
  return res
end

doc"""
    factor_distinct_deg(x::fq_poly)
> Return the distinct degree factorisation of a squarefree polynomial $x$.
"""
function factor_distinct_deg(x::fq_poly)
   R = parent(x)
   F = base_ring(R)
   fac = fq_poly_factor(F)
   tmp = fq_poly_factor(F)
   ccall((:fq_poly_factor_fit_length, :libflint), Void,
         (Ref{fq_poly_factor}, Int, Ref{FqFiniteField}),
         tmp, degree(x), F)
   ccall((:fq_poly_factor_distinct_deg, :libflint), Void,
         (Ref{fq_poly_factor}, Ref{fq_poly}, Ref{Int},
         Ref{FqFiniteField}), fac, x, tmp.exp, F)
   res = Dict{Int, fq_poly}()
   for i in 1:fac.num
      f = R()
      ccall((:fq_poly_factor_get_poly, :libflint), Void,
            (Ref{fq_poly}, Ref{fq_poly_factor}, Int,
            Ref{FqFiniteField}), f, fac, i-1, F)
      d = unsafe_load(tmp.exp,i)
      res[d] = f
   end
   return res
end

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

function zero!(z::fq_poly)
   ccall((:fq_poly_zero, :libflint), Void,
         (Ref{fq_poly}, Ref{FqFiniteField}),
         z, base_ring(parent(z)))
   return z
end

function fit!(z::fq_poly, n::Int)
   ccall((:fq_poly_fit_length, :libflint), Void,
         (Ref{fq_poly}, Int, Ref{FqFiniteField}),
         z, n, base_ring(parent(z)))
   return nothing
end

function setcoeff!(z::fq_poly, n::Int, x::fq)
   ccall((:fq_poly_set_coeff, :libflint), Void,
         (Ref{fq_poly}, Int, Ref{fq}, Ref{FqFiniteField}),
         z, n, x, base_ring(parent(z)))
   return z
end

function mul!(z::fq_poly, x::fq_poly, y::fq_poly)
   ccall((:fq_poly_mul, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
   return z
end

function add!(z::fq_poly, x::fq_poly, y::fq_poly)
   ccall((:fq_poly_add, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
   return z
end

function sub!(z::fq_poly, x::fq_poly, y::fq_poly)
   ccall((:fq_poly_sub, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, x, y, base_ring(parent(x)))
   return z
end

function addeq!(z::fq_poly, x::fq_poly)
   ccall((:fq_poly_add, :libflint), Void,
         (Ref{fq_poly}, Ref{fq_poly}, Ref{fq_poly},
         Ref{FqFiniteField}), z, z, x, base_ring(parent(x)))
   return z
end

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

promote_rule(::Type{fq_poly}, ::Type{V}) where {V <: Integer} = fq_poly

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

promote_rule(::Type{fq_poly}, ::Type{fq}) = fq_poly

###############################################################################
#
#   Polynomial substitution
#
###############################################################################

function (f::fq_poly)(a::fq)
   if parent(a) != base_ring(f)
      return subst(f, a)
   end
   return evaluate(f, a)
end

################################################################################
#
#   Parent object call overloads
#
################################################################################

function (R::FqPolyRing)()
   z = fq_poly()
   z.parent = R
   return z
end

function (R::FqPolyRing)(x::fq)
  z = fq_poly(x)
  z.parent = R
  return z
end

function (R::FqPolyRing)(x::fmpz)
   return R(base_ring(R)(x))
end

function (R::FqPolyRing)(x::Integer)
   return R(fmpz(x))
end

function (R::FqPolyRing)(x::Array{fq, 1})
   length(x) == 0 && error("Array must be non-empty")
   base_ring(R) != parent(x[1]) && error("Coefficient rings must coincide")
   z = fq_poly(x)
   z.parent = R
   return z
end

function (R::FqPolyRing)(x::Array{fmpz, 1})
   length(x) == 0 && error("Array must be non-empty")
   z = fq_poly(x, base_ring(R))
   z.parent = R
   return z
end

function (R::FqPolyRing)(x::Array{T, 1}) where {T <: Integer}
   length(x) == 0 && error("Array must be non-empty")
   return R(map(fmpz, x))
end

function (R::FqPolyRing)(x::fmpz_poly)
   z = fq_poly(x, base_ring(R))
   z.parent = R
   return z
end

function (R::FqPolyRing)(x::fq_poly)
  parent(x) != R && error("Unable to coerce to polynomial")
  return x
end

################################################################################
#
#   PolynomialRing constructor
#
################################################################################

function PolynomialRing(R::FqFiniteField, s::AbstractString; cached = true)
   S = Symbol(s)
   parent_obj = FqPolyRing(R, S, cached)
   return parent_obj, parent_obj([R(0), R(1)])
end


text:
using SnpArrays
plinkname = "/scratch/users/bbchu/ukb_SHAPEIT/subset/ukb.10k.chr10"
xdata = SnpData(plinkname)
@time r, θ, α = fastphase(xdata; n=1000, out="ukb_chr10_n1000")

# for reading aggregated results later
# r, θ, α = process_fastphase_output(pwd(), T=1, extension="ukb_chr10_n1000")

ext: jl


text:
calcMidpoints(edges::AbstractVector) = Float64[0.5 * (edges[i] + edges[i+1]) for i in 1:length(edges)-1]

"Make histogram-like bins of data"
function binData(data, nbins)
  lo, hi = ignorenan_extrema(data)
  edges = collect(linspace(lo, hi, nbins+1))
  midpoints = calcMidpoints(edges)
  buckets = Int[max(2, min(searchsortedfirst(edges, x), length(edges)))-1 for x in data]
  counts = zeros(Int, length(midpoints))
  for b in buckets
    counts[b] += 1
  end
  edges, midpoints, buckets, counts
end

"""
A hacky replacement for a histogram when the backend doesn't support histograms directly.
Convert it into a bar chart with the appropriate x/y values.
"""
function histogramHack(; kw...)
  d = KW(kw)

  # we assume that the y kwarg is set with the data to be binned, and nbins is also defined
  edges, midpoints, buckets, counts = binData(d[:y], d[:bins])
  d[:x] = midpoints
  d[:y] = float(counts)
  d[:seriestype] = :bar
  d[:fillrange] = d[:fillrange] == nothing ? 0.0 : d[:fillrange]
  d
end

"""
A hacky replacement for a bar graph when the backend doesn't support bars directly.
Convert it into a line chart with fillrange set.
"""
function barHack(; kw...)
  d = KW(kw)
  midpoints = d[:x]
  heights = d[:y]
  fillrange = d[:fillrange] == nothing ? 0.0 : d[:fillrange]

  # estimate the edges
  dists = diff(midpoints) * 0.5
  edges = zeros(length(midpoints)+1)
  for i in 1:length(edges)
    if i == 1
      edge = midpoints[1] - dists[1]
    elseif i == length(edges)
      edge = midpoints[i-1] + dists[i-2]
    else
      edge = midpoints[i-1] + dists[i-1]
    end
    edges[i] = edge
  end

  x = Float64[]
  y = Float64[]
  for i in 1:length(heights)
    e1, e2 = edges[i:i+1]
    append!(x, [e1, e1, e2, e2])
    append!(y, [fillrange, heights[i], heights[i], fillrange])
  end

  d[:x] = x
  d[:y] = y
  d[:seriestype] = :path
  d[:fillrange] = fillrange
  d
end

"""
A hacky replacement for a sticks graph when the backend doesn't support sticks directly.
Convert it into a line chart that traces the sticks, and a scatter that sets markers at the points.
"""
function sticksHack(; kw...)
  dLine = KW(kw)
  dScatter = copy(dLine)

  # these are the line vertices
  x = Float64[]
  y = Float64[]
  fillrange = dLine[:fillrange] == nothing ? 0.0 : dLine[:fillrange]

  # calculate the vertices
  yScatter = dScatter[:y]
  for (i,xi) in enumerate(dScatter[:x])
    yi = yScatter[i]
    for j in 1:3 push!(x, xi) end
    append!(y, [fillrange, yScatter[i], fillrange])
  end

  # change the line args
  dLine[:x] = x
  dLine[:y] = y
  dLine[:seriestype] = :path
  dLine[:markershape] = :none
  dLine[:fillrange] = nothing

  # change the scatter args
  dScatter[:seriestype] = :none

  dLine, dScatter
end

function regressionXY(x, y)
  # regress
  β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)

  # make a line segment
  regx = [ignorenan_minimum(x), ignorenan_maximum(x)]
  regy = β * regx + α
  regx, regy
end

function replace_image_with_heatmap(z::Array{T}) where T<:Colorant
    @show T, size(z)
    n, m = size(z)
    # idx = 0
    colors = ColorGradient(vec(z))
    newz = reshape(linspace(0, 1, n*m), n, m)
    newz, colors
    # newz = zeros(n, m)
    # for i=1:n, j=1:m
    #     push!(colors, T(z[i,j]...))
    #     newz[i,j] = idx / (n*m-1)
    #     idx += 1
    # end
    # newz, ColorGradient(colors)
end

function imageHack(d::KW)
    is_seriestype_supported(:heatmap) || error("Neither :image or :heatmap are supported!")
    d[:seriestype] = :heatmap
    d[:z], d[:fillcolor] = replace_image_with_heatmap(d[:z].surf)
end
# ---------------------------------------------------------------

"Build line segments for plotting"
mutable struct Segments{T}
    pts::Vector{T}
end

# Segments() = Segments{Float64}(zeros(0))

Segments() = Segments(Float64)
Segments(::Type{T}) where {T} = Segments(T[])
Segments(p::Int) = Segments(NTuple{2,Float64}[])

# Segments() = Segments(zeros(0))

to_nan(::Type{Float64}) = NaN
to_nan(::Type{NTuple{2,Float64}}) = (NaN, NaN)

coords(segs::Segments{Float64}) = segs.pts
coords(segs::Segments{NTuple{2,Float64}}) = Float64[p[1] for p in segs.pts], Float64[p[2] for p in segs.pts]

function Base.push!(segments::Segments{T}, vs...) where T
    if !isempty(segments.pts)
        push!(segments.pts, to_nan(T))
    end
    for v in vs
        push!(segments.pts, convert(T,v))
    end
    segments
end

function Base.push!(segments::Segments{T}, vs::AVec) where T
    if !isempty(segments.pts)
        push!(segments.pts, to_nan(T))
    end
    for v in vs
        push!(segments.pts, convert(T,v))
    end
    segments
end

# -----------------------------------------------------
# helper to manage NaN-separated segments

mutable struct SegmentsIterator
    args::Tuple
    n::Int
end

function iter_segments(args...)
    tup = Plots.wraptuple(args)
    n = maximum(map(length, tup))
    SegmentsIterator(tup, n)
end

# helpers to figure out if there are NaN values in a list of array types
anynan(i::Int, args::Tuple) = any(a -> !isfinite(_cycle(a,i)), args)
anynan(istart::Int, iend::Int, args::Tuple) = any(i -> anynan(i, args), istart:iend)
allnan(istart::Int, iend::Int, args::Tuple) = all(i -> anynan(i, args), istart:iend)

function Base.start(itr::SegmentsIterator)
    nextidx = 1
    if anynan(1, itr.args)
        _, nextidx = next(itr, 1)
    end
    nextidx
end
Base.done(itr::SegmentsIterator, nextidx::Int) = nextidx > itr.n
function Base.next(itr::SegmentsIterator, nextidx::Int)
    i = istart = iend = nextidx

    # find the next NaN, and iend is the one before
    while i <= itr.n + 1
        if i > itr.n || anynan(i, itr.args)
            # done... array end or found NaN
            iend = i-1
            break
        end
        i += 1
    end

    # find the next non-NaN, and set nextidx
    while i <= itr.n
        if !anynan(i, itr.args)
            break
        end
        i += 1
    end

    istart:iend, i
end

# Find minimal type that can contain NaN and x
# To allow use of NaN separated segments with categorical x axis

float_extended_type(x::AbstractArray{T}) where {T} = Union{T,Float64}
float_extended_type(x::AbstractArray{T}) where {T<:Real} = Float64

# ------------------------------------------------------------------------------------

nop() = nothing
notimpl() = error("This has not been implemented yet")

isnothing(x::Void) = true
isnothing(x) = false

_cycle(wrapper::InputWrapper, idx::Int) = wrapper.obj
_cycle(wrapper::InputWrapper, idx::AVec{Int}) = wrapper.obj

_cycle(v::AVec, idx::Int) = v[mod1(idx, length(v))]
_cycle(v::AMat, idx::Int) = size(v,1) == 1 ? v[1, mod1(idx, size(v,2))] : v[:, mod1(idx, size(v,2))]
_cycle(v, idx::Int)       = v

_cycle(v::AVec, indices::AVec{Int}) = map(i -> _cycle(v,i), indices)
_cycle(v::AMat, indices::AVec{Int}) = map(i -> _cycle(v,i), indices)
_cycle(v, indices::AVec{Int})       = fill(v, length(indices))

_cycle(grad::ColorGradient, idx::Int) = _cycle(grad.colors, idx)
_cycle(grad::ColorGradient, indices::AVec{Int}) = _cycle(grad.colors, indices)

makevec(v::AVec) = v
makevec(v::T) where {T} = T[v]

"duplicate a single value, or pass the 2-tuple through"
maketuple(x::Real)                     = (x,x)
maketuple(x::Tuple{T,S}) where {T,S} = x

mapFuncOrFuncs(f::Function, u::AVec)        = map(f, u)
mapFuncOrFuncs(fs::AVec{F}, u::AVec) where {F<:Function} = [map(f, u) for f in fs]

unzip(xy::AVec{Tuple{X,Y}}) where {X,Y}              = [t[1] for t in xy], [t[2] for t in xy]
unzip(xyz::AVec{Tuple{X,Y,Z}}) where {X,Y,Z}         = [t[1] for t in xyz], [t[2] for t in xyz], [t[3] for t in xyz]
unzip(xyuv::AVec{Tuple{X,Y,U,V}}) where {X,Y,U,V}    = [t[1] for t in xyuv], [t[2] for t in xyuv], [t[3] for t in xyuv], [t[4] for t in xyuv]

unzip(xy::AVec{FixedSizeArrays.Vec{2,T}}) where {T}  = T[t[1] for t in xy], T[t[2] for t in xy]
unzip(xy::FixedSizeArrays.Vec{2,T}) where {T}        = T[xy[1]], T[xy[2]]

unzip(xyz::AVec{FixedSizeArrays.Vec{3,T}}) where {T} = T[t[1] for t in xyz], T[t[2] for t in xyz], T[t[3] for t in xyz]
unzip(xyz::FixedSizeArrays.Vec{3,T}) where {T}       = T[xyz[1]], T[xyz[2]], T[xyz[3]]

unzip(xyuv::AVec{FixedSizeArrays.Vec{4,T}}) where {T} = T[t[1] for t in xyuv], T[t[2] for t in xyuv], T[t[3] for t in xyuv], T[t[4] for t in xyuv]
unzip(xyuv::FixedSizeArrays.Vec{4,T}) where {T}       = T[xyuv[1]], T[xyuv[2]], T[xyuv[3]], T[xyuv[4]]

# given 2-element lims and a vector of data x, widen lims to account for the extrema of x
function _expand_limits(lims, x)
  try
    e1, e2 = ignorenan_extrema(x)
    lims[1] = NaNMath.min(lims[1], e1)
    lims[2] = NaNMath.max(lims[2], e2)
  # catch err
  #   warn(err)
  end
  nothing
end

expand_data(v, n::Integer) = [_cycle(v, i) for i=1:n]

# if the type exists in a list, replace the first occurence.  otherwise add it to the end
function addOrReplace(v::AbstractVector, t::DataType, args...; kw...)
    for (i,vi) in enumerate(v)
        if isa(vi, t)
            v[i] = t(args...; kw...)
            return
        end
    end
    push!(v, t(args...; kw...))
    return
end

function replaceType(vec, val)
  filter!(x -> !isa(x, typeof(val)), vec)
  push!(vec, val)
end

function replaceAlias!(d::KW, k::Symbol, aliases::Dict{Symbol,Symbol})
  if haskey(aliases, k)
    d[aliases[k]] = pop!(d, k)
  end
end

function replaceAliases!(d::KW, aliases::Dict{Symbol,Symbol})
  ks = collect(keys(d))
  for k in ks
      replaceAlias!(d, k, aliases)
  end
end

createSegments(z) = collect(repmat(reshape(z,1,:),2,1))[2:end]

Base.first(c::Colorant) = c
Base.first(x::Symbol) = x

sortedkeys(d::Dict) = sort(collect(keys(d)))

const _scale_base = Dict{Symbol, Real}(
    :log10 => 10,
    :log2 => 2,
    :ln => e,
)

"create an (n+1) list of the outsides of heatmap rectangles"
function heatmap_edges(v::AVec, scale::Symbol = :identity)
  vmin, vmax = ignorenan_extrema(v)
  extra_min = extra_max = 0.5 * (vmax-vmin) / (length(v)-1)
  if scale in _logScales
      vmin > 0 || error("The axis values must be positive for a $scale scale")
      while vmin - extra_min <= 0
          extra_min /= _scale_base[scale]
      end
  end
  vcat(vmin-extra_min, 0.5 * (v[1:end-1] + v[2:end]), vmax+extra_max)
end

function calc_r_extrema(x, y)
    xmin, xmax = ignorenan_extrema(x)
    ymin, ymax = ignorenan_extrema(y)
    r = 0.5 * NaNMath.min(xmax - xmin, ymax - ymin)
    ignorenan_extrema(r)
end

function convert_to_polar(x, y, r_extrema = calc_r_extrema(x, y))
    rmin, rmax = r_extrema
    theta, r = filter_radial_data(x, y, r_extrema)
    r = (r - rmin) / (rmax - rmin)
    x = r.*cos.(theta)
    y = r.*sin.(theta)
    x, y
end

# Filters radial data for points within the axis limits
function filter_radial_data(theta, r, r_extrema::Tuple{Real, Real})
    n = max(length(theta), length(r))
    rmin, rmax = r_extrema
    x, y = zeros(n), zeros(n)
    for i in 1:n
        x[i] = _cycle(theta, i)
        y[i] = _cycle(r, i)
    end
    points = map((a, b) -> (a, b), x, y)
    filter!(a -> a[2] >= rmin && a[2] <= rmax, points)
    x = map(a -> a[1], points)
    y = map(a -> a[2], points)
    x, y
end

function fakedata(sz...)
  y = zeros(sz...)
  for r in 2:size(y,1)
    y[r,:] = 0.95 * vec(y[r-1,:]) + randn(size(y,2))
  end
  y
end

isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited
isatom() = isdefined(Main, :Atom) && Main.Atom.isconnected()

function is_installed(pkgstr::AbstractString)
    try
        Pkg.installed(pkgstr) === nothing ? false : true
    catch
        false
    end
end

istuple(::Tuple) = true
istuple(::Any)   = false
isvector(::AVec) = true
isvector(::Any)  = false
ismatrix(::AMat) = true
ismatrix(::Any)  = false
isscalar(::Real) = true
isscalar(::Any)  = false

is_2tuple(v) = typeof(v) <: Tuple && length(v) == 2

isvertical(d::KW) = get(d, :orientation, :vertical) in (:vertical, :v, :vert)
isvertical(series::Series) = isvertical(series.d)

ticksType(ticks::AVec{T}) where {T<:Real}                      = :ticks
ticksType(ticks::AVec{T}) where {T<:AbstractString}            = :labels
ticksType(ticks::Tuple{T,S}) where {T<:AVec,S<:AVec}  = :ticks_and_labels
ticksType(ticks)                                        = :invalid

limsType(lims::Tuple{T,S}) where {T<:Real,S<:Real}    = :limits
limsType(lims::Symbol)                                  = lims == :auto ? :auto : :invalid
limsType(lims)                                          = :invalid

# axis_Symbol(letter, postfix) = Symbol(letter * postfix)
# axis_symbols(letter, postfix...) = map(s -> axis_Symbol(letter, s), postfix)

Base.convert(::Type{Vector{T}}, rng::Range{T}) where {T<:Real}         = T[x for x in rng]
Base.convert(::Type{Vector{T}}, rng::Range{S}) where {T<:Real,S<:Real} = T[x for x in rng]

Base.merge(a::AbstractVector, b::AbstractVector) = sort(unique(vcat(a,b)))

nanpush!(a::AbstractVector, b) = (push!(a, NaN); push!(a, b))
nanappend!(a::AbstractVector, b) = (push!(a, NaN); append!(a, b))

function nansplit(v::AVec)
    vs = Vector{eltype(v)}[]
    while true
        idx = findfirst(isnan, v)
        if idx <= 0
            # no nans
            push!(vs, v)
            break
        elseif idx > 1
            push!(vs, v[1:idx-1])
        end
        v = v[idx+1:end]
    end
    vs
end

function nanvcat(vs::AVec)
    v_out = zeros(0)
    for v in vs
        nanappend!(v_out, v)
    end
    v_out
end

# given an array of discrete values, turn it into an array of indices of the unique values
# returns the array of indices (znew) and a vector of unique values (vals)
function indices_and_unique_values(z::AbstractArray)
    vals = sort(unique(z))
    vmap = Dict([(v,i) for (i,v) in enumerate(vals)])
    newz = map(zi -> vmap[zi], z)
    newz, vals
end

# this is a helper function to determine whether we need to transpose a surface matrix.
# it depends on whether the backend matches rows to x (transpose_on_match == true) or vice versa
# for example: PyPlot sends rows to y, so transpose_on_match should be true
function transpose_z(d, z, transpose_on_match::Bool = true)
    if d[:match_dimensions] == transpose_on_match
        # z'
        permutedims(z, [2,1])
    else
        z
    end
end

function ok(x::Number, y::Number, z::Number = 0)
    isfinite(x) && isfinite(y) && isfinite(z)
end
ok(tup::Tuple) = ok(tup...)

# compute one side of a fill range from a ribbon
function make_fillrange_side(y, rib)
    frs = zeros(length(y))
    for (i, (yi, ri)) in enumerate(zip(y, Base.Iterators.cycle(rib)))
        frs[i] = yi + ri
    end
    frs
end

# turn a ribbon into a fillrange
function make_fillrange_from_ribbon(kw::KW)
    y, rib = kw[:y], kw[:ribbon]
    rib = wraptuple(rib)
    rib1, rib2 = -first(rib), last(rib)
    # kw[:ribbon] = nothing
    kw[:fillrange] = make_fillrange_side(y, rib1), make_fillrange_side(y, rib2)
    (get(kw, :fillalpha, nothing) == nothing) && (kw[:fillalpha] = 0.5)
end

#turn tuple of fillranges to one path
function concatenate_fillrange(x,y::Tuple)
    rib1, rib2 = first(y), last(y)
    yline = vcat(rib1,(rib2)[end:-1:1])
    xline = vcat(x,x[end:-1:1])
    return xline, yline
end

function get_sp_lims(sp::Subplot, letter::Symbol)
    axis_limits(sp[Symbol(letter, :axis)])
end

"""
    xlims([plt])

Returns the x axis limits of the current plot or subplot
"""
xlims(sp::Subplot) = get_sp_lims(sp, :x)

"""
    ylims([plt])

Returns the y axis limits of the current plot or subplot
"""
ylims(sp::Subplot) = get_sp_lims(sp, :y)

"""
    zlims([plt])

Returns the z axis limits of the current plot or subplot
"""
zlims(sp::Subplot) = get_sp_lims(sp, :z)

xlims(plt::Plot, sp_idx::Int = 1) = xlims(plt[sp_idx])
ylims(plt::Plot, sp_idx::Int = 1) = ylims(plt[sp_idx])
zlims(plt::Plot, sp_idx::Int = 1) = zlims(plt[sp_idx])
xlims(sp_idx::Int = 1) = xlims(current(), sp_idx)
ylims(sp_idx::Int = 1) = ylims(current(), sp_idx)
zlims(sp_idx::Int = 1) = zlims(current(), sp_idx)

function get_clims(sp::Subplot)
    zmin, zmax = Inf, -Inf
    for series in series_list(sp)
        for vals in (series[:z], series[:line_z], series[:marker_z], series[:fill_z])
            if (typeof(vals) <: AbstractSurface) && (eltype(vals.surf) <: Real)
                zmin, zmax = _update_clims(zmin, zmax, ignorenan_extrema(vals.surf)...)
            elseif (vals != nothing) && (eltype(vals) <: Real)
                zmin, zmax = _update_clims(zmin, zmax, ignorenan_extrema(vals)...)
            end
        end
    end
    clims = sp[:clims]
    if is_2tuple(clims)
        isfinite(clims[1]) && (zmin = clims[1])
        isfinite(clims[2]) && (zmax = clims[2])
    end
    return zmin < zmax ? (zmin, zmax) : (0.0, 0.0)
end

_update_clims(zmin, zmax, emin, emax) = min(zmin, emin), max(zmax, emax)

function hascolorbar(series::Series)
    st = series[:seriestype]
    hascbar = st in (:heatmap, :contour)
    if series[:marker_z] != nothing || series[:line_z] != nothing || series[:fill_z] != nothing
        hascbar = true
    end
    # no colorbar if we are creating a surface LightSource
    if xor(st == :surface, series[:fill_z] != nothing)
        hascbar = true
    end
    return hascbar
end

function hascolorbar(sp::Subplot)
    cbar = sp[:colorbar]
    hascbar = false
    if cbar != :none
        for series in series_list(sp)
            if hascolorbar(series)
                hascbar = true
            end
        end
    end
    hascbar
end

function get_linecolor(sp::Subplot, series::Series, i::Int = 1)
    lc = series[:linecolor]
    lz = series[:line_z]
    if lz == nothing
        isa(lc, ColorGradient) ? lc : _cycle(lc, i)
    else
        cmin, cmax = get_clims(sp)
        grad = isa(lc, ColorGradient) ? lc : cgrad()
        grad[clamp((_cycle(lz, i) - cmin) / (cmax - cmin), 0, 1)]
    end
end

function get_fillcolor(sp::Subplot, series::Series, i::Int = 1)
    fc = series[:fillcolor]
    fz = series[:fill_z]
    lz = series[:line_z]
    if fz == nothing && lz == nothing
        isa(fc, ColorGradient) ? fc : _cycle(fc, i)
    else
        cmin, cmax = get_clims(sp)
        grad = isa(fc, ColorGradient) ? fc : cgrad()
        if fz != nothing
            grad[clamp((_cycle(fz, i) - cmin) / (cmax - cmin), 0, 1)]
        elseif lz != nothing
            grad[clamp((_cycle(lz, i) - cmin) / (cmax - cmin), 0, 1)]
        end
    end
end

# ---------------------------------------------------------------

makekw(; kw...) = KW(kw)

wraptuple(x::Tuple) = x
wraptuple(x) = (x,)

trueOrAllTrue(f::Function, x::AbstractArray) = all(f, x)
trueOrAllTrue(f::Function, x) = f(x)

allLineTypes(arg)   = trueOrAllTrue(a -> get(_typeAliases, a, a) in _allTypes, arg)
allStyles(arg)      = trueOrAllTrue(a -> get(_styleAliases, a, a) in _allStyles, arg)
allShapes(arg)      = trueOrAllTrue(a -> is_marker_supported(get(_markerAliases, a, a)), arg) ||
                        trueOrAllTrue(a -> isa(a, Shape), arg)
allAlphas(arg)      = trueOrAllTrue(a -> (typeof(a) <: Real && a > 0 && a < 1) ||
                        (typeof(a) <: AbstractFloat && (a == zero(typeof(a)) || a == one(typeof(a)))), arg)
allReals(arg)       = trueOrAllTrue(a -> typeof(a) <: Real, arg)
allFunctions(arg)   = trueOrAllTrue(a -> isa(a, Function), arg)

# ---------------------------------------------------------------
# ---------------------------------------------------------------

"""
Allows temporary setting of backend and defaults for Plots. Settings apply only for the `do` block.  Example:
```
with(:gr, size=(400,400), type=:histogram) do
  plot(rand(10))
  plot(rand(10))
end
```
"""
function with(f::Function, args...; kw...)
    newdefs = KW(kw)

    if :canvas in args
        newdefs[:xticks] = nothing
        newdefs[:yticks] = nothing
        newdefs[:grid] = false
        newdefs[:legend] = false
    end

  # dict to store old and new keyword args for anything that changes
  olddefs = KW()
  for k in keys(newdefs)
    olddefs[k] = default(k)
  end

  # save the backend
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end
  oldbackend = CURRENT_BACKEND.sym

  for arg in args

    # change backend?
    if arg in backends()
      backend(arg)
    end

    # # TODO: generalize this strategy to allow args as much as possible
    # #       as in:  with(:gadfly, :scatter, :legend, :grid) do; ...; end
    # # TODO: can we generalize this enough to also do something similar in the plot commands??

    # k = :seriestype
    # if arg in _allTypes
    #   olddefs[k] = default(k)
    #   newdefs[k] = arg
    # elseif haskey(_typeAliases, arg)
    #   olddefs[k] = default(k)
    #   newdefs[k] = _typeAliases[arg]
    # end

    k = :legend
    if arg in (k, :leg)
      olddefs[k] = default(k)
      newdefs[k] = true
    end

    k = :grid
    if arg == k
      olddefs[k] = default(k)
      newdefs[k] = true
    end
  end

  # display(olddefs)
  # display(newdefs)

  # now set all those defaults
  default(; newdefs...)

  # call the function
  ret = f()

  # put the defaults back
  default(; olddefs...)

  # revert the backend
  if CURRENT_BACKEND.sym != oldbackend
    backend(oldbackend)
  end

  # return the result of the function
  ret
end

# ---------------------------------------------------------------
# ---------------------------------------------------------------

mutable struct DebugMode
  on::Bool
end
const _debugMode = DebugMode(false)

function debugplots(on = true)
  _debugMode.on = on
end

debugshow(x) = show(x)
debugshow(x::AbstractArray) = print(summary(x))

function dumpdict(d::KW, prefix = "", alwaysshow = false)
  _debugMode.on || alwaysshow || return
  println()
  if prefix != ""
    println(prefix, ":")
  end
  for k in sort(collect(keys(d)))
    @printf("%14s: ", k)
    debugshow(d[k])
    println()
  end
  println()
end
DD(d::KW, prefix = "") = dumpdict(d, prefix, true)

function dumpcallstack()
  error()  # well... you wanted the stacktrace, didn't you?!?
end

# ---------------------------------------------------------------
# ---------------------------------------------------------------
# used in updating an existing series

extendSeriesByOne(v::UnitRange{Int}, n::Int = 1) = isempty(v) ? (1:n) : (minimum(v):maximum(v)+n)
extendSeriesByOne(v::AVec, n::Integer = 1)       = isempty(v) ? (1:n) : vcat(v, (1:n) + ignorenan_maximum(v))
extendSeriesData(v::Range{T}, z::Real) where {T}        = extendSeriesData(float(collect(v)), z)
extendSeriesData(v::Range{T}, z::AVec) where {T}        = extendSeriesData(float(collect(v)), z)
extendSeriesData(v::AVec{T}, z::Real) where {T}         = (push!(v, convert(T, z)); v)
extendSeriesData(v::AVec{T}, z::AVec) where {T}         = (append!(v, convert(Vector{T}, z)); v)

# -------------------------------------------------------
# NOTE: backends should implement the following methods to get/set the x/y/z data objects

tovec(v::AbstractVector) = v
tovec(v::Void) = zeros(0)

function getxy(plt::Plot, i::Integer)
    d = plt.series_list[i].d
    tovec(d[:x]), tovec(d[:y])
end
function getxyz(plt::Plot, i::Integer)
    d = plt.series_list[i].d
    tovec(d[:x]), tovec(d[:y]), tovec(d[:z])
end

function setxy!(plt::Plot, xy::Tuple{X,Y}, i::Integer) where {X,Y}
    series = plt.series_list[i]
    series.d[:x], series.d[:y] = xy
    sp = series.d[:subplot]
    reset_extrema!(sp)
    _series_updated(plt, series)
end
function setxyz!(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer) where {X,Y,Z}
    series = plt.series_list[i]
    series.d[:x], series.d[:y], series.d[:z] = xyz
    sp = series.d[:subplot]
    reset_extrema!(sp)
    _series_updated(plt, series)
end

function setxyz!(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer) where {X,Y,Z<:AbstractMatrix}
    setxyz!(plt, (xyz[1], xyz[2], Surface(xyz[3])), i)
end

# -------------------------------------------------------
# indexing notation

# Base.getindex(plt::Plot, i::Integer) = getxy(plt, i)
Base.setindex!(plt::Plot, xy::Tuple{X,Y}, i::Integer) where {X,Y} = (setxy!(plt, xy, i); plt)
Base.setindex!(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer) where {X,Y,Z} = (setxyz!(plt, xyz, i); plt)

# -------------------------------------------------------

# operate on individual series

function push_x!(series::Series, xi)
    push!(series[:x], xi)
    expand_extrema!(series[:subplot][:xaxis], xi)
    return
end
function push_y!(series::Series, yi)
    push!(series[:y], yi)
    expand_extrema!(series[:subplot][:yaxis], yi)
    return
end
function push_z!(series::Series, zi)
    push!(series[:z], zi)
    expand_extrema!(series[:subplot][:zaxis], zi)
    return
end

function Base.push!(series::Series, yi)
    x = extendSeriesByOne(series[:x])
    expand_extrema!(series[:subplot][:xaxis], x[end])
    series[:x] = x
    push_y!(series, yi)
end
Base.push!(series::Series, xi, yi) = (push_x!(series,xi); push_y!(series,yi))
Base.push!(series::Series, xi, yi, zi) = (push_x!(series,xi); push_y!(series,yi); push_z!(series,zi))

# -------------------------------------------------------

function attr!(series::Series; kw...)
    d = KW(kw)
    preprocessArgs!(d)
    for (k,v) in d
        if haskey(_series_defaults, k)
            series[k] = v
        else
            warn("unused key $k in series attr")
        end
    end
    _series_updated(series[:subplot].plt, series)
    series
end

function attr!(sp::Subplot; kw...)
    d = KW(kw)
    preprocessArgs!(d)
    for (k,v) in d
        if haskey(_subplot_defaults, k)
            sp[k] = v
        else
            warn("unused key $k in subplot attr")
        end
    end
    sp
end

# -------------------------------------------------------
# push/append for one series

# push value to first series
Base.push!(plt::Plot, y::Real) = push!(plt, 1, y)
Base.push!(plt::Plot, x::Real, y::Real) = push!(plt, 1, x, y)
Base.push!(plt::Plot, x::Real, y::Real, z::Real) = push!(plt, 1, x, y, z)

# y only
function Base.push!(plt::Plot, i::Integer, y::Real)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesByOne(xdata), extendSeriesData(ydata, y)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, y::AVec)
    xdata, ydata = plt[i]
    if !isa(xdata, UnitRange{Int})
        error("Expected x is a UnitRange since you're trying to push a y value only")
    end
    plt[i] = (extendSeriesByOne(xdata, length(y)), extendSeriesData(ydata, y))
    plt
end

# x and y
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec)
    @assert length(x) == length(y)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y)), i)
    plt
end

# x, y, and z
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real, z::Real)
    # @show i, x, y, z
    xdata, ydata, zdata = getxyz(plt, i)
    # @show xdata, ydata, zdata
    setxyz!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y), extendSeriesData(zdata, z)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec, z::AVec)
    @assert length(x) == length(y) == length(z)
    xdata, ydata, zdata = getxyz(plt, i)
    setxyz!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y), extendSeriesData(zdata, z)), i)
    plt
end

# tuples
Base.push!(plt::Plot, xy::Tuple{X,Y}) where {X,Y}                  = push!(plt, 1, xy...)
Base.push!(plt::Plot, xyz::Tuple{X,Y,Z}) where {X,Y,Z}             = push!(plt, 1, xyz...)
Base.push!(plt::Plot, i::Integer, xy::Tuple{X,Y}) where {X,Y}      = push!(plt, i, xy...)
Base.push!(plt::Plot, i::Integer, xyz::Tuple{X,Y,Z}) where {X,Y,Z} = push!(plt, i, xyz...)

# -------------------------------------------------------
# push/append for all series

# push y[i] to the ith series
function Base.push!(plt::Plot, y::AVec)
    ny = length(y)
    for i in 1:plt.n
        push!(plt, i, y[mod1(i,ny)])
    end
    plt
end

# push y[i] to the ith series
# same x for each series
function Base.push!(plt::Plot, x::Real, y::AVec)
    push!(plt, [x], y)
end

# push (x[i], y[i]) to the ith series
function Base.push!(plt::Plot, x::AVec, y::AVec)
    nx = length(x)
    ny = length(y)
    for i in 1:plt.n
        push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)])
    end
    plt
end

# push (x[i], y[i], z[i]) to the ith series
function Base.push!(plt::Plot, x::AVec, y::AVec, z::AVec)
    nx = length(x)
    ny = length(y)
    nz = length(z)
    for i in 1:plt.n
        push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)], z[mod1(i,nz)])
    end
    plt
end

# ---------------------------------------------------------------

# Some conversion functions
# note: I borrowed these conversion constants from Compose.jl's Measure

const PX_PER_INCH   = 100
const DPI           = PX_PER_INCH
const MM_PER_INCH   = 25.4
const MM_PER_PX     = MM_PER_INCH / PX_PER_INCH

inch2px(inches::Real)   = float(inches * PX_PER_INCH)
px2inch(px::Real)       = float(px / PX_PER_INCH)
inch2mm(inches::Real)   = float(inches * MM_PER_INCH)
mm2inch(mm::Real)       = float(mm / MM_PER_INCH)
px2mm(px::Real)         = float(px * MM_PER_PX)
mm2px(mm::Real)         = float(px / MM_PER_PX)

"Smallest x in plot"
xmin(plt::Plot) = ignorenan_minimum([ignorenan_minimum(series.d[:x]) for series in plt.series_list])
"Largest x in plot"
xmax(plt::Plot) = ignorenan_maximum([ignorenan_maximum(series.d[:x]) for series in plt.series_list])

"Extrema of x-values in plot"
ignorenan_extrema(plt::Plot) = (xmin(plt), xmax(plt))


text:
using Test
using SparseFFT
using FFTW, Random

include("spfft1.jl")
include("spfft2.jl")


text:
module T4NAFEMS_examples
using FinEtools
using FinEtools.AlgoBaseModule.richextrapol

function T4NAFEMS_T3_algo()
    ## Two-dimensional heat transfer with convection: convergence study
    #

    ## Description
    #
    # Consider a plate of uniform thickness, measuring 0.6 m by 1.0 m. On one
    # short edge the temperature is fixed at 100 °C, and on one long edge the
    # plate is perfectly insulated so that the heat flux is zero through that
    # edge. The other two edges are losing heat via convection to an ambient
    # temperature of 0 °C. The thermal conductivity of the plate is 52.0 W/(m
    # .°K), and the convective heat transfer coefficient is 750 W/(m^2.°K).
    # There is no internal generation of heat. Calculate the temperature 0.2 m
    # along the un-insulated long side, measured from the intersection with the
    # fixed temperature side. The reference result is 18.25 °C.

    ##
    # The reference temperature at the point A  is 18.25 °C according to the
    # NAFEMS publication ( hich cites the book Carslaw, H.S. and J.C. Jaeger,
    # Conduction of Heat in Solids. 1959: Oxford University Press).

    ##
    # The present  tutorial will investigate the reference temperature  and it
    # will attempt to  estimate the  limit value more precisely using a
    # sequence of meshes and Richardson's extrapolation.

    ## Solution
    #

    println("""
    NAFEMS benchmark.
    Two-dimensional heat transfer with convection: convergence study.
    Solution with linear triangles.
    Version: 05/29/2017
    """
    )

    kappa = [52. 0; 0 52.]*phun("W/(M*K)"); # conductivity matrix
    h = 750*phun("W/(M^2*K)");# surface heat transfer coefficient
    Width = 0.6*phun("M");# Geometrical dimensions
    Height = 1.0*phun("M");
    HeightA = 0.2*phun("M");
    Thickness = 0.1*phun("M");
    tolerance  = Width/1000;

    m = MatHeatDiff(kappa)

    modeldata = nothing
    resultsTempA = FFlt[]
    for nref = 1:5
        t0 = time()

        # The mesh is created from two triangles to begin with
        fens,fes = T3blockx([0.0, Width], [0.0, HeightA])
        fens2,fes2 = T3blockx([0.0, Width], [HeightA, Height])
        fens,newfes1,fes2 = mergemeshes(fens, fes, fens2, fes2, tolerance)
        fes = cat(newfes1,fes2)
        # Refine the mesh desired number of times
        for ref = 1:nref
            fens,fes = T3refine(fens,fes);
        end
        bfes = meshboundary(fes)

        # Define boundary conditions

        ##
        # The prescribed temperature is applied along edge 1 (the bottom
        # edge in Figure 1)..

        l1 = selectnode(fens; box=[0. Width 0. 0.], inflate=tolerance)
        essential1 = FDataDict("node_list"=>l1, "temperature"=> 100.);

        ##
        # The convection boundary condition is applied along the edges
        # 2,3,4. The elements along the boundary are quadratic line
        # elements L3. The order-four Gauss quadrature is sufficiently
        # accurate.
        l2 = selectelem(fens, bfes; box=[Width Width  0.0 Height], inflate =tolerance)
        l3 = selectelem(fens, bfes; box=[0.0 Width Height Height], inflate =tolerance)
        cfemm = FEMMHeatDiffSurf(IntegData(subset(bfes,vcat(l2,l3)),
        GaussRule(1, 3), Thickness), h)
        convection1 = FDataDict("femm"=>cfemm, "ambient_temperature"=>0.);

        # The interior
        femm = FEMMHeatDiff(IntegData(fes, TriRule(3), Thickness), m)
        region1 = FDataDict("femm"=>femm)

        # Make the model data
        modeldata = FDataDict("fens"=> fens,
        "regions"=>[region1],
        "essential_bcs"=>[essential1],
        "convection_bcs"=>[convection1]);

        # Call the solver
        modeldata = FinEtools.AlgoHeatDiffModule.steadystate(modeldata)

        println("Total time elapsed = ",time() - t0,"s")

        l4 = selectnode(fens; box=[Width Width HeightA HeightA], inflate =tolerance)

        geom = modeldata["geom"]
        Temp = modeldata["temp"]

        ##
        # Collect the temperature  at the point A  [coordinates
        # (Width,HeightA)].
        println("$(Temp.values[l4][1])")
        push!(resultsTempA, Temp.values[l4][1]);

    end

    ##
    # These are the computed results for the temperature at point A:
    println("$( resultsTempA  )")

    # Postprocessing
    geom = modeldata["geom"]
    Temp = modeldata["temp"]
    regions = modeldata["regions"]
    vtkexportmesh("T4NAFEMS--T3.vtk", connasarray(regions[1]["femm"].integdata.fes),  [geom.values Temp.values/100], FinEtools.MeshExportModule.T3;
    scalars=[("Temperature", Temp.values)])
    vtkexportmesh("T4NAFEMS--T3--base.vtk", connasarray(regions[1]["femm"].integdata.fes), [geom.values 0.0*Temp.values/100], FinEtools.MeshExportModule.T3)

    # ##
    # # Richardson extrapolation is used to estimate the true solution from the
    # # results for the finest three meshes.
    #    [xestim, beta] = richextrapol(results(end-2:end),mesh_sizes(end-2:end));
    #     disp(['Estimated true solution for temperature at A: ' num2str(xestim) ' degrees'])

    # ##
    # # Plot the estimated true error.
    #    figure
    #     loglog(mesh_sizes,abs(results-xestim)/xestim,'bo-','linewidth',3)
    #     grid on
    #      xlabel('log(mesh size)')
    #     ylabel('log(|estimated temperature error|)')
    #     set_graphics_defaults

    # ##
    # # The estimated true error has  a slope of approximately 4 on the log-log
    # scale.
    # ##
    # # Plot the absolute values of the approximate error (differences  of
    # # successive solutions).
    #     figure
    #     loglog(mesh_sizes(2:end),abs(diff(results)),'bo-','linewidth',3)
    #     Thanksgrid on
    #     xlabel('log(mesh size)')
    #     ylabel('log(|approximate temperature error|)')
    #     set_graphics_defaults

    ## Discussion
    #
    ##
    # The last segment  of the approximate error curve is close to the slope of
    # the estimated true error. Nevertheless, it would have been more
    # reassuring if the  three successive approximate errors  were located more
    # closely on a straight line.

    ##
    # The use of uniform mesh-size meshes is sub optimal: it would be more
    # efficient to use graded meshes. The tutorial pub_T4NAFEMS_conv_graded
    # addresses use of graded meshes  in convergence studies.
end # T4NAFEMS_T3_algo

function T4NAFEMS_T6_algo()
    ## Two-dimensional heat transfer with convection: convergence study
    #

    ## Description
    #
    # Consider a plate of uniform thickness, measuring 0.6 m by 1.0 m. On one
    # short edge the temperature is fixed at 100 °C, and on one long edge the
    # plate is perfectly insulated so that the heat flux is zero through that
    # edge. The other two edges are losing heat via convection to an ambient
    # temperature of 0 °C. The thermal conductivity of the plate is 52.0 W/(m
    # .°K), and the convective heat transfer coefficient is 750 W/(m^2.°K).
    # There is no internal generation of heat. Calculate the temperature 0.2 m
    # along the un-insulated long side, measured from the intersection with the
    # fixed temperature side. The reference result is 18.25 °C.

    ##
    # The reference temperature at the point A  is 18.25 °C according to the
    # NAFEMS publication ( hich cites the book Carslaw, H.S. and J.C. Jaeger,
    # Conduction of Heat in Solids. 1959: Oxford University Press).

    ##
    # The present  tutorial will investigate the reference temperature  and it
    # will attempt to  estimate the  limit value more precisely using a
    # sequence of meshes and Richardson's extrapolation.

    ## Solution
    #

    println("""
    NAFEMS benchmark.
    Two-dimensional heat transfer with convection: convergence study.
    Solution with quadratic triangles.
    Version: 05/29/2017
    """
    )

    kappa = [52. 0; 0 52.]*phun("W/(M*K)"); # conductivity matrix
    h = 750*phun("W/(M^2*K)");# surface heat transfer coefficient
    Width = 0.6*phun("M");# Geometrical dimensions
    Height = 1.0*phun("M");
    HeightA = 0.2*phun("M");
    Thickness = 0.1*phun("M");
    tolerance  = Width/1000;

    m = MatHeatDiff(kappa)

    modeldata = nothing
    resultsTempA = FFlt[]; params = FFlt[];
    for nref = 1:5
        t0 = time()

        # The mesh is created from two triangles to begin with
        fens,fes = T3blockx([0.0, Width], [0.0, HeightA])
        fens2,fes2 = T3blockx([0.0, Width], [HeightA, Height])
        fens,newfes1,fes2 = mergemeshes(fens, fes, fens2, fes2, tolerance)
        fes = cat(newfes1,fes2)
        # Refine the mesh desired number of times
        for ref = 1:nref
            fens,fes = T3refine(fens,fes);
        end
        fens, fes = T3toT6(fens,fes);
        bfes = meshboundary(fes)

        # Define boundary conditions

        ##
        # The prescribed temperature is applied along edge 1 (the bottom
        # edge in Figure 1)..

        l1 = selectnode(fens; box=[0. Width 0. 0.], inflate=tolerance)
        essential1 = FDataDict("node_list"=>l1, "temperature"=> 100.);

        ##
        # The convection boundary condition is applied along the edges
        # 2,3,4. The elements along the boundary are quadratic line
        # elements L3. The order-four Gauss quadrature is sufficiently
        # accurate.
        l2 = selectelem(fens, bfes; box=[Width Width  0.0 Height], inflate =tolerance)
        l3 = selectelem(fens, bfes; box=[0.0 Width Height Height], inflate =tolerance)
        cfemm = FEMMHeatDiffSurf(IntegData(subset(bfes,vcat(l2,l3)),
        GaussRule(1, 3), Thickness), h)
        convection1 = FDataDict("femm"=>cfemm, "ambient_temperature"=>0.);

        # The interior
        femm = FEMMHeatDiff(IntegData(fes, TriRule(3), Thickness), m)
        region1 = FDataDict("femm"=>femm)

        # Make the model data
        modeldata = FDataDict("fens"=> fens,
        "regions"=>[region1],
        "essential_bcs"=>[essential1],
        "convection_bcs"=>[convection1]);

        # Call the solver
        modeldata = FinEtools.AlgoHeatDiffModule.steadystate(modeldata)

        println("Total time elapsed = ",time() - t0,"s")

        l4 = selectnode(fens; box=[Width Width HeightA HeightA], inflate =tolerance)

        geom = modeldata["geom"]
        Temp = modeldata["temp"]

        ##
        # Collect the temperature  at the point A  [coordinates
        # (Width,HeightA)].
        push!(resultsTempA, Temp.values[l4][1]);
        push!(params, 1.0/2^nref);

    end

    ##
    # These are the computed results for the temperature at point A:
    println("$( resultsTempA  )")

    solnestim, beta, c, residual = FinEtools.AlgoBaseModule.richextrapol(resultsTempA, params)
    println("Solution estimate = $(solnestim)")
    println("Convergence rate estimate  = $(beta )")

    # Postprocessing
    geom = modeldata["geom"]
    Temp = modeldata["temp"]
    regions = modeldata["regions"]
    vtkexportmesh("T4NAFEMS--T6.vtk", connasarray(regions[1]["femm"].integdata.fes), [geom.values Temp.values/100], FinEtools.MeshExportModule.T6;
    scalars=[("Temperature", Temp.values)])
    vtkexportmesh("T4NAFEMS--T6--base.vtk", connasarray(regions[1]["femm"].integdata.fes), [geom.values 0.0*Temp.values/100], FinEtools.MeshExportModule.T6)

    # ##
    # # Richardson extrapolation is used to estimate the true solution from the
    # # results for the finest three meshes.
    #    [xestim, beta] = richextrapol(results(end-2:end),mesh_sizes(end-2:end));
    #     disp(['Estimated true solution for temperature at A: ' num2str(xestim) ' degrees'])

    # ##
    # # Plot the estimated true error.
    #    figure
    #     loglog(mesh_sizes,abs(results-xestim)/xestim,'bo-','linewidth',3)
    #     grid on
    #      xlabel('log(mesh size)')
    #     ylabel('log(|estimated temperature error|)')
    #     set_graphics_defaults

    # ##
    # # The estimated true error has  a slope of approximately 4 on the log-log
    # scale.
    # ##
    # # Plot the absolute values of the approximate error (differences  of
    # # successive solutions).
    #     figure
    #     loglog(mesh_sizes(2:end),abs(diff(results)),'bo-','linewidth',3)
    #     Thanksgrid on
    #     xlabel('log(mesh size)')
    #     ylabel('log(|approximate temperature error|)')
    #     set_graphics_defaults

    ## Discussion
    #
    ##
    # The last segment  of the approximate error curve is close to the slope of
    # the estimated true error. Nevertheless, it would have been more
    # reassuring if the  three successive approximate errors  were located more
    # closely on a straight line.

    ##
    # The use of uniform mesh-size meshes is sub optimal: it would be more
    # efficient to use graded meshes. The tutorial pub_T4NAFEMS_conv_graded
    # addresses use of graded meshes  in convergence studies.
end # T4NAFEMS_T6_algo

function allrun()
    println("#####################################################")
    println("# T4NAFEMS_T3_algo ")
    T4NAFEMS_T3_algo()
    println("#####################################################")
    println("# T4NAFEMS_T6_algo ")
    T4NAFEMS_T6_algo()
end # function allrun

end # module T4NAFEMS_examples


text:
using JSON
using Genie, HTTP
using Genie.Router, Genie.Responses

const PORT = 8000
const HOST = "127.0.0.1"
const PROTOCAL = "http"
const BASE_URL = "$PROTOCAL://$HOST:$PORT"
const ASSETS_PATH = joinpath(dirname(dirname(@__FILE__)), "test", "assets")
const DEFAULT_SWAGGER_JSON = "https://petstore.swagger.io/v2/swagger.json"

serve_assets(ASSETS_PATH)
Genie.AppServer.startup(PORT, HOST; open_browser = false, verbose = true)

# basic settings
swagger_document = JSON.parsefile(joinpath(ASSETS_PATH, "swagger.json"))
swagger_options = Dict{String, Any}()
options = Options()

swagger_options["url"] = "https://petstore.swagger.io/v2/swagger.json"
options.swagger_options = swagger_options

@testset "Render Assets Tests" begin
    for (root, dirs, files) in walkdir(ASSETS_PATH)
        for file in files
            r = HTTP.request("GET", "$BASE_URL/$file")
            local_file = read(open(joinpath(ASSETS_PATH, file)), String)
            @test r.status == 200
            @test local_file == String(r.body)
        end
    end
end

@testset "Basic Swagger Tests" begin
    route_name = "docs"
    route("/$route_name", method = GET) do
        render_swagger(nothing, options=options)
    end

    r = HTTP.request("GET", "$BASE_URL/$route_name")
    html_string = String(r.body)
    opts = Dict{String, Any}()
    opts["customOptions"] = swagger_options

    @test r.status == 200
    @test occursin("<title>$(options.custom_site_title)</title>", html_string)
    @test occursin("const options = $(JSON.json(opts))", html_string)
end

@testset "Customizations Tests" begin
    route_name = "docs_custom"
    custom_options = Options()
    custom_swagger_options = Dict{String, Any}()
    custom_swagger_options["url"] = DEFAULT_SWAGGER_JSON
    custom_options.swagger_options = custom_swagger_options

    # custom css
    css_string = ".swagger-ui .topbar { display: none }"
    custom_options.custom_css = css_string

    # custom stylesheet
    stylesheet_path = joinpath(ASSETS_PATH, "custom.css")
    custom_options.custom_css_url = stylesheet_path

    # custom favicon
    favicon_path = joinpath(ASSETS_PATH, "favicon.ico")
    custom_options.custom_favicon = favicon_path

    # custom site title
    title = "Jimmy's new API"
    custom_options.custom_site_title = title

    route("/$route_name", method = GET) do
        render_swagger(nothing, options=custom_options)
    end

    r = HTTP.request("GET", "$BASE_URL/$route_name")
    html_string = replace(String(r.body), "\n" => "")

    @test r.status == 200
    @test occursin(css_string, html_string)
    @test occursin("<link rel='stylesheet' type='text/css' href=$(stylesheet_path) />", html_string)
    @test occursin("<title>$(title)</title>", html_string)
    @test occursin("<link rel='icon' href=$(favicon_path) />", html_string)
end


text:
module TestComposite

using Test
using Mimi

import Mimi:
    ComponentId, ComponentPath, DatumReference, ComponentDef, AbstractComponentDef, CompositeComponentDef,
    Binding, ModelDef, build, time_labels, compdef, find_comp

@defcomp Comp1 begin
    par_1_1 = Parameter(index=[time])      # external input
    var_1_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_1_1[t] = p.par_1_1[t]
    end
end

@defcomp Comp2 begin
    par_2_1 = Parameter(index=[time])      # connected to Comp1.var_1_1
    par_2_2 = Parameter(index=[time])      # external input
    var_2_1 = Variable(index=[time])       # computed
    foo = Parameter()

    function run_timestep(p, v, d, t)
        v.var_2_1[t] = p.par_2_1[t] + p.foo * p.par_2_2[t]
    end
end

@defcomp Comp3 begin
    par_3_1 = Parameter(index=[time])      # connected to Comp2.var_2_1
    var_3_1 = Variable(index=[time])       # external output
    foo = Parameter(default=30)

    function run_timestep(p, v, d, t)
        # @info "Comp3 run_timestep"
        v.var_3_1[t] = p.par_3_1[t] * 2
    end
end

@defcomp Comp4 begin
    par_4_1 = Parameter(index=[time])      # connected to Comp2.var_2_1
    var_4_1 = Variable(index=[time])        # external output
    foo = Parameter(default=300)

    function run_timestep(p, v, d, t)
        # @info "Comp4 run_timestep"
        v.var_4_1[t] = p.par_4_1[t] * 2
    end
end

m = Model()
set_dimension!(m, :time, 2005:2020)

@defcomposite A begin
    component(Comp1)
    component(Comp2)

    foo1 = Comp1.foo
    foo2 = Comp2.foo

    # Should accomplish the same as calling
    #   connect_param!(m, "/top/A/Comp2:par_2_1", "/top/A/Comp1:var_1_1")
    # after the `@defcomposite top ...`
    Comp2.par_2_1 = Comp1.var_1_1
    Comp2.par_2_2 = Comp1.var_1_1
end

@defcomposite B begin
    component(Comp3) # bindings=[foo => bar, baz => [1 2 3; 4 5 6]])
    component(Comp4)

    foo3 = Comp3.foo
    foo4 = Comp4.foo
end

@defcomposite top begin
    component(A)

    fooA1 = A.foo1
    fooA2 = A.foo2

    # TBD: component B isn't getting added to mi
    component(B)
    foo3 = B.foo3
    foo4 = B.foo4
end

# We have created the following composite structure:
#
#          top
#        /    \
#       A       B
#     /  \     /  \
#    1    2   3    4

top_ref = add_comp!(m, top, nameof(top))
top_comp = compdef(top_ref)

md = m.md

@test find_comp(md, :top) == top_comp

#
# Test various ways to access sub-components
#
c1 = find_comp(md, ComponentPath(:top, :A, :Comp1))
@test c1.comp_id == Comp1.comp_id

c2 = md[:top][:A][:Comp2]
@test c2.comp_id == Comp2.comp_id

c3 = find_comp(md, "/top/B/Comp3")
@test c3.comp_id == Comp3.comp_id

set_param!(m, "/top/A/Comp1:foo", 1)
set_param!(m, "/top/A/Comp2:foo", 2)

# TBD: default values set in @defcomp are not working...
# Also, external_parameters are stored in the parent, so both of the
# following set parameter :foo in "/top/B", with 2nd overwriting 1st.
set_param!(m, "/top/B/Comp3:foo", 10)
set_param!(m, "/top/B/Comp4:foo", 20)

set_param!(m, "/top/A/Comp1", :par_1_1, collect(1:length(time_labels(md))))

# connect_param!(m, "/top/A/Comp2:par_2_1", "/top/A/Comp1:var_1_1")
# connect_param!(m, "/top/A/Comp2:par_2_2", "/top/A/Comp1:var_1_1")
connect_param!(m, "/top/B/Comp3:par_3_1", "/top/A/Comp2:var_2_1")
connect_param!(m, "/top/B/Comp4:par_4_1", "/top/B/Comp3:var_3_1")

build(m)

run(m)

#
# TBD
#
# 1. Create parallel structure of exported vars/pars in Instance hierarchy?
#    - Perhaps just a dict mapping local name to a component path under mi, to where the var actually exists
# 2. Be able to connect to the leaf version of vars/pars or by specifying exported version below the compdef
#    given as first arg to connect_param!().
# 3. set_param!() should work with relative path from any compdef.
# 4. set_param!() stores external_parameters in the parent object, creating namespace conflicts between comps.
#    Either store these in the leaf or store them with a key (comp_name, param_name)

mi = m.mi

@test mi[:top][:A][:Comp2, :par_2_2] == collect(1.0:16.0)
@test mi["/top/A/Comp2", :par_2_2] == collect(1.0:16.0)

@test mi["/top/A/Comp2", :var_2_1] == collect(3.0:3:48.0)
@test mi["/top/A/Comp1", :var_1_1] == collect(1.0:16.0)
@test mi["/top/B/Comp4", :par_4_1] == collect(6.0:6:96.0)

end # module

m = TestComposite.m
md = m.md
top   = Mimi.find_comp(md, :top)
A     = Mimi.find_comp(top, :A)
comp1 = Mimi.find_comp(A, :Comp1)

nothing


text:
### visualize a solution.

function visualizefiltersolution(𝓧, xq, wfunc, dfunc,f)
    𝑤_new = get𝑤scaled(𝓧)
    h_new = geth(𝑤_new, wfunc, dfunc, 𝓧)

    f_itp = xx->getp(xx, 𝑤_new, h_new, wfunc, dfunc, 𝓧)
    f_itp_xq = collect( f_itp(xq[i]) for i = 1:length(xq) )
    f_itp_𝓧 = collect( f_itp(𝓧[i]) for i = 1:length(𝓧) )

    𝑒 = xx->(f(xx)-f_itp(xx))
    𝑒_xq = collect( 𝑒(xq[i]) for i = 1:length(xq) )
    𝑒_𝓧 = collect( 𝑒(𝓧[i]) for i = 1:length(𝓧) )

    return f_itp_xq, f_itp_𝓧, 𝑒_xq, 𝑒_𝓧
end

function computeDTFTviaformula(h::AbstractArray{T}, ω::T)::Complex{T} where T <: Real
    N = length(h)
    return sum( h[n+1]*exp(-im*ω*n) for n = 0:N-1 )
end

function computeDTFTviaformula(h::AbstractArray{T}, ω::Vector{T}) where T

    return collect( computeDTFTviaformula(h, ω[i]) for i = 1:length(ω))
end

function getfreqrsp(h::Vector{Float64}, resolution_multiple::Int = 20)
    N_samples = length(h)

    ω_set_fft = collect( LinRange(0,2*π-2*π/N_samples,N_samples))
    DFT_evals = fft(h)

    ω_set = collect( LinRange(0,2*π,N_samples*resolution_multiple) )
    DTFT_evals = collect( computeDTFTviaformula(h, ω_set[i]) for i = 1:length(ω_set) )

    return ω_set_fft, DFT_evals, ω_set, DTFT_evals
end

function plotmagnitudersp(h::Vector{Float64}, fig_num::Int, title_string::String = "Magnitude response")

    ω_set_fft, DFT_evals, ω_set, DTFT_evals = getfreqrsp(h)

    # visualize.
    mag_rsp_G = abs.(DTFT_evals)
    mag_rsp_fft = abs.(DFT_evals)

    PyPlot.figure(fig_num)
    fig_num += 1

    PyPlot.plot(ω_set_fft, mag_rsp_fft, ".", label = "DFT")
    PyPlot.plot(ω_set, mag_rsp_G, label = "DTFT")

    PyPlot.title(title_string)
    PyPlot.legend()

    return fig_num
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
testdir = joinpath(dirname(pathof(CombinedParsers)),"../test/")
docdir = joinpath(dirname(pathof(CombinedParsers)),"../docs/src/")
include(joinpath(testdir, "pcretest-parser.jl"))

tests_string=read(joinpath(testdir,"testoutput1"),String);
tests = pcre_tests(tests_string)[1];

tt=tests[14]
@pcre_testset tt[2] true

ignore_idx = Dict{Int,String}(
    ## 14 => "unicode escape in test parser needs to support \"\\u81\".",
    69 => "because of very long compile time. The complex pattern parses email adresses",
    70 => "because of very long compile time. The complex pattern parses email adresses",
    89 => "CombinedParsers brackets use unicode character ranges. Avoid such patterns for now, contributions very welcome.",
    90 => "CombinedParsers brackets use unicode character ranges. Avoid such patterns for now, contributions very welcome.",
    135 => "escaped characters needs investigation.",
    664 => "because of very long matching time for `Repeat(Repeat1('a'))` when there is no match. Avoid such patterns, contributions optimizing these cases are also very welcome.", ## slooow
    887 => raw"""do I misunderstand [recursive backreferences](https://www.pcre.org/original/doc/html/pcrepattern.html#SEC19)?

```@repl session
pattern = "(?P<abn>(?P=abn)xxx|)+"
match(Regex(pattern),"xxx")
match(Regcomb(pattern),"xxx")
```

IMO the pcre behaviour is confusing after considering the logic of expanding the repeat
```@repl session
pattern = "(?P<abn>(?P=abn)xxx|)((?P=abn)xxx|)"
match(Regex(pattern),"xxx")
match(Regcomb(pattern),"xxx")
```

```@repl session
parse(re"(\1xxx|)+","xxx")
match(r"(\1xxx|)+","xxx")
parse(re"(|\1xxx)+", "xxx")
parse(r"(|\1xxx)+", "xxx")
```

""",
    ## 104 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    ## 260 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    ## 261 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    ## 562 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    524 => "Captures in Lookaheads need refactoring. Low prio because instead of using backreference, `FlatMap` is recommended in `CombinedParsers.jl`.",
    527 => "Capture backtracking needs refactoring. Low prio because instead of using captures, `parse` is recommended in `CombinedParsers.jl`.",
    528 => "Capture backtracking needs refactoring. Low prio because instead of using captures, `parse` is recommended in `CombinedParsers.jl`.",
    625 => "\n!!! note\n    Lookahead in Lookbehind needs refactoring.",
    ## 1146 => "case-ignoring backreferences still need to be done, contributions are very welcome,",
    1173 => "because of very long compile time. The complex pattern parses PCRE regex syntax."
)

optimize_idx = []

import Test: record, finish
using Test: AbstractTestSet, Result, Pass, Fail, Error
using Test: get_testset_depth, get_testset
struct CustomTestSet <: Test.AbstractTestSet
    description::AbstractString
    results::Vector
    # constructor takes a description string and options keyword arguments
    CustomTestSet(desc) = new(desc, [])
end

record(ts::CustomTestSet, child::AbstractTestSet) = push!(ts.results, child)
record(ts::CustomTestSet, res::Result) = push!(ts.results, res)
function finish(ts::CustomTestSet)
    # just record if we're not the top-level parent
    if get_testset_depth() > 0
        record(get_testset(), ts)
    end
    ts
end

unsupported=Dict{Any,Any}()
errors=Any[]
testresults = @testset CustomTestSet "pcre testset 1" begin
    for (i,tt) in enumerate(tests)
        tr = @testset "$i" begin
            if !haskey(ignore_idx,i) && !in(i,optimize_idx)
                nam = string(tt[2].pattern)[1:min(end,20)]
                try
                    println(i)
                    @pcre_testset tt[2] false
                catch e
                    ##sleep(.1)
                    if e isa UnsupportedError
                        @warn "unsupported" e.message tt[2].pattern
                        ##readline()
                        global unsupported
                        push!(get!(()->Any[],unsupported,e.message),(i,tt[2]))
                        (results =[],)
                    else
                        @warn "error in $i" tt[2].pattern exception=e
                        ##readline()
                        global errors
                        push!(errors,(i,tt[2]))
                    end
                end
            end
        end
        nfailed = sum([ !(t isa Pass) for t in tr.results ])
        if nfailed>0
            @warn "error in $i $(tt[2].pattern)"
            global err_idx,err_test=i,tt[2]
            global errors
            push!(errors,(i,tt[2]))
            ##readline()
        end
    end
end

unsupportednos = Set(
    vcat([ first.(v) for v in values(unsupported) ]...)
)
successes = Dict()
failures = Dict()
for ts in testresults.results
    i = parse(Int,ts.description)
    nsuccess = sum([ t isa Pass for t in ts.results ])
    nfailed = length(ts.results) - nsuccess
    if nfailed == 0
        if nsuccess > 0
            successes[i] = (nsuccess, length(ts.results))
        end
    elseif !(i in unsupportednos)
        failures[i] = (nfailed,length(ts.results))
    end
end

n_successes, n_failed = sum(getindex.(values(successes),1)), sum(getindex.(values(failures),1))
n_patterns = (
    success = length(successes),
    failed = length(failures),
    skipped = length(ignore_idx),
    unsupported = sum(length.(values(unsupported))))
io=stdout

benchmarks = [
    # (range_Regcomb = (33.0, 28515.5), range_Regex = (83.0, 682.5), proportion_better = 0.5093632958801498, mean_ratio = 2.7263015758568434, ratio_mean = 1.2951982070322912)
    "benchmark-2020-06-29_10h29-4304c235feb678ca0ea7902a1a296dd2dace0ae4.json" =>
    (range_Regcomb = (34.0, 19704.0), range_Regex = (86.0, 521.0), proportion_better = 0.5917602996254682, mean_ratio = 1.2631184139535836, ratio_mean = 0.8598574284168594)
][end].second

function print_testset(io,prefix,t,testdef; pad="   ")
    pat,opt = t.pattern
    pat = replace(pad*"/"*pat, "\n" => "\n"*pad)*"/"*opt
    println(io,"\n$pad```\n$pat\n$pad```\n$pad$prefix")
    testseqs = vcat([ "`"*s.sequence*"`" for s in t.test],"`".*t.tests_nomatch .* "` no match")
    for (j,tr) in enumerate(testdef.results)
        if (j>1)
            s = Base.escape_string(testseqs[j-1])
            result = tr isa Pass ? "✅" : "❌"
            println(io,pad*"- $result $s")
        end
    end
end

open(joinpath(docdir,"man","pcre-compliance.md"),"w") do io
  println(io,"""
# Compliance with the PCRE test set
!!! note
    PCRE features supported by `@re_str`
    - ✅ sequences, alternations (`|`), repetitions (`*`,`+`,`{n}`, `{min,}`, `{min,max}`), optional matches (`?`)
    - ✅ escaped characters and generic character types
    - ✅ character ranges (`[]`)
    - ✅ non-capturing groups
    - ✅ capturing groups, backreferences, subroutines (all by index, relative index and name)
    - ✅ simple assertions (`\\A`, `\\z`, `\\Z`, `\\b`, `\\B`, `^`, `\$`)
    - ✅ lookaheads and lookbehinds
    - ✅ atomic groups
    - ✅ lazy repetitions
    - ✅ conditional expressions
    - ✅ internal and pattern options setting
    - ✅ comments
!!! warning
    PCRE functionality that is currently not supported:
    - ❌ Capture groups in lookbehinds.
    - ❌ Lookaheads within lookbehinds.
    - ❌ ACCEPT, SKIP, COMMIT, THEN, PRUNE, \\K
```@setup session
using CombinedParsers
using CombinedParsers.Regexp
```
CombinedParsers.jl is tested and benchmarked against the PCRE C library testset.
The PCRE test output is downloaded from
[the PCRE source repository](https://github.com/rurban/pcre/blob/master/testdata/testoutput1),
parsed with
[a `CombinedParser`](https://github.com/gkappler/CombinedParsers.jl/blob/master/test/pcretest-parser.jl), to run tests benchmarks on `Base.Regex` and `CombinedParsers.Regexp.Regcomb`.
(Note: tests are relaxed for some cases allowing empty captures (`""`) for unset captures (`nothing`).
## Test Overview
$n_successes successful tests on $(n_patterns.success) patterns
(See [list of compliant patterns](pcre-compliance-succeeded.md)).\n
$n_failed failed tests on $(n_patterns.failed) patterns
(See [list of failed patterns](pcre-compliance-failed.md)).
### Performance Overview:
The PCRE C backend of `@r_str` has arrived at a widely optimized codebase after decades of improvements.
C PCRE2 optimized is among the fastest regex libraries ([second behind Rust](https://github.com/mariomka/regex-benchmark/tree/optimized), running [mariomka](https://github.com/mariomka)'s benchmark will position CombinedParser among its competition).

Although CombinedParsers.jl is a very young package that will be optimized further,
`@re_str` pure Julia Regcomb is often competitive with PCRE `@r_str` Regex.\n\n
PCRE benchmarks have a range between $(benchmarks.range_Regex[1])ns to $(benchmarks.range_Regex[2])ns.
CombinedParsers benchmarks range between $(benchmarks.range_Regcomb[1])ns to $(benchmarks.range_Regcomb[2])ns.
$(round(benchmarks.proportion_better*100))% of benchmarks are faster with CombinedParsers compared to PCRE.
The average ratio of `time_Recomb/time_Regex` is $(round(benchmarks.mean_ratio,digits=2)).\n\n
These benchmarkin results are for the first 100 test patterns in the PCRE test set, comparing `match(Regex(pattern,flags),s)` with `_iterate(Regcomb(pattern,flags),s)`.
![](log_btimes.png)\n\n
Benchmark timings for regular expression construction and matching comparing `Regex` (x axis) and `Regcomb` (y axis), both on a log10 scale.
Points represent an individual benchmark for a pattern construction or match.
Cases with `CombinedParsers` being faster than the C library PCRE are paint green, slower cases are red.
### Benchmark ratios histogram:
![](log_btime_ratio_histogram.svg)\n\n

The histograms of ratios of `time_Regcomb/time_Regex` on a log scale demonstrate that `CombinedParser` implementation is competitive.
Worst cases are investigated for further optimization [in this IJulia notebook](https://github.com/gkappler/CombinedParsers.jl/blob/master/benchmark/benchmarks.ipynb).

Next steps in optimization are
- caching codeunit lengths of matches for backtracking.
- memoization of sub-parsings.
## Unsupported
$(n_patterns.unsupported) unsupported patterns were omitted for the following reasons:
""")
    for (n,r) in sort([ s.first => length(s.second) for s in unsupported ])
        println(io,"- `$n` excluded $r patterns.")
    end
end

open(joinpath(docdir,"man","pcre-compliance-failed.md"),"w") do io
    println(io,"""\n
    # Failed PCRE Tests
    $n_failed failed tests on $(n_patterns.failed) patterns.
    """)
    println(io,"""\n
## Skipped
$(n_patterns.skipped) patterns were skipped for the following reasons:
""")
    for (i,r) in sort(collect(ignore_idx))
        t = tests[i][2]
        pat,opt = t.pattern
        @show lastindex(pat)
        println(io,"\n---")
        if lastindex(pat)>100
            print(io,"\n\n(no $i) skipped, $r\n")
        else
            print_testset(io,"\n(no $i) skipped, $r\n",t,testresults.results[i], pad="")
        end
    end
    println(io,"""\n
## Failed tests
""")
    for (i,(r,n)) in sort(collect(failures))
        if n>0
            t = tests[i][2]
            println(io,"\n---")
            print_testset(io,r==0 ? "\n(no $i) failed to compile" : "(no $i) failed  $(r-1) of $(n-1) times:",t,testresults.results[i], pad="")
        end
    end
end

open(joinpath(docdir,"man","pcre-compliance-succeeded.md"),"w") do io
    println(io,"""
# PCRE Compliance
$n_successes successful tests on $(n_patterns.success) patterns.\n
""")
    for (i,(r,n)) in sort(collect(successes))
        t = tests[i][2]
        println(io,"\n---")
        print_testset(io,"(no $i) succeeded $(r-1) of $(n-1) times:\n",t,testresults.results[i],pad="")
    end
end


text:
## run.jl : entry point method for simulating an Estimator
#
# Runs the Estimator for the given tolerance, and returns a History object.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

"""
run(estimator::Estimator, ε::Real)
run(estimator::Estimator, ε::Vector{<:Real})

Run the estimator and compute the expected value of the quantity of interest up to the given tolerance(s) `ε`. Returns a [`History`](@ref)-object that contains usefull diagnostics about the simulation.

# Examples

An example detailing how to use MultilevelEstimators is provided in the [example in the documentation](@ref Example).

See also: [`Estimator`](@ref), [`History`](@ref)
"""
function run(estimator::Estimator, tols::AbstractVector{<:Real})

    # input checking
    check_larger_than("run", tols, "supplied tolerance(s)", 0)

    # make history
    history = History()

    # run the sequence of tolerances
    for tol in tols
        clear(estimator) # prepare new run
        t = @elapsed _run(estimator,tol)
        update_history!(history, estimator, tol, t) # log the results in history
    end

    return history
end

run(estimator::Estimator, tol::Real) = run(estimator, get_tols(estimator, tol))

## Main routine ##
function _run(estimator::Estimator{T1, T2}, ϵ::Real) where {T1<:AbstractIndexSet, T2<:AbstractSampleMethod}

    # print status
    estimator[:verbose] && print_header(estimator, ϵ)

    # start "level" is 0
    L = 0

    # initial MSE splitting parameter
    θ = T1 <: SL ? 0.5 : estimator[:min_splitting]

    # main loop
    is_converged = false
    while !is_converged

        # print level
        estimator[:verbose] && print_level(estimator, L)

        # update index set
        index_set = new_index_set(estimator, L)
        estimator[:verbose] && print_index_set(estimator, index_set)

        # obtain initial variance estimate
        ns = regress_nb_of_samples(estimator, index_set, ϵ, θ, L)
        for index in index_set
            if !has_samples_at_index(estimator, index)
                sample!(estimator, index, ns[index])
            end
        end

        # add new indices to the index set
        for index in index_set
            push!(estimator, index)
        end
        set_sz(estimator, L)

        # print status
        estimator[:verbose] && print_status(estimator)

        # update value of the MSE splitting parameter
        θ = T1 <: SL ? 0.5 : estimator[:do_mse_splitting] ? compute_splitting(estimator, ϵ) : estimator[:min_splitting]

        if T2 <: MC

            # evaluate optimal number of samples
            n_opt = Dict(τ => optimal_nb_of_samples(estimator, τ, ϵ, θ) for τ in keys(estimator))

            # print optimal number of samples
            estimator[:verbose] && print_optimal_nb_of_samples(estimator, n_opt)

            # take additional samples if required
            update_samples(estimator, n_opt)

        else

            while varest(estimator) > θ*ϵ^2
                max_index = find_index_with_max_var_over_cost(estimator)

                # increase the number of samples already taken
                n_opt = next_number_of_samples(estimator, max_index)

                # print optimal number of samples
                estimator[:verbose] && print_optimal_nb_of_samples(estimator, n_opt)

                # take additional samples
                update_samples(estimator, n_opt)

                # recompute splitting parameter
                θ = T1 <: SL ? 0.5 : estimator[:do_mse_splitting] ? compute_splitting(estimator, ϵ) : estimator[:min_splitting]

                # check next iteration
                estimator[:verbose] && print_qmc_convergence(estimator, ϵ, θ)
            end

        end

        # check for convergence
        if T1 <: SL
            is_converged = true
        else
            # show status
            estimator[:verbose] && print_status(estimator)
            estimator[:verbose] && L ≥ 2 && print_mse_analysis(estimator, ϵ, θ)

            # check if the new level exceeds the maximum level
			is_converged = L > estimator[:min_index_set_param] && converged(estimator, ϵ, θ)
            if !is_converged && max_level_exceeded(estimator)
                estimator[:verbose] && warn_max_level(estimator)
                is_converged = true
            else
                # update level
                L += 1
            end
        end
    end

    # update boundary in case of AD
    T1 <: AD && update_boundary(estimator)

    # print convergence status
    estimator[:verbose] && print_convergence(estimator, T1 <: SL ? true : converged(estimator, ϵ, θ))
end

## Unbiased estimator routine ##
function _run(estimator::Estimator{<:U, <:AbstractSampleMethod}, ϵ::Real)

    # print status
    estimator[:verbose] && print_header(estimator, ϵ)

    # update index set
    index_set = collect(keys(pmf(estimator)))
    estimator[:verbose] && print_index_set(estimator, index_set)

    # add new indices to the index set
    for index in index_set
        index ∈ current_index_set(estimator) || push!(estimator, index)
    end

    # main loop
    is_converged = varest(estimator) ≤ ϵ^2
    while !is_converged

        # print status
        estimator[:verbose] && print_status(estimator)

        # increase the number of samples already taken
        n_opt = next_number_of_samples(estimator)

        # print optimal number of samples
        estimator[:verbose] && print_optimal_nb_of_samples(estimator, n_opt)

        # take additional samples
        update_samples(estimator, n_opt)

        # update pmf
        update_pmf(estimator)
        estimator[:verbose] && print_pmf(estimator)

        # check next iteration
        estimator[:verbose] && print_unbiased_convergence(estimator, ϵ)

        is_converged = varest(estimator) ≤ ϵ^2

        # print index set
        estimator[:verbose] && print_index_set(estimator, index_set)
    end

    # print convergence status
    estimator[:verbose] && print_convergence(estimator, true)
end