text:
#= basic support for sparse matrices, more Magma style than Julia
  a sparse matrix is a collection (Array) of sparse rows (SmatRow)
  A sparse matrix cannot have 0 rows!
  A SmatRow is two Arrays, one (pos) containing the columns, values
  contains the values.
  A[i,j] = A.rows[i].values[A.rows[i].pos[j]]
  to be formal

  The upper_triangular stuff performs very elementary transformations
  until the matrix becomes dense. At this point, the dense bit is extraced and
  converted to an fmpz_mat which is then hnf'ed.

  Missing:
   full HNF, Howell, modular and not
   better elemination strategy
=#

import Base.push!, Base.max, Nemo.nbits, Base.sparse, Base.Array

export upper_triangular, vcat!, show, sub, Smat, SmatRow, random_SmatSLP,
       fmpz_mat, rows, cols, copy, push!, mul, mul!, abs_max, toNemo, sparse,
       valence_mc, swap_rows!

################################################################################
#
#  Base rings for sparse matrices
#
################################################################################

# This doesn't work for empty sparse matrices.
# The ring should actually be put into the type.
base_ring(A::Smat) = parent(A.rows[1].values[1])

################################################################################
#
#  SmatRow
#
################################################################################

function show{T}(io::IO, A::SmatRow{T})
  print(io, "sparse row with positions $(A.pos) and values $(A.values)\n")
end

################################################################################
#
#  Smat
#
################################################################################

@doc """
  Smat(A::fmpz_mat) -> Smat{fmpz}

  Constructs the Smat (Hecke-sparse matrix) with coefficients fmpz
  corresponding to A.

""" ->
function Smat(A::fmpz_mat)
  m = Smat{fmpz}()
  m.c = cols(A)
  m.r = 0
  for i=1:rows(A)
    if is_zero_row(A, i)
      continue
    end
    r = SmatRow{fmpz}()
    for j =1:cols(A)
      if A[i,j] != 0
        m.nnz += 1
        push!(r.values, fmpz(A[i,j]))
        push!(r.pos, j)
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

@doc """
  Smat{T}(A::Array{T, 2}) -> Smat{T}

  Constructs the Smat (Hecke-sparse matrix) with coefficients of
  type T corresponding to A.

""" ->
function Smat{T}(A::Array{T, 2})
  m = Smat{T}()
  m.c = Base.size(A, 2)
  m.r = 0
  for i=1:Base.size(A, 1)
    if is_zero_row(A, i)
      continue
    end
    r = SmatRow{T}()
    for j =1:Base.size(A, 2)
      if A[i,j] != 0
        m.nnz += 1
        push!(r.values, T(A[i,j]))
        push!(r.pos, j)
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

################################################################################
#
#  Sparse matrix from dense (with more options)
#
################################################################################

# The the entries of the sparse matrix will be elements of the ring R
# It defaults to the base ring of the input matrix.
# If zerorows is set, then zero rows will not be remove

function _Smat{T <: MatElem, S <: Ring}(A::T; R::S = base_ring(A),
                                              zerorows::Bool = false)

  m = Smat{elem_type(R)}()
  m.c = cols(A)
  m.r = 0

  for i=1:rows(A)
    if is_zero_row(A, i)
      if !zerorows
        continue
      else
        r = SmatRow{elem_type(R)}()
        #push!(m.rows, SmatRow{elem_type(R)}())
      end
    else
      r = SmatRow{elem_type(R)}()
      for j = 1:cols(A)
        t = A[i, j]
        if t != 0
          m.nnz += 1
          push!(r.values, R(t))
          push!(r.pos, j)
        end
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

# a faster version for nmod_mat -> Smat{T}
# it avoids the creation of elements in ResidueRing(ZZ, n)
function Smat{S <: Ring}(A::nmod_mat; R::S = base_ring(A), zerorows::Bool = false)
  if R == base_ring(A)
    return _Smat(A, R)
  end

  m = Smat{elem_type(R)}()
  m.c = cols(A)
  m.r = 0

  for i=1:rows(A)
    if is_zero_row(A, i)
      if !zerorows
        continue
      else
        #push!(m.rows, SmatRow{elem_type(R)}())
        r = SmatRow{elem_type(R)}()
      end
    else
      r = SmatRow{elem_type(R)}()
      for j =1:cols(A)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &A, i - 1, j - 1)
        if t != 0
          m.nnz += 1
          push!(r.values, R(t))
          push!(r.pos, j)
        end
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

# Smat{fmpz} -> Smat{UIntMod}
function mod(A::Smat{fmpz}, n::Int)
  z = Smat{UIntMod}()
  R = ZZModUInt(UInt(n))
  z.r = A.r
  z.c = A.c
  for r in A.rows
    rz = SmatRow{UIntMod}()
    for (i, v) in zip(r.pos, r.values)
      if mod(v, n) == 0
        continue
      else
        push!(rz.pos, i)
        push!(rz.values, R(v))
      end
    end
    if length(rz.pos) != 0
      push!(z.rows, rz)
      z.nnz += length(rz.pos)
    end
  end
  return z
end

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

function show{T}(io::IO, A::Smat{T})
  print(io, "Sparse ", A.r, " x ", A.c, " matrix with ")
  print(io, A.nnz, " non-zero entries\n")
end

################################################################################
#
#  Conversion
#
################################################################################

@doc """
  toNemo(io::IOStream, A::Smat; name = "A")

  Prints the Smat as a julia-program into the file corresponding to io.
  The file can be included to get the matrix.
  `name` controls the variable name of the matrix.
"""->
function toNemo(io::IOStream, A::Smat; name = "A")
  T = typeof(A.rows[1].values[1])
  println(io, name, " = Smat{$T}()")
  for i=A.rows
    print(io, "push!($(name), SmatRow{$T}(Tuple{Int, $T}[");
    for j=1:length(i.pos)-1
      print(io, "($(i.pos[j]), $(i.values[j])), ")
    end
    println(io, "($(i.pos[end]), $(i.values[end]))]))")
  end
end

@doc """
  toNemo(io::ASCIIString, A::Smat; name = "A")

  Prints the Smat as a julia-program into the file named io.
  The file can be included to get the matrix.
  `name` controls the variable name of the matrix.
"""->
function toNemo(f::ASCIIString, A::Smat; name = "A")
  io = open(f, "w")
  toNemo(io, A, name=name)
  close(io)
end

################################################################################
#
#  Transpose
#
################################################################################

@doc """
  transpose{T}(A::Smat{T}) -> Smat{T}

  The transpose of A
""" ->
function transpose{T}(A::Smat{T})
  B = Smat{T}()
  n = rows(A)
  m = cols(A)
  B.rows = Array(SmatRow{T}, m)
  for i=1:m
    B.rows[i] = SmatRow{T}()
  end
  for i = 1:length(A.rows)
    for j = 1:length(A.rows[i].pos)
      r = A.rows[i].pos[j]
      push!(B.rows[r].pos, i)
      push!(B.rows[r].values, A.rows[i].values[j])
    end
  end
  B.c = n
  B.r = m
  B.nnz = A.nnz

  return B
end

################################################################################
#
# multiplication of Smat * Dense
# for Dense = {Array{T, 1}, Array{T, 2}, fmpz_mat}
# both mul and mul! versions
#
################################################################################
function mul!{T}(c::Array{T, 1}, A::Smat{T}, b::Array{T, 1})
  assert( length(b) == cols(A))
  assert( length(c) == rows(A))
  for i = 1:length(A.rows)
    s = 0
    I = A.rows[i]
    for j=1:length(I.pos)
      s += I.values[j]*b[I.pos[j]]
    end
    c[i] = s
  end
  return c
end

function mul{T}(A::Smat{T}, b::Array{T, 1})
  assert( length(b) == cols(A))
  c = Array(T,rows(A))
  mul!(c, A, b)
  return c
end

function mul_mod_big!{S, T}(c::Array{S, 1}, A::Smat{T}, b::Array{S, 1}, mod::S)
  assert( length(b) == cols(A))
  assert( length(c) == rows(A))
  for i = 1:length(A.rows)
    s = 0
    I = A.rows[i]
    for j=1:length(I.pos)
      s += S(I.values[j] *b[I.pos[j]] % mod)
    end
    c[i] = s % mod
  end
  return c
end

function mul_mod!{S, T}(c::Array{S, 1}, A::Smat{T}, b::Array{S, 1}, mod::S)
  assert( length(b) == cols(A))
  assert( length(c) == rows(A))
  for i = 1:length(A.rows)
    s = 0
    I = A.rows[i]
    for j=1:length(I.pos)
      s += S(I.values[j]) *b[I.pos[j]]
    end
    c[i] = s % mod
  end
  return c
end

function mul!{T}(c::Array{T, 2}, A::Smat{T}, b::Array{T, 2})
  sz = size(b)
  assert( sz[1] == cols(A))
  tz = size(c)
  assert( tz[1] == rows(A))
  assert( tz[2] == sz[2])
  for m = 1:sz[2]
    for i = 1:length(A.rows)
      s = 0
      for j=1:length(A.rows[i].pos)
        s += A.rows[i].values[j]*b[A.rows[i].pos[j], m]
      end
      c[i, m] = s
    end
  end
  return c
end

function mul{T}(A::Smat{T}, b::Array{T, 2})
  sz = size(b)
  assert( sz[1] == cols(A))
  c = Array(T, sz[1], sz[2])
  return mul!(c, A, b)
end

function mul!{T}(c::fmpz_mat, A::Smat{T}, b::fmpz_mat)
  assert( rows(b) == cols(A))
  assert( rows(c) == rows(A))
  assert( cols(c) == cols(b))
  for m = 1:cols(b)
    for i = 1:length(A.rows)
      s = 0
      for j=1:length(A.rows[i].pos)
        s += A.rows[i].values[j]*b[A.rows[i].pos[j], m]
      end
      c[i, m] = s
    end
  end
  return c
end

function mul{T}(A::Smat{T}, b::fmpz_mat)
  assert( rows(b) == cols(A))
  c = MatrixSpace(ZZ, rows(A), cols(b))()
  return mul!(c, A, b)
end

function SLP_AddRow{T}(i::Int, j::Int, v::T)
  assert(v != 0)
  slp = SmatSLP_add_row(i, j, v)
  return slp
end

function SLP_SwapRows(i::Int, j::Int)
  slp = SmatSLP_swap_row(i, j)
  return slp
end

function mul!{T}(a::Array{T, 1}, s::SmatSLP_swap_row)
  t = a[s.row]
  a[s.row] = a[s.col]
  a[s.col] = t
end

function mul!{T}(a::Array{T, 1}, s::SmatSLP_add_row{T})
  a[s.row] = a[s.row]*s.val + a[s.col]
end

function mul_t!{T}(a::Array{T, 1}, s::SmatSLP_swap_row)
  t = a[s.row]
  a[s.row] = a[s.col]
  a[s.col] = t
end

function mul_t!{T}(a::Array{T, 1}, s::SmatSLP_add_row{T})
  a[s.col] = a[s.col]*s.val + a[s.row]
end

function apply!{T}(a::Array{T}, b::Array{SmatSLP_add_row{T}, 1})
  for i=length(b):-1:1
    mul!(a, b[i])
  end
  return a
end

function apply_t!{T}(a::Array{T}, b::Array{SmatSLP_add_row{T}, 1})
  for i=1:length(b)
    mul_t!(a, b[i])
  end
  return a
end

function random_SmatSLP{T}(A::Smat{T}, i::Int, v::UnitRange)
  a = Array(SmatSLP_add_row{Int}, i)
  for j=1:i
    c = rand(v)
    while c==0
      c = rand(v)
    end
    i1 = rand(1:rows(A))
    i2 = rand(1:rows(A))
    while i1==i2
      i1 = rand(1:rows(A))
      i2 = rand(1:rows(A))
    end
    a[j] = SLP_AddRow(i1, i2, Int(c))
  end
  return a
end

@doc """
  valence_mc{T}(A::Smat{T}; extra_prime = 2, trans = Array{SmatSLP_add_row{T}, 1}()) -> T

  Uses a Monte-Carlo alorithm to  compute the valence of A. The valence is the vaence of the minimal polynomial f of A'*A, thus the last non-zero coefficient,
  typically f(0).

  The valence is computed modulo various primes until the computation
  stabilises for extra_prime many.

  trans, if given, is  a SLP (straight-line-program) in GL(n, Z). Then
  the valence of trans * A  is computed instead.
""" ->
function valence_mc{T}(A::Smat{T}; extra_prime = 2, trans = Array{SmatSLP_add_row{T}, 1}())
  # we work in At * A (or A * At) where we choose the smaller of the 2
  # matrices
  if false && cols(A) > rows(A)
    At = A
    A = transpose(A)
  else
    At = transpose(A)
  end
  if abs_max(A) > 2^20
    mm = mul_mod_big!
    println("mul big case")
  else
    mm = mul_mod!
    println("mul small case")
  end
  c1 = Array(Int, cols(A))
  c2 = Array(Int, rows(A))

  for i=1:cols(A)
    c1[i] = Int(rand(-10:10))
  end

  c = copy(c1) ## need the starting value for the other primes as well

  p = next_prime(2^30)
  k = FiniteField(p)
  d = 10
  v = Array{typeof(k(1)), 1}()
  push!(v, k(c1[1]))
  while true
    while length(v) <= d
      mm(c2, A, c1, p)
      if length(trans) >0
        apply_t!(c2, trans)
        apply!(c2, trans)
      end
      mm(c1, At, c2, p)
      push!(v, k(c1[1]))
    end
    fl, f = berlekamp_massey(v)
    if !fl || degree(f) > div(d, 2)-1
      d += 10
      continue
    end
    df = degree(f)
    println("Poly degree is $df, dims $(rows(A)) x $(cols(A))")

    V = fmpz(leading_coefficient(f))
    pp = fmpz(p)

    v = Array(typeof(k(1)), 2*degree(f)+1)
    while true
      p = next_prime(p)
      println(p)
      k = FiniteField(p)
      copy!(c1, c)
      v[1] = k(c1[1])
      @time for i=1:2*degree(f)
        mm(c2, A, c1, p)
        if length(trans) >0
          apply_t!(c2, trans)
          apply!(c2, trans)
        end
        mm(c1, At, c2, p)
        v[i+1] = k(c1[1])
      end
      @time fl, f = berlekamp_massey(v)
      if !fl || degree(f) > df
        println("initial guess was wrong...")
        d += 10
        return -1
      end

      Vn = crt(V, pp, fmpz(leading_coefficient(f)), fmpz(p))
      pp *= p
      if 2*Vn > pp
        Vn = Vn - pp
      end
      if Vn == V
        extra_prime -= 1
        if extra_prime == 0
          return V
        end
      end
      V = Vn
    end
  end
end

################################################################################
#
# Basics: sub
#         vcat!
#         push! (add a SmatRow to it)
#         getindex -> A[i,j] for reading
#         rows, cols
#         copy
#
################################################################################

function sub{T}(A::Smat{T}, r::UnitRange, c::UnitRange)
  B = Smat{T}()
  B.nnz = 0
  for i=r
    rw = SmatRow{T}()
    ra = A.rows[i]
    for j=1:length(ra.values)
      if ra.pos[j] in c
        push!(rw.values, ra.values[j])
        push!(rw.pos, ra.pos[j]-c.start+1)
      end
    end
    if length(rw.pos)>0
      push!(B.rows, rw)
      B.nnz += length(rw.pos)
    end
  end
  B.r = length(r)
  B.c = length(c)
  return B
end

function copy{T}(A::Smat{T})
  return sub(A, 1:rows(A), 1:cols(A))
end

function rows(A::Smat)
  @assert A.r == length(A.rows)
  return A.r
end

function cols(A::Smat)
  return A.c
end

function getindex{T}(A::Smat{T}, i::Int, j::Int)
  if i in 1:A.r
    ra = A.rows[i]
    p = findfirst(x->x==j, ra.pos)
    if p!= 0
      return ra.values[p]
    end
  end
  return T(0)
end

function getindex{T}(A::Smat{T}, i::Int)
  if i in 1:A.r
    return A.rows[i]
  end
  return SmatRow{T}()
end

function vcat!{T}(A::Smat{T}, B::Smat{T})
  A.r += B.r
  A.c = max(A.c, B.c)
  A.nnz += B.nnz
  A.rows = vcat(A.rows, B.rows)
  @assert length(A.rows) == A.r
end

function push!{T}(A::Smat{T}, B::SmatRow{T})
  if length(B.pos)>0
    push!(A.rows, B)
    A.r += 1
    @assert length(A.rows) == A.r
    A.nnz += length(B.pos)
    A.c = max(A.c, B.pos[end])
  end
end

################################################################################
#
# convert to fmpz_mat
#
################################################################################
@doc """
  fmpz_mat{T <: Integer}(A::Smat{T})

  The same matix A, but as an fmpz_mat.
  Requires a conversion from type T to fmpz.
""" ->
function fmpz_mat{T <: Integer}(A::Smat{T})
  B = MatrixSpace(FlintZZ, A.r, A.c)()
  for i = 1:length(A.rows)
    ra = A.rows[i]
    for j = 1:length(ra.pos)
      B[i, ra.pos[j]] = ra.values[j]
    end
  end
  return B
end

@doc """
  fmpz_mat(A::Smat{fmpz})

  The same matix A, but as an fmpz_mat.
""" ->
function fmpz_mat(A::Smat{fmpz})
  B = MatrixSpace(FlintZZ, A.r, A.c)()
  for i = 1:length(A.rows)
    ra = A.rows[i]
    for j = 1:length(ra.pos)
      B[i, ra.pos[j]] = ra.values[j]
    end
  end
  return B
end

################################################################################
#
# convert to Julia
#
################################################################################
@doc """
  sparse{T}(A::Smat{T}) -> sparse{T}

  The same matrix, but as a julia-sparse matrix
""" ->
function sparse{T}(A::Smat{T})
  I = Array(Int, A.nnz)
  J = Array(Int, A.nnz)
  V = Array(T, A.nnz)
  i = 1
  for r = 1:rows(A)
    for j=1:length(A.rows[r].pos)
      I[i] = r
      J[i] = A.rows[r].pos[j]
      V[i] = A.rows[r].values[j]
      i += 1
    end
  end
  return sparse(I, J, V)
end

@doc """
  Array{T}(A::Smat{T}) -> Array{T, 2}

  The same matrix, but as a two-dimensional julia-Array.
""" ->
function Array{T}(A::Smat{T})
  R = Array(T, A.r, A.c)
  for i=1:rows(A)
    for j=1:length(A.rows[i].pos)
      R[i,A.rows[i].pos[j]] = A.rows[i].values[j]
    end
  end
  return R
end

################################################################################
#
# Echelonisation (transform to upper triangular)
#
# inplace operation, uses
#    add_scaled_row!
#    transform_row!
#
################################################################################

@doc """
  swap_rows!{T}(A::Smat{T}, i::Int, j::Int)

  swaps, inplace, the i-th row and the j-th
""" ->
function swap_rows!{T}(A::Smat{T}, i::Int, j::Int)
  a = A.rows[i]
  A.rows[i] = A.rows[j]
  A.rows[j] = a
end

# rows j -> row i*c + row j
@doc """
  add_scaled_row!{T}(A::Smat{T}, i::Int, j::Int, c::T)

  adds, inplace, the c*i-th row to the j-th
""" ->
function add_scaled_row!{T}(A::Smat{T}, i::Int, j::Int, c::T)
  sr = SmatRow{T}()
  pi = 1
  pj = 1
  @assert c != 0
  while pi <= length(A.rows[i].pos) && pj <= length(A.rows[j].pos)
    if A.rows[i].pos[pi] < A.rows[j].pos[pj]
      push!(sr.pos, A.rows[i].pos[pi])
      push!(sr.values, c*A.rows[i].values[pi])
      pi += 1
    elseif A.rows[i].pos[pi] > A.rows[j].pos[pj]
      push!(sr.pos, A.rows[j].pos[pj])
      push!(sr.values, A.rows[j].values[pj])
      pj += 1
    else
      n = c*A.rows[i].values[pi] + A.rows[j].values[pj]
      if n != 0
        push!(sr.pos, A.rows[i].pos[pi])
        push!(sr.values, n)
      end
      pi += 1
      pj += 1
    end
  end
  while pi <= length(A.rows[i].pos)
    push!(sr.pos, A.rows[i].pos[pi])
    push!(sr.values, c*A.rows[i].values[pi])
    pi += 1
  end
  while pj <= length(A.rows[j].pos)
    push!(sr.pos, A.rows[j].pos[pj])
    push!(sr.values, A.rows[j].values[pj])
    pj += 1
  end
  A.nnz = A.nnz - length(A.rows[j].pos) + length(sr.pos)
  A.rows[j] = sr
  return A
end

# row i -> a*row i + b * row j
# row j -> c*row i + d * row j
@doc """
  transform_row!{T}(A::Smat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)

  Inplace, replaces the i-th row and the j-th row by
  [a,b; c,d] * [i-th-row ; j-th row]
""" ->
function transform_row!{T}(A::Smat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)
  sr = SmatRow{T}()
  tr = SmatRow{T}()
  pi = 1
  pj = 1
  while pi <= length(A.rows[i].pos) && pj <= length(A.rows[j].pos)
    if A.rows[i].pos[pi] < A.rows[j].pos[pj]
      if a != 0
        push!(sr.pos, A.rows[i].pos[pi])
        push!(sr.values, a*A.rows[i].values[pi])
      end
      if c != 0
        push!(tr.pos, A.rows[i].pos[pi])
        push!(tr.values, c*A.rows[i].values[pi])
      end
      pi += 1
    elseif A.rows[i].pos[pi] > A.rows[j].pos[pj]
      if b != 0
        push!(sr.pos, A.rows[j].pos[pj])
        push!(sr.values, b*A.rows[j].values[pj])
      end
      if d != 0
        push!(tr.pos, A.rows[j].pos[pj])
        push!(tr.values, d*A.rows[j].values[pj])
      end
      pj += 1
    else
      m = a*A.rows[i].values[pi] + b*A.rows[j].values[pj]
      n = c*A.rows[i].values[pi] + d*A.rows[j].values[pj]
      if m != 0
        push!(sr.pos, A.rows[i].pos[pi])
        push!(sr.values, m)
      end
      if n != 0
        push!(tr.pos, A.rows[i].pos[pi])
        push!(tr.values, n)
      end
      pi += 1
      pj += 1
    end
  end
  while pi <= length(A.rows[i].pos)
    if a != 0
      push!(sr.pos, A.rows[i].pos[pi])
      push!(sr.values, a*A.rows[i].values[pi])
    end
    if c != 0
      push!(tr.pos, A.rows[i].pos[pi])
      push!(tr.values, c*A.rows[i].values[pi])
    end
    pi += 1
  end
  while pj <= length(A.rows[j].pos)
    if b != 0
      push!(sr.pos, A.rows[j].pos[pj])
      push!(sr.values, b*A.rows[j].values[pj])
    end
    if d != 0
      push!(tr.pos, A.rows[j].pos[pj])
      push!(tr.values, d*A.rows[j].values[pj])
    end
    pj += 1
  end
  A.nnz = A.nnz - length(A.rows[j].pos) + length(sr.pos) - length(A.rows[i].pos) + length(tr.pos)
  A.rows[i] = sr
  A.rows[j] = tr
  @assert i<j

  return A
end

function abs_max(A::Smat{Int})
  m = abs(A.rows[1].values[1])
  for i in A.rows
    for j in i.values
      m = max(m, abs(j))
    end
  end
  return m
end

function abs_max(A::Smat{BigInt})
  m = abs(A.rows[1].values[1])
  for i in A.rows
    for j in i.values
      if ccall((:__gmpz_cmpabs, :libgmp), Int, (Ptr{BigInt}, Ptr{BigInt}),
        &m, &j) < 0
        m = j
      end
    end
  end
  return abs(m)
end

@doc """
  abs_max(A::Smat{Int}) -> Int
  abs_max(A::Smat{BigInt}) -> BigInt
  abs_max(A::Smat{fmpz}) -> fmpz

  Finds the largest, in absolute value, entry of A
""" ->
function abs_max(A::Smat{fmpz})
  m = abs(A.rows[1].values[1])
  for i in A.rows
    for j in i.values
      if cmpabs(m, j) < 0
        m = j
      end
    end
  end
  return abs(m)
end

function max(A::Smat{Int})
  m = A.rows[1].values[1]
  for i in A.rows
    for j in i.values
      m = max(m, j)
    end
  end
  return m
end

function max(A::Smat{BigInt})
  m = A.rows[1].values[1]
  for i in A.rows
    for j in i.values
      if ccall((:__gmpz_cmp, :libgmp), Int, (Ptr{BigInt}, Ptr{BigInt}),
        &m, &j) < 0
        m = j
      end
    end
  end
  return m
end

@doc """
  max(A::Smat{Int}) -> Int
  max(A::Smat{BigInt}) -> BigInt
  max(A::Smat{fmpz}) -> fmpz

  Finds the largest entry of A
""" ->
function max(A::Smat{fmpz})
  m = A.rows[1].values[1]
  for i in A.rows
    for j in i.values
      if cmp(m, j) < 0
        m = j
      end
    end
  end
  return m
end

function min(A::Smat{Int})
  m = 0
  for i in A.rows
    for j in i.values
      m = min(m, j)
    end
  end
  return m
end

function min(A::Smat{BigInt})
  m = BigInt(0)
  for i in A.rows
    for j in i.values
      if ccall((:__gmpz_cmp, :libgmp), Int, (Ptr{BigInt}, Ptr{BigInt}),
        &m, &j) > 0
        m = j
      end
    end
  end
  return m
end

@doc """
  min(A::Smat{Int}) -> Int
  min(A::Smat{BigInt}) -> BigInt
  min(A::Smat{fmpz}) -> fmpz

  Finds the smallest entry of A
""" ->
function min(A::Smat{fmpz})
  m = fmpz(0)
  for i in A.rows
    for j in i.values
      if cmp(m, j) > 0
        m = j
      end
    end
  end
  return m
end

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

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

  Returns the number of bits necessary to represent a
""" ->
function nbits(a::Int)
  a==0 && return 0
  return Int(ceil(log(abs(a))/log(2)))
end

function one_step{T}(A::Smat{T}, sr = 1)
  i = sr
  assert(i>0)
  all_r = Array(Int, 0)
  min = A.c
  while i <= length(A.rows)
    @assert length(A.rows[i].pos)>0
    @assert A.rows[i].pos[1] >= sr
    if A.rows[i].pos[1] < min
      min = A.rows[i].pos[1]
      all_r = [i]
    elseif A.rows[i].pos[1] == min
      push!(all_r, i)
    end
    i += 1
  end

  if length(all_r) == 0
    return A.r+1
  end
#  println("found multiple $(length(all_r)) possibilities of weight $([length(A.rows[_i].pos) for _i in all_r])")
  #sorting
  sort!(all_r, lt=function(a,b) return length(A.rows[a].pos) < length(A.rows[b].pos)end)
#  println("found multiple $(length(all_r)) possibilities of weight $([length(A.rows[_i].pos) for _i in all_r])")
  # we need to have rows[all_r[1]] == sr, so that the new pivot will be in
  # row sr

  if !(sr in all_r)
    q = A.rows[sr]
    A.rows[sr] = A.rows[all_r[1]]
    A.rows[all_r[1]] = q
    all_r[1] = sr
  else
    p = findfirst(all_r, sr)
    if p > 1
      q = A.rows[sr]
      A.rows[sr] = A.rows[all_r[1]]
      A.rows[all_r[1]] = q
      all_r[p] = all_r[1]
      all_r[1] = sr
    end
  end
  if length(all_r) == 1
    return sr+1
  end

  for j=2:length(all_r)
    x = A.rows[sr].values[1]
    y = A.rows[all_r[j]].values[1]
    g = x
    @assert x!= 0
    if y % x == 0
      c = -div(y, x)
#      @assert x*c == -y
      add_scaled_row!(A, sr, all_r[j], c)
    else
      g,a,b = gcdx(x, y)
#      @assert g == a*x + b*y
      c = -div(y, g)
#      @assert y == -c*g
      d = div(x, g)
#      @assert x == d*g
      transform_row!(A, sr, all_r[j], a, b, c, d)
    end

#    @assert A.rows[sr].entry[1]valuesval == g
#    @assert A.rows[sr].entry[1].val != 0
#    if length(A.rows[all_r[j]].entry) == 0 ||
#             A.rows[all_r[j]].entry[1].col > min
#    else
#      println("bad", x, " and ", y, " in ", sr, ", ", all_r[j])
#      println("min: ", min)
#    end

#    @assert length(A.rows[all_r[j]].entry) == 0 ||
#             A.rows[all_r[j]].entry[1].col > min
#  println("in one step: ilog2(max) now ", nbits(max(A)), " j:", j, " length: ", length(all_r))
  end
  sort!(all_r)
  for j=length(all_r):-1:2
    if length(A.rows[all_r[j]].pos) == 0
      deleteat!(A.rows, all_r[j])
      A.r -= 1
    end
  end
  return sr+1
end

@doc """
  upper_triangular{T}(A::Smat{T}; mod = 0)

  Inplace: transform A into an upper triangular matrix. If mod
  is non-zero, reduce entries modulo mod during the computation.
""" ->
function upper_triangular{T}(A::Smat{T}; mod = 0)
  for i = 1:min(rows(A), cols(A))
    x = one_step(A, i)
#    println("after one step: ilog2(max) now ", nbits(max(A)))
    if x>A.r
      return
    end
    if A.nnz > (A.r-i) * (A.c-i) /2 || nbits(abs_max(A)) > 200
      #println("calling  at level ", i, " bits: ", nbits(abs_max(A)), "nnz: ", A.nnz)
      h = sub(A, i:A.r, i:A.c)
      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r
      h = fmpz_mat(h)
      #println("calling dense hnf on a ", rows(h), " by ", cols(h), " matrix")
      if mod==0
        h = hnf(h)
      else
        h = nmod_mat(mod, h)
        rref!(h)
        h = Array(lift(h))
      end
      h = Smat(h)
      for j in h.rows
        rw = SmatRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i-1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      return
    end
  end
  return
end

function sparsity{T}(A::Smat{T})
  return A.nnz/(A.r * A.c), nbits(abs_max(A))
end

#include("/home/claus/bigRel")

# New stuff from tommy

################################################################################
#
#  Application of a transformation on the left side of a sparse matrix
#
################################################################################

function apply_left!{T}(A::Smat{T}, t::TrafoSwap{T})
  swap_rows!(A, t.i, t.j)
  return nothing
end

function apply_left!{T}(A::Smat{T}, t::TrafoAddScaled{T})
  add_scaled_row!(A, t.i, t.j, t.s)
  return nothing
end

function apply_left!{T}(A::Smat{T}, t::TrafoParaAddScaled{T})
  transform_row!(A, t.i, t.j, t.a, t.b, t.c, t.d)
  return nothing
end

function apply_left!{T}(A::Smat{T}, t::TrafoDeleteZero{T})
  deleteat!(A.rows, t.i)
  A.r -= 1
end

function apply_left!{T, S}(A::Smat{T}, t::TrafoPartialDense{S})
  i = t.i
  h = sub(A, t.rows, t.cols)
  deleteat!(A.rows, t.rows)
  A.r -= length(t.rows)
  @assert length(A.rows) == A.r
  # make dense matrix
  hdense = typeof(t.U)(h)

  #println(t.U)

  hdense = t.U * hdense

  h = _Smat(hdense, R = parent(A.rows[1].values[1]), zerorows = true)

  for k in 1:length(h.rows)
    j = h.rows[k]
    rw = SmatRow{T}()
    for e in 1:length(j.pos)
      push!(rw.pos, j.pos[e] + i-1)
      push!(rw.values, j.values[e])
    end
    insert!(A.rows, i + k - 1, rw)
    A.r += 1
  end
  return nothing
end

################################################################################
#
#  Application of a transformation on the right hand side of a row
#
################################################################################

function apply_right!{T}(x::Array{T, 1}, t::TrafoAddScaled{T})
  x[t.i] = x[t.i] + x[t.j]*t.s
  return nothing
end

function apply_right!{T}(x::Array{T, 1}, t::TrafoSwap{T})
  r = x[t.i]
  x[t.i] = x[t.j]
  x[t.j] = r
  return nothing
end

function apply_right!{T}(x::Array{T, 1}, t::TrafoParaAddScaled{T})
  r = t.a * x[t.i] + t.c * x[t.j]
  s = t.b * x[t.i] + t.d * x[t.j]
  x[t.i] = r
  x[t.j] = s
  return nothing
end

function apply_right!{T}(x::Array{T, 1}, t::TrafoDeleteZero{T})
  # move ith position to the back
  for j in length(x):-1:t.i+1
    r = x[j]
    x[j] = x[j - 1]
    x[j - 1] = r
  end
end

function apply_right!{T}(x::Array{T, 1}, t::TrafoPartialDense)
  s = MatrixSpace(parent(x[1]), 1, rows(t.U))(x[t.rows])
  #println("s :$s")
  s = s*t.U
  for (i, j) in zip(t.rows,1:rows(t.U))
    x[i] = s[1, j]
  end
  return nothing
end

################################################################################
#
#  One echelonization step w/o transformations
#
################################################################################

function _one_step{T}(A::Smat{T}, sr = 1)
  nr, trafo = _one_step_with_trafo(A, sr)
  return nr
end

function _one_step_with_trafo{T}(A::Smat{T}, sr = 1)
  trafos = Trafo[]
  i = sr
  assert(i>0)
  all_r = Array(Int, 0)
  min = A.c

  while i <= length(A.rows)
    @assert length(A.rows[i].pos)>0
    @assert A.rows[i].pos[1] >= sr
    if A.rows[i].pos[1] < min
      min = A.rows[i].pos[1]
      all_r = [i]
    elseif A.rows[i].pos[1] == min
      push!(all_r, i)
    end
    i += 1
  end

  if length(all_r) == 0
    return A.r+1, trafos
  end

  #sorting
  sort!(all_r, lt = (a,b) -> (length(A.rows[a].pos) < length(A.rows[b].pos)))
  # we need to have rows[all_r[1]] == sr, so that the new pivot will be in
  # row sr

  if !(sr in all_r)
    push!(trafos, TrafoSwap{T}(sr, all_r[1]))
    q = A.rows[sr]
    A.rows[sr] = A.rows[all_r[1]]
    A.rows[all_r[1]] = q
    all_r[1] = sr
  else
    p = findfirst(all_r, sr)
    if p > 1
      push!(trafos, TrafoSwap{T}(sr, all_r[1]))

      q = A.rows[sr]
      A.rows[sr] = A.rows[all_r[1]]
      A.rows[all_r[1]] = q

      all_r[p] = all_r[1]
      all_r[1] = sr
    end
  end

  if length(all_r) == 1
    return sr + 1, trafos
  end

  for j = 2:length(all_r)
    x = A.rows[sr].values[1]
    y = A.rows[all_r[j]].values[1]
    t = echelonize_basecase!(x, y, sr, all_r[j], A)
    push!(trafos, t)
  end

  sort!(all_r)

  for j = length(all_r):-1:2
    if length(A.rows[all_r[j]].pos) == 0
      deleteat!(A.rows, all_r[j])
      push!(trafos, TrafoDeleteZero{T}(all_r[j]))
      A.r -= 1
    end
  end
  return sr + 1, trafos
end

function nmod_mat(A::Smat{UIntMod})
  R = parent(A.rows[1].values[1])
  B = nmod_mat(A.r, A.c, R.mod.n)
  B.parent = NmodMatSpace(ResidueRing(FlintZZ, R.mod.n), A.r, A.c)

  for i = 1:length(A.rows)
    ra = A.rows[i]
    for j = 1:length(ra.pos)
      # alternatively, @inline the set_entry function in Nemo/src/flint/nmod_mat
      ccall((:nmod_mat_set_entry, :libflint), Void,
          (Ptr{nmod_mat}, Int, Int, UInt), &B, i - 1, ra.pos[j] - 1, ra.values[j].m)

      #set_entry!(B, i, ra.pos[j], ra.values[j].m)
      # B[i, ra.pos[j]] = ra.values[j].m
    end
  end
  return B
end

################################################################################
#
#  Echelonization via dense matrix
#
################################################################################

# Tommy: They do not remove zero rows! Don't let them remove zero rows

function echelonize_via_dense(h::Smat{UIntMod})
  R = parent(h.rows[1].values[1])
  # turn into dense structure
  hdense = nmod_mat(h)
  # echelonize
  # I need to the following trick since we do not have the transformation
  rref!(hdense)
  # put back into sparse matrix over R (!= base_ring h)
  return _Smat(hdense, R = R)
end

function echelonize_via_dense_with_trafo(h::Smat{fmpz})
  hdense = fmpz_mat(h)
  # echelonize
  hdense, U = hnf_with_transform(hdense)
  # put back into sparse matrix
  h = _Smat(hdense, zerorows = true)
  return h, U
end

function echelonize_via_dense(h::Smat{fmpz})
  hdense = fmpz_mat(h)
  # echelonize
  hdense = hnf(hdense)
  # put back into sparse matrix
  h = _Smat(hdense)
  return h
end

################################################################################
#
#  Base case of echelonization
#
################################################################################

# Given two rows i, j with entries x, y (x != 0) in the same column, this
# function must produce a zero at entry y.

function echelonize_basecase!{T}(x::T, y::T, i::Int, j::Int, A::Hecke.Smat{T})
  c = -divexact(y, x)
  add_scaled_row!(A, i, j, c)
  return TrafoAddScaled(i, j, c)
end

function echelonize_basecase!(x::fmpz, y::fmpz, i::Int, j::Int,
                              A::Hecke.Smat{fmpz})
  if y % x == 0
    c = -div(y, x)
    add_scaled_row!(A, i, j, c)
    return TrafoAddScaled(i, j, c)
  else
    g,a,b = gcdx(x, y)
    c = -div(y, g)
    d = div(x, g)
    Hecke.transform_row!(A, i, j, a, b, c, d)
    return TrafoParaAddScaled(i, j, a, b, c, d)
  end
end

################################################################################
#
#  Echelonization for sparse matrices
#
################################################################################

# there is some density limit at level i, at which dense echelonization
# will be invoked
function upper_triangular_with_trafo!(M::Smat, density_limit::Float64 = 0.5)
  f = (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit)
  return _upper_triangular_with_trafo!(M, f)
end

function upper_triangular_with_trafo!(M::Smat{fmpz},
                                      density_limit::Float64 = 0.5,
                                      size_limit::Int = 200)
  f =  (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit || nbits(abs_max(A)) > size_limit)
  return _upper_triangular_with_trafo!(M, f)
end

# is_dense_enough is a function (::Smat{T}, i::Int) -> Bool
# At each level i, is_dense_enough(A, i) is called.
# If it evaluates to true, then dense echelonization will be called.
function _upper_triangular_with_trafo!{T}(A::Smat{T}, is_dense_enough::Function)
  trafo = Trafo[]

  for i = 1:min(rows(A), cols(A))
    x, t = _one_step_with_trafo(A, i)
    append!(trafo, t)

    if x > A.r
      return trafo
    end

    if is_dense_enough(A, i)
      h = sub(A, i:A.r, i:A.c)

      h, U = echelonize_via_dense_with_trafo(h)

      # h will have zero rows

      push!(trafo, TrafoPartialDense(i, i:A.r, i:A.c, U))

      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r

      # append all rows of h to A
      # including zero rows

      for k in 1:length(h.rows)
        j = h.rows[k]
        if length(j.pos) == 0
          push!(A.rows, SmatRow{T}())
          A.r += 1
          continue
        end
        rw = SmatRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i-1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end

      k = A.r

      # Now delete trailing zero rows.
      while length(A.rows[k].pos) == 0
        deleteat!(A.rows, k)
        A.r -= 1
        push!(trafo, TrafoDeleteZero{T}(k))
        k -= 1
      end

      return trafo
    end
  end
  return trafo
end

# now the same without transformation

function upper_triangular!(M::Smat, density_limit::Float64 = 0.5)
  f = (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit)
  return _upper_triangular!(M, f)
end

function upper_triangular!(M::Smat{fmpz}, density_limit::Float64 = 0.5, size_limit::Int = 200)
  f =  (A, i) -> (A.nnz > (A.r-i) * (A.c-i) * density_limit || nbits(abs_max(A)) > size_limit)
  return _upper_triangular!(M, f)
end

function _upper_triangular!{T}(A::Smat{T}, is_dense_enough)
  for i = 1:min(rows(A), cols(A))
    x = _one_step(A, i)

    if x>A.r
      return nothing
    end

    if is_dense_enough(A, i)

      h = sub(A, i:A.r, i:A.c)
      deleteat!(A.rows, i:A.r)
      A.r -= length(i:A.r)
      @assert length(A.rows) == A.r
      # make dense matrix
      h = echelonize_via_dense(h)

      for j in h.rows
        rw = SmatRow{T}()
        for e in 1:length(j.pos)
          push!(rw.pos, j.pos[e] + i-1)
          push!(rw.values, j.values[e])
        end
        push!(A.rows, rw)
        A.r += 1
      end
      return nothing
    end
  end
  return nothing
end


text:
using CSV, DataFrames, wgregseq, CairoMakie

##
# Open Gene information
dir = @__DIR__
home_dir = joinpath(split(dir, '/')[1:end-2]...)
s = open("/$home_dir/data/ecocyc/genes.dat") do file
    read(file, String)
end

# Drop Comments
gene_list = split(s, "//")[3:end];
# Separate genes
gene_list = [split(x, '\n') for x in gene_list]
# Separate information
gene_list = [[occursin(" - ", x) ? split(x, " - ") : SubString{String}[] for x in gene] for gene in gene_list]
# Drop empty entries
gene_list = [gene[.~isempty.(gene)] for gene in gene_list]

##
# Make DataFrame

ID_list = String[]
name_list = String[]
TU_list = []
direction_list = String[]
position_list = Float64[]
synonym_list = []
accession_list = String[]

for x in gene_list
    name = filter(x -> x[1] == "UNIQUE-ID", x)
    if ~isempty(name)
        push!(ID_list, name[1][2])
    else
        push!(ID_list, "None")
    end

    name = filter(x -> x[1] == "COMMON-NAME", x)
    if ~isempty(name)
        push!(name_list, name[1][2])
    else
        push!(name_list, "None")
    end

    synonyms = filter(x -> x[1] == "SYNONYMS", x)
    _synonym_list = []
    if ~isempty(synonyms)
        for synonym in synonyms
            push!(_synonym_list, synonym[2])
        end
    else
        push!(_synonym_list, "none")
    end
    push!(synonym_list, _synonym_list)

    direction = filter(x -> x[1] == "TRANSCRIPTION-DIRECTION", x)
    if ~isempty(direction)
        push!(direction_list, direction[1][2])
    else
        push!(direction_list, "none")
    end

    accession = filter(x -> x[1] == "ACCESSION-1", x)
    if ~isempty(accession)
        push!(accession_list, accession[1][2])
    else
        push!(accession_list, "none")
    end

    components = filter(x -> x[1] == "COMPONENT-OF", x)
    tu_list_gene = []
    if ~isempty(components)
        for component in components
            if occursin("TU", component[2])
                push!(tu_list_gene, component[2])
            end
        end
    else
        push!(tu_list_gene, "none")
    end

    push!(TU_list, tu_list_gene)
    if direction_list[end] == "+"
        position = filter(x -> x[1] == "LEFT-END-POSITION", x)
        if ~isempty(position)
            push!(position_list, parse(Float64, position[1][2]))
        else
            push!(position_list, NaN)
        end
    elseif direction_list[end] == "-"
        position = filter(x -> x[1] == "RIGHT-END-POSITION", x)
        if ~isempty(position)
            push!(position_list, parse(Float64, position[1][2]))
        else
            push!(position_list, NaN)
        end
    else
        push!(position_list, NaN)
    end
end
df_genes = DataFrames.DataFrame(
    ID=ID_list,
    gene=name_list,
    direction=direction_list,
    transcription_units=TU_list,
    gene_position=position_list,
    synonyms=synonym_list,
    accession=accession_list
    )
CSV.write("/$home_dir/data/all_genes_table.csv", df_genes)
##
# Open Transcription Units
s = open("/$home_dir/data/ecocyc/transunits.dat") do file
    read(file, String)
end

# Drop Comments
tu_list = split(s, "//")[3:end]
# Separate promoters
tu_list = [split(x, '\n') for x in tu_list]
# Separate information
tu_list = [[occursin(" - ", x) ? split(x, " - ") : SubString{String}[] for x in tu] for tu in tu_list]
# Drop empty entries
tu_list = [tu[.~isempty.(tu)] for tu in tu_list]

##
# Write to DataFrame
ID_list = String[]
promoter_list = []

for x in tu_list
    name = filter(x -> x[1] == "UNIQUE-ID", x)
    if ~isempty(name)
        push!(ID_list, name[1][2])
    else
        push!(ID_list, "None")
    end

    components = filter(x -> ((x[1] == "COMPONENTS") && (occursin("PM", x[2]))), x)
    if ~isempty(components)
        for component in components
            if occursin("PM", component[2])
                push!(promoter_list, component[2])
            end
        end
    else
        push!(promoter_list, "none")
    end
end

df_tu = DataFrames.DataFrame(TU_ID=ID_list, promoter_ID=promoter_list)

tu_genes = Vector{String}[]
tu_gene_positions = Vector{Float64}[]
tu_direction = String[]
# Combine Genes and Transcription Units
for tu in df_tu.TU_ID
    _df = df_genes[map(x -> tu in x["transcription_units"], eachrow(df_genes)), ["gene", "direction", "gene_position"]]
    push!(tu_genes, _df.gene)
    push!(tu_gene_positions, _df.gene_position)
    direction = _df.direction |> unique
    if length(direction) > 1
        throw(ErrorException("Found more than one direction for transcription unit $tu."))
    elseif length(direction) == 0
        push!(tu_direction, "")
    else
        push!(tu_direction, direction[1])
    end
end
insertcols!(df_tu, 3, :genes=>tu_genes)
insertcols!(df_tu, 4, :direction=>tu_direction)
insertcols!(df_tu, 5, :gene_position=>tu_gene_positions)
##
# Open Promoter

s = open("/$home_dir/data/ecocyc/promoters.dat") do file
    read(file, String)
end

# Drop Comments
promoter_list = split(s, "//")[3:end]
# Separate promoters
promoter_list = [split(x, '\n') for x in promoter_list]
# Separate information
promoter_list = [[occursin(" - ", x) ? split(x, " - ") : SubString{String}[] for x in promoter] for promoter in promoter_list]
# Drop empty entries
promoter_list = [promoter[.~isempty.(promoter)] for promoter in promoter_list]

##
# Make DataFrame
ID_list = String[]
name_list = String[]
TSS_list = Float64[]
evidence_list = Vector{String}[]

for x in promoter_list
    name = filter(x -> x[1] == "UNIQUE-ID", x)
    if ~isempty(name)
        push!(ID_list, name[1][2])
    else
        push!(ID_list, "None")
    end

    name = filter(x -> x[1] == "COMMON-NAME", x)
    if ~isempty(name)
        push!(name_list, name[1][2])
    else
        push!(name_list, "None")
    end

    TSS = filter(x -> x[1] == "ABSOLUTE-PLUS-1-POS", x)
    if ~isempty(TSS)
        push!(TSS_list, parse(Float64, TSS[1][2]))
    else
        push!(TSS_list, NaN)
    end

    cits = filter(x -> x[1] == "CITATIONS", x)
    if ~isempty(cits)
        cit_list = String[]
        for cit in cits
            if occursin("EV-EXP", cit[2])
                push!(cit_list, "EXP")
            elseif occursin("EV-COMP", cit[2])
                push!(cit_list, "COMP")
            end
            if length(cit_list) == 0
                push!(cit_list, "none")
            end
        end
        push!(evidence_list, cit_list)
    else
        push!(evidence_list, String["none"])
    end

end
df_tss = DataFrames.DataFrame(promoter_ID=ID_list, promoter=name_list, tss=TSS_list, evidence=evidence_list)

df_tss
"rspAp" in name_list
##
# Join DataFrames
df_joint = outerjoin(df_tu, df_tss, on = :promoter_ID) |> unique

# Replace missing values
df_joint.TU_ID = coalesce.(df_joint.TU_ID, "None")
df_joint.genes = coalesce.(df_joint.genes, [["None"]])
df_joint.direction = coalesce.(df_joint.direction, "0")
df_joint.tss = coalesce.(df_joint.tss, NaN)
df_joint.evidence = coalesce.(df_joint.evidence, [["None"]])
println("First 20 rows of list:")
println(first(df_joint, 20))
df_joint
##
gdf = groupby(df_joint, "promoter_ID")
temp_df = DataFrame()
for _df in gdf
    ind = argmax(maximum(length.(_df.genes)))
    append!(temp_df, _df[ind:ind, :])
end
df_joint = temp_df

# Split DataFrame into promoters with TUs and without
df_joint_prom = df_joint[(df_joint.direction .!= "0") .& (.~ isnan.(df_joint.tss)), :]
CSV.write("/$home_dir/data/promoter_list_ecocyc.csv", df_joint_prom[:, ["promoter", "genes", "gene_position", "direction", "tss", "evidence"]])
##
genes_w_promoters = unique(vcat(df_joint_prom.genes...))
genes_wo_promoters = filter(x -> x ∉ genes_w_promoters, df_genes.gene)

df_no_prom = DataFrames.DataFrame()
for gene in genes_wo_promoters
    append!(df_no_prom, df_tu[map(x -> gene in x, df_tu.genes), :])
end
df_no_prom
CSV.write("/$home_dir/data/operons_without_promoters.csv", df_no_prom[:, ["genes", "direction", "gene_position"]])

ext: jl


text:
#--- seminaive formula uses partitions and reccurence, but does not use blocks

"""

  part(n::Int)

Returns Vector{Vector{Vector}} that includes all partitions of set [1, 2, ..., m]
into subests of size > 1 and < m

```jldoctest
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> part(4)
3-element Array{Array{Array{Int64,1},1},1}:
 Array{Int64,1}[[1,2],[3,4]]
 Array{Int64,1}[[1,3],[2,4]]
 Array{Int64,1}[[1,4],[2,3]]
```
"""
function part(m::Int)
    parts = Vector{Vector{Int}}[]
    for part in partitions(1:m)
      subsetslen = map(length, part)
      if !((1 in subsetslen) | (m in subsetslen))
        @inbounds push!(parts, part)
      end
    end
    parts
end

"""
    pyramidmoment(data::Matrix, m::Int)

Returns Array{Float, m}, the m'th moment tensor

```jldoctest
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> pyramidmoment(M, 3)
2×2×2 Array{Float64,3}:
[:, :, 1] =
 -0.523092   0.142552
  0.142552  -0.0407653

[:, :, 2] =
  0.142552   -0.0407653
 -0.0407653   0.0120729
```
"""
function pyramidmoment(data::Matrix{T}, m::Int) where T<: AbstractFloat
    n = size(data,2)
    ret = zeros(T, fill(n, m)...)
    for ind in pyramidindices(m, n)
      @inbounds temp = momel(data, ind)
      for per in collect(permutations([ind...]))
        @inbounds ret[per...] = temp
      end
    end
    ret
end

"""
    mixedel(cum::Vector{Array{T}}, mulind::Tuple, parts::Vector{Vector{Vector{Int}}})

Returns Float, element of the sum of products of lower cumulants at given
 multi-index and set of its partitions
"""
function mixedel(cum::Vector{Array{T}}, mulind::Tuple,
                  parts::Vector{Vector{Vector{Int}}}) where T<: AbstractFloat
    sum = 0.
    for k = 1:length(parts)
        prod = 1.
        for el in parts[k]
            @inbounds prod*= cum[size(el,1)][mulind[el]...]
        end
        sum += prod
    end
    sum
end

"""

  mixedarr(cumulants::Vector{Array}, m::Int)

Returns Array{Float, m}, the mixed array (sum of products of lower cumulants)
"""
function mixedarr(cumulants::Vector{Array{T}}, m::Int) where T<: AbstractFloat
    n = size(cumulants[1], 1)
    sumofprod = zeros(fill(n, m)...)
    parts = part(m)
    for ind in pyramidindices(m, n)
      pyramid = mixedel(cumulants, ind, parts)
      for per in collect(permutations([ind...]))
        @inbounds sumofprod[per...] = pyramid
      end
    end
    sumofprod
end

"""

  pyramidcumulants(X::Matrix{Float}, m::Int)

Returns [Array{Float, 2}, ..., Array{Float, m}], vector co cumulants tensors of order
2, .., m

```
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> pyramidcumulants(M, 3)[2]
2×2×2 Array{Float64,3}:
[:, :, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0
```
"""
function pyramidcumulants(X::Matrix{T}, m::Int = 4) where T<: AbstractFloat
    cumulants = Array{T}[]
    push!(cumulants, pyramidmoment(X, 1))
    X = X .- mean(X, dims=1)
    for i in 2:m
      if i < 4
        push!(cumulants, pyramidmoment(X, i))
      else
        cumulant = pyramidmoment(X, i) - mixedarr(cumulants, i)
        push!(cumulants, cumulant)
      end
    end
    cumulants
end


text:
module CT
    #using Pkg
    #Pkg.add("GenericSVD")
    #Pkg.add("JLD")
    #Pkg.add("LinearAlgebra")
    #Pkg.add("Plots")
    #Pkg.add("NLopt")
    #Pkg.add("Random")
    #Pkg.add("JuMP")
    #Pkg.add("Printf")
    #Pkg.add("FileIO")
    #Pkg.add("ImageMagick")
    #Pkg.add("Images")
    #Pkg.add("FixedPointNumbers")
    #Pkg.add("ImageDraw")
    using GenericSVD
    using JLD
    using LinearAlgebra
    using Plots
    using NLopt
    using Random
    using JuMP
    using Printf
    using FileIO
    using ImageMagick
    using Images
    using FixedPointNumbers
    using ImageDraw
    # AB 1
    include("ueberpruefeAbsorptionskoeffizienten2x2.jl");
    include("ueberpruefeLGS.jl")
    include("ueberpruefeMatrixVektor2x2.jl");
    include("ueberpruefeEindeutigkeit2x2.jl");
    include("ueberpruefeMatrixVektora.jl");
    # AB 2
    include("Dateneinlesen.jl");
    include("Startpunkte.jl");
    include("distance.jl");
    include("getij.jl");
    include("getwallpositions.jl");
    include("createStartingPositionsmein.jl");
    include("ueberpruefeStartpunkt.jl");
    include("ueberpruefeberechneteSchnittpunkte.jl");
    include("Schnittpunkte.jl");
    include("Pixelberechner.jl");
    include("almostequal.jl");
    include("ueberpruefePixelkanten.jl");
    include("ueberpruefePara.jl");
    include("berechneallgemeineSchnittpunkte.jl");
    include("resize.jl");
    include("sortieren.jl");
    include("finde.jl");
    include("ueberpruefeaundb.jl");
    include("ueberpruefeZeile.jl");
    include("ueberpruefePixel.jl");
    include("Pixelberechner.jl");
    include("ueberpruefeAbstand.jl");
    include("ueberpruefeMittelpunkte.jl");
    include("bestimmePixelnummer.jl");
    include("ueberpruefePixelStrahl.jl");
    include("erstelleMatrix.jl");
    include("ueberpruefeSatzvervollstaendigung1.jl");
    include("ueberpruefeLaengeStrahl.jl");
    include("ueberpruefeLaengePixel.jl");
    include("ueberpruefeMatrix.jl")
    # AB 3
    include("ueberpruefeeindeutigeLoesung2x2.jl")
    include("ueberpruefeeindeutigeLoesung3x2.jl")
    include("ueberpruefekeineLoesung.jl")
    include("ueberpruefenormalenVektor.jl")
    include("ueberpruefeParameter.jl")
    include("ueberpruefeSatzvervollstaendigung2.jl")
    include("ueberpruefenaechsterPunkt.jl")
    include("echteCTDaten.jl");
    include("findenaechstesBild.jl")
    include("Bilderrekonstruktion.jl")
    include("createStartingPositionsmein2.jl")
    include("distance2.jl")
    include("getij2.jl")
    include("getwallpositions2.jl")
    include("erstelleA.jl")
    include("genTikhonov.jl")
    # Zusatz AB 3
    include("ueberpruefeFormel.jl");
    include("createStartingPositionsmeinsneu.jl");
    include("MartinStrahlen.jl");
    include("ueberpruefeSatzvervollstaendigung.jl");
    include("ueberpruefeneueFormel.jl");
    include("ueberpruefeAbleitung.jl");
    include("ueberpruefef_nah.jl");
    include("ueberpruefeerror.jl");
    include("erzeugeFehlerbehafteteDaten.jl");
    include("ueberpruefeSchuld.jl");
    include("ueberpruefep.jl");
    include("ueberpruefef_nahfehler.jl")
    include("ueberpruefef_nahfehler.jl")
end


text:
using Test
using ModelConstructors, SMC, DSGE
using LinearAlgebra, PDMats, Distributions, SparseArrays
using Printf, Distributed, Random, HDF5, FileIO, JLD2

import SparseArrays.SparseMatrixCSC
import ModelConstructors.Setting

my_tests = [
            "smc",
            "helpers",
            "initialization",
            "resample",
            "util",
            "mutation"
            ]

for test in my_tests
    test_file = string("$test.jl")
    @printf " * %s\n" test_file
    include(test_file)
end


text:
# Run the benchmarks for Julia (run.jl) and Python (run.py). Copy the results in the vectors below and plot

using VegaLite
using DataFrames

# Agents.jl benchmark results
jlresults = [
0.124292399,
0.261890796,
0.382272875,
0.526156528,
0.643057459,
0.795138978,
1.012378349,
1.180964479,
1.248486157,
1.369958317]

# Mesa benchmark results
pyresults = [
0.9389688929659314,
1.9137998179649003,
2.9129630780080333,
4.070444974990096,
5.039975131978281,
6.156719513994176,
7.4113737950101495,
8.584513186011463,
10.217067367048003,
11.377028140006587
]

size_range = (100:100:1000) .* 100

dd = DataFrame(runtime=vcat(jlresults, pyresults), lang=vcat(["Agents.jl" for i in 1:10], ["Mesa" for i in 1:10]), gridsize=vcat(size_range, size_range));

@vlplot(data=dd,
  mark={:line,
    size=4,
    point={filled=false, fill="white", size=60}
  },
  x={:gridsize,
    type="ordinal",
    title="Grid size",
    axis={titleFontSize=16, labelFontSize=14},
    grid=true,
  },
  y={:runtime,
    type="quantitative",
    title="Run time (seconds)",
    axis={titleFontSize=13,
    labelFontSize=14},
    grid=true
  },
  color={:lang,
    type="nominal",
    title="Framework",
    legend={titleFontSize=15, labelFontSize=14},
    scale={range=["#EA98D2", "#659CCA"]},
  },
  height=400,
  width=500,
)

dd2 = DataFrame(ratio=pyresults./jlresults, gridsize=size_range);

@vlplot(data=dd2,
  mark={:bar},
  x={:gridsize,
    type="ordinal",
    title="Grid size",
    axis={titleFontSize=16, labelFontSize=14},
  },
  y={:ratio,
    type="quantitative",
    title="Mesa/Agents run time",
    axis={titleFontSize=13,
    labelFontSize=14}
  },
  height=400,
  width=500,
)


text:
module Impl
export Mode, checkmode

"""
    struct Mode
\n `Mode` could be seen as an extension of `Val` and is intended to be used as
a type for a function argument. Its semantics is a set of symbols (flags) whose
value is taken into account during method dispatch. A specialization of `Mode`
allows to declare a list of symbols which a function method would accept for
the argument. The symbols could be accompanied with additional parameters with
prescribed types.

Here is a more detailed description of the mechanics of the type. A specialized
`Mode` type **M=Mode{[s₁⇒t₁, s₂⇒t₂, ...]}** is determined by a collection of
symbols **s₁**, **s₂**, ...  of type `Symbol` and corresponding types **t₁**,
**t₂**, ... (of type `DataType` or `Union` of `DataType`s). An instance
**m::Mode** of `Mode` additionally contains values **vᵢ** for each **sᵢ** of
type **tᵢ**. Let **param(x)** denote the collection **{sᵢ=>tᵢ}ᵢ** of **x**
(here **x** is an instance or a specialized type of `Mode`). Then **m** is of
type **M** only if **param(m) ⊆ param(M)**. This allows the dispatch to choose
an appropriate method of a function based on **param(m)**.

See also: [`checkmode`](@ref)

# Constructor of the type specialization
    Mode[ s₁ [=> t₁] [, s₂ [=> t₂]]... ]
\n Construct a specialization **M=Mode{[s₁⇒t₁, s₂⇒t₂, ...]}** to be used as a
type for a argument in a function method declaration. The argument with type
**M** would accept only instances **m::Mode** with **param(m) ⊆ param(M)**.
Some or all **tᵢ** may absent in the type declaration which defaults to
`Nothing`. An example: `Mode[:a, :b => Int, :c => Tuple{Int, String}]`.

Symbol `==` could also be added as the first argument in a constructor call to indicate that the concrete type is wanted (not an `UnionAll` as above). This
option is added only for auxiliary purposes and generally should not be useful.

Note that such nonconventional syntax with brackets `[`,`]` is used to make the
code for a type declaration look similar to the presentation of the type name
in a function signature.

# Constructors of an instance
    Mode( s₁[=> v₁] [, s₂ [=> v₂]]... )
\n Construct an instance **m::Mode{[s₁⇒typeof(v₁), s₂⇒typeof(v₂), ...]}** with
values **v₁, v₂, ...**. Some or all of **vᵢ** might be omited which defaults to
`nothing`. An example: `Mode(:a => 25, :b, :c => "Hello world!")`.

    Mode(s)
\n Construct an instance **m::Mode{[s⇒Nothing]}** with `nothing` value. The
type and value could be set by a subsequent call **m**`=>`**v**. Several
instances could be joined with `~` operator. For example, `Mode(:a)=> 25 ~ Mode
(:b) ~ Mode(:c)=> "Hello world!"` is equivalent to `Mode(:a => 25, :b, :c =>
"Hello world!")`.

# Operations on an instance
Let **m, m₁, m₂::Mode**. Then
- `m => v` given `m::Mode{[s => Nothing]}`: return a new instance with symbol
  **s** having value and type of `v`.
- `m₁ ~ m₂`: join `m₁` and `m₂`. Throws an ArgumentError if `m₁` and `m₂`
  contain the same symbols.
- `keys(m), values(m), pairs(m)`: return symbols / values / pairs of symbols
  and values.
- `m[s₁ [, s₂]... ]` for `sᵢ::Symbol`: return the value of `s₁` / tuple of
  values of `s₁, s₂, ...`.
- `m[]`: if `m` contains only one value, return it; throw an ArgumentError
  otherwise
- `checkmode(m, ...)`: check if `m` is a `Mode` and contains the prescribed
  symbols, and, possibly, do an action. See [`ArgumentModes.checkmode`](@ref)
  for detailed description.

# Examples

### 1. Usage in function declarations and calls
```julia-repl
julia> f(x) = println("Anything: \$x")
       f(x::Mode[:iterator => Any]) = println("Values from iterator: \$(x[]...)")
       f(x::Mode[:fromargs], y...) = println("Fromargs: \$(y...)")
       function f(x::Union{Int, Mode[:a, :b, :c => Int]})
           checkmode(x, :a) do _;  println(":a")  end
           checkmode(x, :c) do c;  println(":c = \$c")  end
           checkmode(x, :b) && println(":b")
           if x isa Int;  println("x = \$x") end
       end
f (generic function with 4 methods)

julia> methods(f)
# 4 methods for generic function "f":
[1] f(x::Mode[:iterator => Any]) in Main at REPL[34]:1
[2] f(x::Mode[:fromargs], y...) in Main at REPL[36]:1
[3] f(x::Union{Int64, Mode[:c => Int64, :a, :b]}) in Main at REPL[38]:1
[4] f(x) in Main at REPL[33]:1

julia> f(25.0)
Anything: 25.0

julia> f(Mode(:iterator)=> 1:5)
Values from iterator: 12345

julia> f(Mode(:fromargs), 1, 2, 3, 4, 5)
Fromargs: 12345

julia> f(1)
x = 1

julia> f(Mode(:a, :c => 125))
:a
:c = 125

julia> f(Mode(:fromargs, :iterator => 1:5), 2)
ERROR: MethodError: no method matching f(::Mode[==, :iterator => UnitRange, :fromargs], ::Int64)
Closest candidates are:
  f(::Mode[:fromargs], ::Any...) at REPL[36]:1
  f(::Any) at REPL[33]:1
Stacktrace:
 [1] top-level scope
   @ REPL[59]:1

```

### 2. Construction of types and instances, operations on an instance
```julia-repl
julia> M = Mode[:a => Real, :b => Union{Number, String}, :c]
Mode[:a => Real, :b => Union{Number, String}, :c]

julia> m1 = Mode(:a => 2.5, :c)
Mode[==, :a => Float64, :c]((a = 2.5, c = nothing))

julia> m1 isa M
true

julia> m2 = Mode(:a)=> 2.5 ~ Mode(:c)
Mode[==, :a => Float64, :c]((a = 2.5, c = nothing))

julia> m1 === m2
true

julia> Mode(:a => 2.5, :c, :d) isa M
false

julia> Mode(:a => "2", :c) isa M
false

julia> Mode(:a => 2.5, :c) ~ Mode(:b)=> "Hello world!"
Mode[==, :a => Float64, :c, :b => String]((a = 2.5, b = "Hello world!", c = nothing))

julia> Mode(:a => 2.5, :c) ~ Mode(:b)=> "Hello world!" ~ Mode(:a)=> 2
ERROR: ArgumentError: Duplicated keys :a

julia> checkmode(m, :a)
true

julia> checkmode(m, &, :a, :b, :g)
false

julia> checkmode(m, |, :a, :b, :g)
true

julia> m[:a, :c]
(2.5, nothing)

julia> m2 = Mode(:a => 2.5)
Mode[:a => Float64]((a = 2.5,))

julia> m2[]
2.5
```

# Notes

- `show` is redefined for the `DataType` of `Mode` for more compact and
  pleasent presentation of **s⇒t** pairs declared in the type. The actual
  signature of the type is `struct ArgnmentModes.Mode{NameType,
  ValuesNamedTuple}` where `NameType` containes all pairs **s⇒t** and
  `ValuesNamedTuple` contains `NamedTuple` for storing the actual values
  associated with the symbols.
"""
struct Mode{NameType, ValuesNamedTuple}
  values::ValuesNamedTuple
end

@inline Mode() = Mode{Union{}, typeof((;))}((;))

# using @generated to reduce compilation/inference time
@generated function _Mode_instance_g(symbs, kwargs)
  #@nospecialize
  n = length(kwargs.parameters)
  types = Vector{DataType}(undef, n)
  for (i,e) in enumerate(kwargs.parameters)
    if e == Symbol;  types[i] = Nothing;
    else   types[i] = e.parameters[2]   end
  end
  symbs = Symbol[symbs.parameters[1]...]

  if n == 1;   si = [1]
  else
    si = sortperm(symbs)
    symbs = symbs[si]
    types = types[si]
  end

  for i in 2:length(symbs)
    if symbs[i-1] == symbs[i]
      dubls = unique(symbs[j]  for j in i:length(symbs)
        if symbs[j-1] == symbs[j])
      msg = "Duplicated keys $(join(QuoteNode.(dubls), ", "))"
      return quote
        $(Expr(:meta, :inline))
        throw(ArgumentError($msg))
      end
    end
  end

  tpls = [:(  Tuple{$(QuoteNode(symbs[i])),$(types[i])}  )  for i in 1:n]
  vals = [types[i] == Nothing ? :(nothing) : :(kwargs[$(si[i])][2])
    for i in 1:n]
  ntvals = [Expr(:kw, s, v)  for (s,v) in zip(symbs, vals)]
  qsymbs = [QuoteNode(s)  for s in symbs]

  if n == 1
    return quote
      $(Expr(:meta, :inline))
      Mode{
          $(tpls...),
          NamedTuple{($(qsymbs...),), Tuple{$(types...)}}
        }((; $(ntvals...)))
    end
  else
    return quote
      $(Expr(:meta, :inline))
      Mode{
          Union{$(tpls...)},
          NamedTuple{($(qsymbs...),), Tuple{$(types...)}}
        }((; $(ntvals...)))
    end
  end
end
@generated function _extract_symbols_from_args_g(kwargs)
  #@nospecialize
  vals = Expr[
    x == Symbol ? :(kwargs[$i]) : :(kwargs[$i][1])
    for (i,x) in enumerate(kwargs.parameters)
  ]
  quote
    $(Expr(:meta, :inline))
    Val{($(vals...),)}()
  end
end

@inline function Mode(kwargs::Union{<:Pair{Symbol}, Symbol}...)
  _Mode_instance_g(_extract_symbols_from_args_g(kwargs), kwargs)
end

@inline Base.:(=>)(::Mode{Tuple{Name, Nothing}}, val) where
  Name = Mode(Name => val)

@generated function _Combined_namedtuples_to_pairs_g(nt1, nt2)
  #@nospecialize
  n1 = length(nt1.names);   n2 = length(nt2.names)
  prs = Vector{Expr}(undef, n1+n2)
  for i in 1:n1;   prs[i] = :($(QuoteNode(nt1.names[i])) => nt1[$i])   end
  for i in 1:n2;   prs[n1+i] = :($(QuoteNode(nt2.names[i])) => nt2[$i])   end
  quote
    $(Expr(:meta, :inline))
    ($(prs...),)
  end
end

@inline function Base.:(~)(mode1::Mode, mode2::Mode)
  Mode(_Combined_namedtuples_to_pairs_g(mode1.values, mode2.values)...)
end

@generated function _getindex_g(symbs, args, concrete)
  #@nospecialize
  n = length(args.parameters)
  symbs = Symbol[symbs.parameters[1]...]
  si = sortperm(symbs)
  symbs = symbs[si]
  tpls = Vector{Expr}(undef, n)
  if concrete == Val{true}
    types = Vector{Union{Expr, Symbol}}(undef, n)
  end
  for i in 1:n
    oi = si[i]
    if args.parameters[oi] == Symbol
      tpls[i] = :(Tuple{$(QuoteNode(symbs[i])), Nothing})
      if concrete == Val{true};   types[i] = :(Nothing)   end
    else
      tpls[i] = :(Tuple{$(QuoteNode(symbs[i])), args[$oi][2]})
      if concrete == Val{true};   types[i] = :(args[$oi][2])   end
    end
  end

  if concrete == Val{true}
    symbs = QuoteNode[QuoteNode(x)  for x in symbs]
    t1 = n == 1 ? tpls[1] : :(Union{$(tpls...)})
    return quote
      $(Expr(:meta, :inline))
      Mode{$t1, NamedTuple{($(symbs...),), Tuple{$(types...)}}}
    end
  else
    t1 = n == 1 ? tpls[1] : :(Union{$(tpls...)})
    return quote
      $(Expr(:meta, :inline))
      Mode{<:$t1}
    end
  end
end
@inline Base.getindex(::Type{Mode}, args::Union{Pair{Symbol, <:Union{DataType,
    Union}}, Symbol}...) =
  _getindex_g(_extract_symbols_from_args_g(args), args, Val(false))
@inline Base.getindex(::Type{Mode}, ::typeof(==), args::Union{Pair{Symbol,
    <:Union{DataType, Union}}, Symbol}...) =
  _getindex_g(_extract_symbols_from_args_g(args), args, Val(true))
@inline Base.getindex(::Type{Mode}) = typeof(Mode())
@inline Base.getindex(::Type{Mode}, ::typeof(==)) = typeof(Mode())

@inline Base.keys(m::Mode) = keys(m.values)
@inline Base.values(m::Mode) = values(m.values)
@inline Base.pairs(m::Mode) = pairs(m.values)
@inline Base.getindex(m::Mode, k::Symbol, kk::Symbol...) =
  getindex.((m.values,), (k, kk...))
@inline Base.getindex(m::Mode, k::Symbol) = m.values[k]
@inline function Base.getindex(m::Mode)
  if length(m.values) == 1
    m.values[1]
  else
    throw(ArgumentError("Expected exactly one value, but got keys $(keys(m))"))
  end
end

function text_for_type_desc(io, tpls, concrete)
  Base.show_type_name(io, Base.typename(Mode))
  print(io, "[")
  if concrete;   print(io, "==, ")   end

  function f(x)
    print(io, ":", x[1])
    if x[2] != Nothing
      print(io, " => ")
      if x[2] isa Union
        print(io, "Union{")
        types = Base.uniontypes(x[2])
        Base.show_type_name(io, Base.typename(types[1]))
        for t in Iterators.drop(types, 1)
          print(io, ", ")
          Base.show_type_name(io, Base.typename(t))
        end
        print(io, "}")
      else
        Base.show_type_name(io, Base.typename(x[2]))
      end
    end
  end

  pairs = sort(
      [y[1] => y[2]  for y in (x.parameters  for x in tpls)  if length(y) == 2];
      by= x->x[1]
    )

  if length(pairs) > 0
    f(pairs[1])
    for e in Iterators.drop(pairs, 1)
      print(io, ", ")
      f(e)
    end
  end

  print(io, "]")
end

function extract_list_of_elements(m::Type{<:Mode})
  body = m;   while body isa UnionAll;   body = body.body   end
  args = body.parameters[1];   if args isa TypeVar;   args = args.ub   end
  Base.uniontypes(args)
end

function Base.show(io::IO, m::Type{<:Mode})
  if m == typeof(Mode())
    Base.show_type_name(io, Base.typename(Mode))
    print(io, "[]")
  elseif m == Mode
    Base.show_type_name(io, Base.typename(Mode))
  else
    elems = extract_list_of_elements(m)
    text_for_type_desc(io, elems, !(m isa UnionAll))
  end
  nothing
end

"""
    checkmode(m, s::Symbol)
    checkmode(m, [&, |, ==, or !], s₁::Symbol[, s₂:Symbol]...)
    checkmode(f, m, [&, |, ==, or !], s₁::Symbol[, s₂:Symbol]...)
\n Tests if `m` is a `Mode` and contains the prescribed symbols.

- In the 1st form the function tests for the presence of `s` in `m`.
- In the 2nd form the function tests either (for `&` or if any operation symbol
  is omited) that all symbols `s₁`, `s₂`, ... are in `m`, or (for `|`) that at
  least one of them in `m`, or (for `==`) that `m` contains exactly the
  prescribed collection of the symbols, or (for `!`) that `m` does not contain
  any of the prescribed symbols.
- In the 3rd form the function tests whether all of the symbols `s₁`, `s₂`, ...
  are in `m` (also if `&` is provided). If they are --- return `f(m[s₁], m[s₂],
  ...)`; if not --- return `nothing`. If `==` is also added in the function
  call, the test will be positive only when `m` contains all and only the
  prescribed symbols. If `!` is added in the function call, the test will be
  positive only when `m` contains no any of the prescribed symbols; if that is
  true, `f` is called with no arguments. Similarly for `|` -- if test is true,
  `f` is called with no arguments.

See also: [`ArgumentModes.Mode`](@ref)

# Examples
```julia-repl
julia> m = Mode(:a => 1, :b => 2.5, :c)
Mode[==, :b => Float64, :a => Int64, :c]((a = 1, b = 2.5, c = nothing))

julia> checkmode(m, :a)
true

julia> checkmode(m, :d)
false

julia> checkmode(m, |, :a, :d)
true

julia> checkmode(m, &, :a, :d)
false

julia> checkmode(m, ==, :a, :b)
false

julia> checkmode(m, ==, :a, :b, :c)
true

julia> checkmode(m, :a) do a;  @show a   end
a = 1
1

julia> checkmode(m, :a, :b) do a, b;  @show a, b   end
(a, b) = (1, 2.5)
(1, 2.5)

julia> checkmode(m, :a, :e) do x, y;  @show x, y   end

julia> checkmode(12, :a)
false
```
"""
function checkmode end
@inline checkmode(m, s::Symbol) = m isa Mode && s in keys(m.values)
@inline checkmode(m, ::typeof(&), s::Symbol...) = m isa Mode &&
  foldl(&, s .∈ (keys(m.values),))
@inline checkmode(m, s::Symbol, s1::Symbol...) = checkmode(m, &, s, s1...)
@inline checkmode(m, ::typeof(|), s::Symbol...) = m isa Mode &&
  foldl(|, s .∈ (keys(m.values),))
@inline checkmode(m, ::typeof(==), s::Symbol...) = m isa Mode &&
  foldl(&, s .∈ (keys(m.values),)) && foldl(&, keys(m.values) .∈ (s,))
@inline checkmode(m, ::typeof(!), s::Symbol...) = m isa Mode &&
  !foldl(|, s .∈ (keys(m.values),))
@inline function checkmode(f, m, s::Symbol)
  if checkmode(m, s)
    return f(m.values[s])
  else
    return nothing
  end
end
@inline function checkmode(f, m, ::typeof(&), s::Symbol, ss::Symbol...)
  if checkmode(m, &, s, ss...)
    return f(m.values[s], getindex.((m.values,), ss)...)
  else
    return nothing
  end
end
@inline checkmode(f, m, s::Symbol, ss::Symbol...) =
  checkmode(f, m, &, s, ss...)
@inline function checkmode(f, m, ::typeof(==), s::Symbol, ss::Symbol...)
  if checkmode(m, ==, s, ss...)
    return f(m.values[s], getindex.((m.values,), ss)...)
  else
    return nothing
  end
end
@inline function checkmode(f, m, ::typeof(!), s::Symbol, ss::Symbol...)
  if checkmode(m, !, s, ss...)
    return f()
  else
    return nothing
  end
end
@inline function checkmode(f, m, ::typeof(|), s::Symbol, ss::Symbol...)
  if checkmode(m, |, s, ss...)
    return f()
  else
    return nothing
  end
end

end # module


text:
function _http_get_call_with_url(url::String)::Some

    try

        # should we check if this string is formatted as a URL?
        if (occursin("https://", url) == false)
            throw(ArgumentError("url $(url) is not properly formatted"))
        end

        # ok, so we are going to make a HTTP GET call with the URL that was passed in -
        # we want to handle the errors on our own, so do NOT have HTTP.j throw an excpetion -
        response = HTTP.request("GET", url; status_exception = false)

        # return the body -
        return Some(String(response.body))
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function _process_polygon_response(model::Type{T},
    response::String)::Tuple where T<:AbstractPolygonEndpointModel

    # setup type handler map -> could we put this in a config file to register new handlers?
    type_handler_dict = Dict{Any,Function}()
    type_handler_dict[PolygonAggregatesEndpointModel] = _process_aggregates_polygon_call_response
    type_handler_dict[PolygonOptionsContractReferenceEndpoint] = _process_options_reference_call_response
    type_handler_dict[PolygonPreviousCloseEndpointModel] = _process_previous_close_polygon_call_response
    type_handler_dict[PolygonGroupedDailyEndpointModel] = _process_aggregates_polygon_call_response
    type_handler_dict[PolygonDailyOpenCloseEndpointModel] = _process_daily_open_close_call_response
    type_handler_dict[PolygonTickerNewsEndpointModel] = _process_ticker_news_call_response
    type_handler_dict[PolygonTickerDetailsEndpointModel] = _process_ticker_details_call_response
    type_handler_dict[PolygonMarketHolidaysEndpointModel] = _process_market_holidays_call_response #ycpan
    type_handler_dict[PolygonExchangesEndpointModel] = _process_exchanges_call_response #ycpan
    type_handler_dict[PolygonStockSplitsEndpointModel] = _process_stock_splits_call_response #ycpan
    type_handler_dict[PolygonMarketStatusEndpointModel] = _process_market_status_call_response #ycpan
    type_handler_dict[PolygonDividendsEndpointModel] = _process_dividends_call_response #ycpan
    type_handler_dict[PolygonTickersEndpointModel] = _process_tickers_call_response #ycpan
    type_handler_dict[PolygonConditionsEndpointModel] = _process_conditions_call_response #ycpan
    type_handler_dict[PolygonStockFinancialsEndpointModel] = _process_stock_financials_call_response #ycpan
    type_handler_dict[PolygonTickerTypesEndpointModel] = _process_ticker_types_call_response #ycpan

    # lookup the type -
    if (haskey(type_handler_dict, model) == true)
        handler_function = type_handler_dict[model]
        return handler_function(response);
    end

    # default -
    return nothing
end

function api(model::Type{T}, complete_url_string::String;
    handler::Function = _process_polygon_response)::Tuple where T<:AbstractPolygonEndpointModel

    # execute -
    result_string = _http_get_call_with_url(complete_url_string) |> check

    # process and return -
    return handler(model, result_string)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libCZI_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libCZI")
JLLWrappers.@generate_main_file("libCZI", UUID("1ce5c631-45c7-5a5d-932a-8ac578e282f3"))
end  # module libCZI_jll


text:
# NE - Run carbon cases
push!(LOAD_PATH, joinpath(ENV["REPO_PATH"], "Sesi.jl/src"))

using CSV
using JuMP
using IterTools
using JSON
using DataFrames

import Sesi # import vs using so I can reload this from the REPL

reload("Sesi")

include("jl_utils.jl")

save = false
verb = true

if verb
    println(repeat("#",80))
    println("Starting carbon analysis")
    println(repeat("#",80))
end

folderIn = joinpath(ENV["REPO_PATH"], "dataR1/jl_in/")

results = Dict()
prices = Dict()
costs = Dict()
reductions = Dict()
#fileNms = ["jlInput_x0.csv","jlInput_2018.csv", "jlInput_x2.csv"]
fileNms = ["jlInput_x0.csv", "jlInput_x1.csv", "jlInput_x2.csv",
            "jlInput_x4.csv","jlInput_x6.csv", "jlInput_2018.csv"]

for dc_ratio in [1., 0.5, 0.1]
if verb
println("Dealing with demand ratio " * string(dc_ratio))
end

folderOut = joinpath(ENV["REPO_PATH"],
                     "dataR1/jl_out/demand_" * string(dc_ratio))

for sc in fileNms
    prices[sc] = []
    costs[sc] = []
    reductions[sc] = []

    println("Dealing with " * sc)
    results[sc] = Dict()
    fileNm = folderIn * sc

    carbonPrice_lst = [0., 50., 75., 100., 150., 200., 300., 500., 1000.,
        10000., 100000.]
    for carbonPrice in carbonPrice_lst
        sesi, p, f = build_model(fileNm; verb=false, startDate=Date(2016,1,1),
            endDate=Date(2016,12,31), carbonPrice=carbonPrice, dc_ratio=dc_ratio)
        status = solve(sesi)
        println(status)
        results[sc][carbonPrice] = Sesi.campus_bill(sesi, f, p, verb=0)
        if verb
            @printf("carbon price: %d - %.2f\n",carbonPrice,
                carbon_cost(results, sc, 0., carbonPrice))
            @printf("carbon reduction: %d - %.2f\n",carbonPrice,
                carbon_reduction(results, sc, 0., carbonPrice))
            @printf("sanity : %d - %.2f\n",carbonPrice,
                sanity(results, sc, 0., carbonPrice))
        end
        push!(prices[sc], carbonPrice)
        push!(costs[sc], carbon_cost(results, sc, 0., carbonPrice))
        push!(reductions[sc], carbon_reduction(results, sc, 0., carbonPrice))
        if save
            Sesi.save_data(sesi, f, p, folderOut * sc[1:end-4] * "_"
                * string(carbonPrice))
        end
    end
end

json_data = json([results, prices, costs, reductions]);
open(folderOut * "ne_carbon_json_data", "w") do fw
    write(fw, json_data)
end
end

max_stars_repo_name: jdechalendar/sesi

max_issues_repo_name: jdechalendar/sesi

max_forks_repo_name: jdechalendar/sesi


text:
using Documenter, Clang

makedocs(
    modules = [Clang],
    clean = false,
    sitename = "Clang.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "API Reference" => "api.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaInterop/Clang.jl.git",
    target = "build",
)


text:
"""
    PDFLaTeX

The `pdflatex` engine.

Command line options `-output-format=dvi`, `-quiet`, and `-halt-on-error` are passed to the `pdflatex` executable, and both the input and output directories are also handled. If you would like to pass additional command line arguments, use the [`latexsvg`](@ref) function.
"""
struct PDFLaTeX <: LaTeXEngine end

function runlatex(
    ::PDFLaTeX,
    input_file::AbstractString,
    output_path::AbstractString;
    extra_args::Vector{String}=String[],
    runcommand::Bool=true
)
    command = `pdflatex $extra_args -output-format=dvi -output-directory="$output_path" -quiet -halt-on-error "$input_file"`
    if runcommand
        run(command)
    end
    return command
end

dviext(::PDFLaTeX) = "dvi"


text:
# This data structure is a doubly linked list (i.e. each node has a reference to the item before and the item after).
# The value lies in efficient insertion/deletion from the middle of the list, at the cost of more storage/complexity at each insert/delete.
# We represent this structure using:
#     HeadNode <--> first <--> second <--> ... <--> last <--> TailNode

type DLLNode{T}
  # typ::Type{T}
  ishead::Bool
  istail::Bool
  prevnode::DLLNode{T}
  nextnode::DLLNode{T}
  data::T

  DLLNode{T}(::Type{T}, ishead::Bool, istail::Bool) = new(ishead, istail)
  DLLNode{T}(ishead::Bool, istail::Bool, prevnode::DLLNode{T}, nextnode::DLLNode{T}, data::T) = new(ishead, istail, prevnode, nextnode, data)
end

Base.string{T}(x::DLLNode{T}) = "DLLNode{ishead=$(x.ishead), istail=$(x.istail), data=$(isdefined(x,:data) ? x.data : "")}"
Base.print{T}(io::IO, x::DLLNode{T}) = print(io, string(x))
Base.show{T}(io::IO, x::DLLNode{T}) = print(io, string(x))

type DoublyLinkedList{T}
  n::Int
  headnode::DLLNode{T}
  tailnode::DLLNode{T}

  # constructor creates an empty list with head/tail pointing to each other
  function DoublyLinkedList{T}(::Type{T})
    headnode = DLLNode{T}(T, true, false)
    tailnode = DLLNode{T}(T, false, true)
    headnode.nextnode = tailnode
    tailnode.prevnode = headnode
    new(0, headnode, tailnode)
  end
end

DoublyLinkedList{T}(::Type{T}) = DoublyLinkedList{T}(T)

head{T}(l::DoublyLinkedList{T}) = l.headnode
tail{T}(l::DoublyLinkedList{T}) = l.tailnode
peek{T}(l::SinglyLinkedList{T}) = l.headnode.nextnode.data

Base.first{T}(l::DoublyLinkedList{T}) = l.headnode.nextnode.data #note: error if empty!!
Base.last{T}(l::DoublyLinkedList{T}) = l.tailnode.prevnode.data #note: error if empty!!

# Note: iterating like "for x in X; (do something); end" translates to:
# state = start(X)
# while !done(X, state)
#   (x, state) = next(X, state)
#   (do something)
# end

# iterate with the syntax "for x in list"
Base.start{T}(l::DoublyLinkedList{T}) = l.headnode.nextnode # start with the first item (will be the tail if list is empty)
Base.done{T}(l::DoublyLinkedList{T}, state::DLLNode{T}) = state.istail
Base.next{T}(l::DoublyLinkedList{T}, state::DLLNode{T}) = (state, state.nextnode)

# allow backwards iteration with the syntax "for x in list.tailnode"
Base.start{T}(tailnode::DLLNode{T}) = tailnode.prevnode
Base.done{T}(tailnode::DLLNode{T}, state::DLLNode{T}) = state.ishead
Base.next{T}(tailnode::DLLNode{T}, state::DLLNode{T}) = (state, state.prevnode)

# insert data before node and return new node
function insertBefore!{T}(l::DoublyLinkedList{T}, node::DLLNode{T}, newnode::DLLNode{T})

  if node.ishead
    error("DoublyLinkedList: Can't insert before the head node")
  end

  before = node.prevnode
  node.prevnode = newnode
  before.nextnode = newnode
  newnode.nextnode = node
  newnode.prevnode = before
  l.n += 1
  newnode
end

# insert data before node and return new node
function insertBefore!{T}(l::DoublyLinkedList{T}, node::DLLNode{T}, data::T)

  if node.ishead
    error("DoublyLinkedList: Can't insert before the head node")
  end

  insertBefore!(l, node, DLLNode{T}(false, false, node.prevnode, node, data))
end

function insertAfter!{T}(l::DoublyLinkedList{T}, node::DLLNode{T}, newnode::DLLNode{T})
  if node.istail
    error("DoublyLinkedList: Can't insert after the tail node")
  end

  after = node.nextnode
  node.nextnode = newnode
  after.prevnode = newnode
  newnode.nextnode = after
  newnode.prevnode = node
  l.n += 1
  newnode
end

# insert data after node and return new node
function insertAfter!{T}(l::DoublyLinkedList{T}, node::DLLNode{T}, data::T)

  if node.istail
    error("DoublyLinkedList: Can't insert after the tail node")
  end

  insertAfter!(l, node, DLLNode{T}(false, false, node, node.nextnode, data))
end

function insertSorted!{T}(l::DoublyLinkedList{T}, data::T, compareFunction::Function, ascending::Bool)
  insertSorted!(l, DLLNode{T}(false, false, l.headnode, l.headnode, data), compareFunction, ascending)
end

function insertSorted!{T}(l::DoublyLinkedList{T}, newnode::DLLNode{T}, compareFunction::Function, ascending::Bool)

  # NOTE: it is expected that compareFunction is defined for 2 arguments of type T, and that
  # if the first arg is <, ==, > than the second, then it returns an int: neg, 0, pos

  for node in l
    # println("$compareFunction $newnode $node")
    cmpResult = compareFunction(newnode.data, node.data)

    if cmpResult != 0 && (cmpResult < 0) == ascending
      # add it here
      return insertBefore!(l, node, newnode)
    end
  end

  # add to end
  return insertBefore!(l, l.tailnode, newnode)
end

function Base.deleteat!{T}(l::DoublyLinkedList{T}, node::DLLNode{T})

  if node.ishead || node.istail
    return
  end

  node.prevnode.nextnode = node.nextnode
  node.nextnode.prevnode = node.prevnode
  l.n -= 1

  # clear this node's pointers
  node.prevnode = node
  node.nextnode = node
  return
end

function Base.deleteat!{T}(l::DoublyLinkedList{T}, data::T)
  warn("using slow version. l=$l data=$data")  #TODO:remove this?
  deleteat!(l, find(l, data))
end

function moveToEnd!{T}(l::DoublyLinkedList{T}, node::DLLNode{T})
  before = node.prevnode
  after = node.nextnode
  before.nextnode = after
  after.prevnode = before

  beforetail = l.tailnode.prevnode
  beforetail.nextnode = node
  node.prevnode = beforetail
  l.tailnode.prevnode = node
  node.nextnode = l.tailnode
end

# add to front or back and return the node
Base.push!{T}(l::DoublyLinkedList{T}, data::T) = insertBefore!(l, l.tailnode, data)
Base.unshift!{T}(l::DoublyLinkedList{T}, data::T) = insertAfter!(l, l.headnode, data)
Base.append!{T}(l::DoublyLinkedList{T}, data::T) = push!(l, data)
Base.prepend!{T}(l::DoublyLinkedList{T}, data::T) = unshift!(l, data)

# add to front or back and return the node
Base.push!{T}(l::DoublyLinkedList{T}, newnode::DLLNode{T}) = insertBefore!(l, l.tailnode, newnode)
Base.unshift!{T}(l::DoublyLinkedList{T}, newnode::DLLNode{T}) = insertAfter!(l, l.headnode, newnode)
Base.append!{T}(l::DoublyLinkedList{T}, newnode::DLLNode{T}) = push!(l, newnode)
Base.prepend!{T}(l::DoublyLinkedList{T}, newnode::DLLNode{T}) = unshift!(l, newnode)

Base.length(l::DoublyLinkedList) = l.n
Base.isempty(l::DoublyLinkedList) = length(l) == 0

function Base.filter{T}(f::Function, l::DoublyLinkedList{T})
  newl = DoublyLinkedList(T)
  for x in l
    if f(x.data)
      push!(newl, x.data)
    end
  end
  newl
end

function Base.find{T}(l::DoublyLinkedList{T}, data::T)
  for x in l
    if x.data == data
      return x
    end
  end
  return tail(l)
end

function Base.print{T}(io::IO, l::DoublyLinkedList{T})
  print(io, "{")
  print(io, join([x.data for x in l], ", "))
  print(io, "}")
end

function Base.show{T}(io::IO, l::DoublyLinkedList{T})
  print(io, "DoublyLinkedList(n=$(l.n))")
  print(io, l)
end

toarray{T}(l::DoublyLinkedList{T}) = T[node.data for node in l]


text:
using GAIO
using Test
using SafeTestsets

@testset "GAIO.jl" begin
    @safetestset "Box" begin
        include("box.jl")
    end
    @safetestset "BoxMap" begin
        include("boxmap.jl")
    end
    @safetestset "BoxSet" begin
        include("boxset.jl")
    end
    @safetestset "BoxPartition" begin
        include("partition_regular.jl")
    end
    # @safetestset "TreePartition" begin
    #     include("partition_tree.jl")
    # end
    @safetestset "Algorithms" begin
        include("algorithms.jl")
    end
end


text:
function Fgamma(m::Int64,x::Float64,SMALL::Float64=1e-12)
    #println("Fgamma($m,$x)")
    x = max(x,SMALL) # Evidently needs underflow protection
    return 0.5*x^(-m-0.5)*gammainc(m+0.5,x)
end

function gammainc(a::Float64,x::Float64)
    # This is the series version of gamma from pyquante. For reasons I don't get, it
    # doesn't work around a=1. This works alright, but is only a stopgap solution
    # until Julia gets an incomplete gamma function programmed
    if abs(a-1) < 1e-3
        println("Warning: gammainc_series is known to have problems for a ~ 1")
    end
    if x < (a+1.0)
        #Use the series representation
        gam,gln = gser(a,x)
    else
        #Use continued fractions
        gamc,gln = gcf(a,x)
        gam = 1-gamc
    end
    return exp(gln)*gam
end

function gser(a::Float64,x::Float64,ITMAX::Int64=100,EPS::Float64=3e-9)
    # Series representation of Gamma. NumRec sect 6.1.
    gln=lgamma(a)
    if x == 0
        return 0,gln
    end
    ap = a
    delt = s = 1/a
    for i in 1:ITMAX
        ap += 1
        delt *= (x/ap)
        s += delt
        if abs(delt) < abs(s)*EPS
            break
        end
    end
    return s*exp(-x+a*log(x)-gln),gln
end

function gcf(a::Float64,x::Float64,ITMAX::Int64=200,EPS::Float64=3e-9,FPMIN::Float64=1e-30)
    #Continued fraction representation of Gamma. NumRec sect 6.1"
    gln=lgamma(a)
    b=x+1.0-a
    c=1.0/FPMIN
    d=1.0/b
    h=d
    for i in 1:ITMAX
        an=-i*(i-a)
        b=b+2.0
        d=an*d+b
        if abs(d) < FPMIN
            d=FPMIN
        end
        c=b+an/c
        if abs(c) < FPMIN
            c=FPMIN
        end
        d=1.0/d
        delt=d*c
        h=h*delt
        if abs(delt-1.0) < EPS
            break
        end
    end
    gammcf = exp(-x+a*log(x)-gln)*h
    return gammcf,gln
end

function test_gamma()
    # gammainc test functions. Test values taken from Mathematica
    @testset "test_gamma" begin

    @test maximum([gammainc(0.5,float(x)) for x in 0:10]
            -[0, 1.49365, 1.69181, 1.7471, 1.76416, 1.76968,
                1.77151, 1.77213, 1.77234, 1.77241, 1.77244]) < 1e-5

    @test maximum([gammainc(1.5,float(x)) for x in 0:10]
            -[0, 1.49365, 1.69181, 1.7471, 1.76416, 1.76968,
                1.77151, 1.77213, 1.77234, 1.77241, 1.77244]) < 1e-5
    @test maximum([gammainc(2.5,float(x)) for x in 0:10]
            -[0, 0.200538, 0.59898, 0.922271, 1.12165, 1.22933,
                1.2831, 1.30859, 1.32024, 1.32542, 1.32768]) < 1e-5
    end #testset
end

function test_fgamma()
    @testset "test_fgamma" begin
    for (x,res) in [(0.,1),
                    (30.,0.161802159),
                    (60.,0.114411404),
                    (90.,0.0934165203),
                    (120.,0.08090108),
                    (300.,0.051166336)]
        @test res ≈ Fgamma(0,x)
    end
    end # testset
end


text:
## ICD10CM #####
"""
    ICD10CM

An ICD-10-CM code

Fields are `level`, indicating the number of digits in the code, and `data` containing the code.
"""
struct ICD10CM <: AbstractICD10
  level::Int8
  data::NTuple{7,UInt8}
end

"""
    ICD10CM(str::String, validateinput = false)

Create an ICD10CM code from a string.

ICD10CM codes have the form `ANC[.][C[C[C[C]]]]` where `A` is a letter A-Z,
`N` is a digit 0-9, and `C` is a letter or digit. Parts in brackets are optional.

If `validateinput=true` then the input string is checked for valid format using
a regex.
"""
function ICD10CM(str::AbstractString, validateinput = false)
  punct = occursin(".", str) ? true : false
  if validateinput
    validICD10CMinput(str) || throw(
      DomainError(
        str,
        "ICD-10-CM codes should have format `ANC[.][C[C[C[C]]]]` where `A` is a letter A-Z, `N` is a digit 0-9, and `C` is a letter or digit, and parts in brackets are optional",
      ),
    )
  end
  level = punct ? length(str) - 1 : length(str)
  ch4 = punct ? 5 : 4

  if level == 3
    return ICD10CM(
      level,
      (NTuple{3,UInt8}.(str)..., UInt8(0), UInt8(0), UInt8(0), UInt8(0)),
    )
  elseif level == 4
    return ICD10CM(
      level,
      (NTuple{3,UInt8}.(str)..., UInt8(str[ch4]), UInt8(0), UInt8(0), UInt8(0)),
    )
  elseif level == 5
    return ICD10CM(
      level,
      (NTuple{3,UInt8}.(str)..., UInt8(str[ch4]), UInt8(str[ch4+1]), UInt8(0), UInt8(0)),
    )
  elseif level == 6
    return ICD10CM(
      level,
      (
        NTuple{3,UInt8}.(str)...,
        UInt8(str[ch4]),
        UInt8(str[ch4+1]),
        UInt8(str[ch4+2]),
        UInt8(0),
      ),
    )
  elseif level == 7
    return ICD10CM(
      level,
      (
        NTuple{3,UInt8}.(str)...,
        UInt8(str[ch4]),
        UInt8(str[ch4+1]),
        UInt8(str[ch4+2]),
        UInt8(str[ch4+3]),
      ),
    )
  else
    throw(DomainError(str, "ICD-10-CM codes should be 3-7 characters excluding the period"))
  end
end

"""
    ICD10CM(icd::T) where T <: AbstractICD10

Create an ICD-10-CM code from another type of ICD-10 code.

Note this just changes the type. It does not translate concepts between versions.

Most ICD-10 versions have up to 5 digits, while ICD-10-CM has up to 7, so this conversion cannot be exact.
"""
ICD10CM(icd::T) where {T<:AbstractICD10} = ICD10(icd)

"""
    ICD10(icd::ICD10CM)

Convert from ICD10CM to ICD10.

Note that ICD-10-CM codes may have up to 7 digits, while other ICD-10 versions
only have up to 5.

**This function truncates longer codes.**

Note this just changes the type. It does not translate concepts between versions.
"""
ICD10(icd::ICD10CM) = ICD10(min(icd.level, 5), icd.data[1:5])

"""
    ICD10AM(icd::ICD10CM)

Convert from ICD10CM to ICD10AM.

Note that ICD-10-CM codes may have up to 7 digits, while other ICD-10 versions
only have up to 5.

**This function truncates longer codes.**

Note this just changes the type. It does not translate concepts between versions.
"""
ICD10AM(icd::ICD10CM) = ICD10AM(min(icd.level, 5), icd.data[1:5])

"""
    ICD10CA(icd::ICD10CM)

Convert from ICD10CM to ICD10CA.

Note that ICD-10-CM codes may have up to 7 digits, while other ICD-10 versions
only have up to 5.

**This function truncates longer codes.**

Note this just changes the type. It does not translate concepts between versions.
"""
ICD10CA(icd::ICD10CM) = ICD10CA(min(icd.level, 5), icd.data[1:5])

"""
    ICD10GM(icd::ICD10CM)

Convert from ICD10CM to ICD10GM.

Note that ICD-10-CM codes may have up to 7 digits, while other ICD-10 versions
only have up to 5.

**This function truncates longer codes.**

Note this just changes the type. It does not translate concepts between versions.
"""
ICD10GM(icd::ICD10CM) = ICD10GM(min(icd.level, 5), icd.data[1:5])

function validICD10CMinput(str)
  icdfmt = r"^[[:upper:]][[:digit:]][[:upper:][:digit:]](\.){0,1}([[:upper:][:digit:]]{1,4}){0,1}$"
  return occursin(icdfmt, str)
end


text:
"""
    FBMKernel(; h::Real=0.5)

Fractional Brownian motion kernel with Hurst index `h` from (0,1) given by
```
    κ(x,y) =  ( |x|²ʰ + |y|²ʰ - |x-y|²ʰ ) / 2
```

For `h=1/2`, this is the Wiener Kernel, for `h>1/2`, the increments are
positively correlated and for `h<1/2` the increments are negatively correlated.
"""
struct FBMKernel{T<:Real} <: Kernel
    h::Vector{T}
    function FBMKernel(; h::T=0.5) where {T<:Real}
        @assert 0.0 <= h <= 1.0 "FBMKernel: Given Hurst index h is invalid."
        return new{T}([h])
    end
end

@functor FBMKernel

function (κ::FBMKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    modX = sum(abs2, x)
    modY = sum(abs2, y)
    modXY = sqeuclidean(x, y)
    h = first(κ.h)
    return (modX^h + modY^h - modXY^h)/2
end

(κ::FBMKernel)(x::Real, y::Real) = (abs2(x)^first(κ.h) + abs2(y)^first(κ.h) - abs2(x - y)^first(κ.h)) / 2

Base.show(io::IO, κ::FBMKernel) = print(io, "Fractional Brownian Motion Kernel (h = ", first(κ.h), ")")

_fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h)/2

_mod(x::AbstractVector{<:Real}) = abs2.(x)
_mod(x::ColVecs) = vec(sum(abs2, x.X; dims=1))
_mod(x::RowVecs) = vec(sum(abs2, x.X; dims=2))

function kernelmatrix(κ::FBMKernel, x::AbstractVector)
    modx = _mod(x)
    modxx = pairwise(SqEuclidean(), x)
    return _fbm.(modx, modx', modxx, κ.h)
end

function kernelmatrix!(K::AbstractMatrix, κ::FBMKernel, x::AbstractVector)
    modx = _mod(x)
    pairwise!(K, SqEuclidean(), x)
    K .= _fbm.(modx, modx', K, κ.h)
    return K
end

function kernelmatrix(κ::FBMKernel, x::AbstractVector, y::AbstractVector)
    modxy = pairwise(SqEuclidean(), x, y)
    return _fbm.(_mod(x), _mod(y)', modxy, κ.h)
end

function kernelmatrix!(
    K::AbstractMatrix,
    κ::FBMKernel,
    x::AbstractVector,
    y::AbstractVector,
)
    pairwise!(K, SqEuclidean(), x, y)
    K .= _fbm.(_mod(x), _mod(y)', K, κ.h)
    return K
end


text:
using Bokeh
using Base.Test

# check js and css files are there (in case they've been downloaded)
bkfiles = readdir(joinpath(dirname(@__FILE__), "..", "deps", "bokehjs"))
@test length(filter(x -> endswith(x, ".js"), bkfiles))  == 2
@test length(filter(x -> endswith(x, ".css"), bkfiles)) == 2

# try setting all possible params
# some lime autoopen actually need to be set for travis
Bokeh.debug(true)
Bokeh.debug(false)
w = width()
width(w)
h = height()
height(h)
dft_glyphs = glyphs()
glyphs(dft_glyphs)
dft_size = glyphsize()
glyphsize(dft_size)
t = title()
title(t)
Bokeh.noshow(true)
Bokeh.filewarnings(false)

# check they're correct:
# Bokeh.parse_spec("r--") |> println

# Bokeh.parse_spec("r:") |> println

exdir = Pkg.dir("Bokeh", "test", "test-cases")
testingdir = "/tmp/bokeh_js_testing"
mkpath(testingdir)
cd(testingdir)
for ex in filter(f -> endswith(f, ".jl"), readdir(exdir))
	evalfile(joinpath(exdir, ex))
	hold(false)
end


text:
struct FakeNaturalDiffWithIsApprox  # For testing overloading isapprox(::Composite) works:
    x
end
function Base.isapprox(c::Composite, d::FakeNaturalDiffWithIsApprox; kwargs...)
    return isapprox(c.x, d.x, kwargs...)
end
function Base.isapprox(d::FakeNaturalDiffWithIsApprox, c::Composite; kwargs...)
    return isapprox(c.x, d.x, kwargs...)
end

@testset "check_result.jl" begin
    @testset "_check_add!!_behavour" begin
        check = ChainRulesTestUtils._check_add!!_behavour

        check(10.0, 2.0)
        check(11.0, Zero())
        check([10.0, 20.0],  @thunk([2.0, 0.0]))

        check(12.0, InplaceableThunk(@thunk(2.0), X̄ -> error("Should not have in-placed")))

        check([10.0, 20.0], InplaceableThunk(
            @thunk([2.0, 0.0]),
            X̄ -> (X̄[1] += 2.0; X̄)
        ))
        @test fails(()->check([10.0, 20.0], InplaceableThunk(
            @thunk([2.0, 0.0]),
            X̄ -> (X̄[1] += 3.0; X̄),
        )))
    end

    @testset "check_equal" begin

        @testset "possive cases" begin
            check_equal(1.0, 1.0)
            check_equal(1.0 + im, 1.0 + im)
            check_equal(1.0, 1.0+1e-10)  # isapprox _behaviour
            check_equal((1.5, 2.5, 3.5), (1.5, 2.5, 3.5 + 1e-10))

            check_equal(Zero(), 0.0)

            check_equal([1.0, 2.0], [1.0, 2.0])
            check_equal([[1.0], [2.0]], [[1.0], [2.0]])

            check_equal(@thunk(10*0.1*[[1.0], [2.0]]), [[1.0], [2.0]])

            check_equal(
                Composite{Tuple{Float64, Float64}}(1.0, 2.0),
                Composite{Tuple{Float64, Float64}}(1.0, 2.0)
            )

            diag_eg = Diagonal(randn(5))
            check_equal( # Structual == Structural
                Composite{typeof(diag_eg)}(diag=diag_eg.diag),
                Composite{typeof(diag_eg)}(diag=diag_eg.diag)
            )
            check_equal( # Structural == Natural
                Composite{typeof(diag_eg)}(diag=diag_eg.diag),
                diag_eg
            )

            T = (a=1.0, b=2.0)
            check_equal(
                Composite{typeof(T)}(a=1.0),
                Composite{typeof(T)}(a=1.0, b=Zero())
            )
            check_equal(
                Composite{typeof(T)}(a=1.0),
                Composite{typeof(T)}(a=1.0+1e-10, b=Zero())
            )

            check_equal(
                Composite{FakeNaturalDiffWithIsApprox}(; x=1.4),
                FakeNaturalDiffWithIsApprox(1.4)
            )
            check_equal(
                FakeNaturalDiffWithIsApprox(1.4),
                Composite{FakeNaturalDiffWithIsApprox}(; x=1.4)
            )
        end
        @testset "negative case" begin
            @test fails(()->check_equal(1.0, 2.0))
            @test fails(()->check_equal(1.0 + im, 1.0 - im))
            @test fails(()->check_equal((1.5, 2.5, 3.5), (1.5, 2.5, 4.5)))

            @test fails(()->check_equal(Zero(), 20.0))
            @test fails(()->check_equal(10.0, Zero()))

            @test fails(()->check_equal([1.0, 2.0], [1.0, 3.9]))
            @test fails(()->check_equal([[1.0], [2.0]], [[1.1], [2.0]]))

            @test fails(()->check_equal(@thunk(10*[[1.0], [2.0]]), [[1.0], [2.0]]))
        end
        @testset "type negative" begin
            @test fails() do  # these have different primals so should not be equal
                check_equal(
                    Composite{Tuple{Float32, Float32}}(1f0, 2f0),
                    Composite{Tuple{Float64, Float64}}(1.0, 2.0)
                )
            end
            @test fails() do
                check_equal((1.0, 2.0), Composite{Tuple{Float64, Float64}}(1.0, 2.0))
            end
        end

        @testset "TestIterator" begin
            data = randn(3)
            iter1 = TestIterator(data, Base.HasLength(), Base.HasEltype())
            iter2 = TestIterator(data, Base.HasLength(), Base.EltypeUnknown())
            check_equal(iter2, iter1)

            iter3 = TestIterator(data .+ 1e-10, Base.HasLength(), Base.HasEltype())
            check_equal(iter3, iter1)

            iter_bad = TestIterator(data .+ 010, Base.HasLength(), Base.HasEltype())
            @test fails(()->check_equal(iter_bad, iter1))
        end
    end
end