text:
#=

EGeoHiSSE equations

Ignacio Quintero Mächler

t(-_-t)

Created 18 03 2019
=#

"""
    exp_expr(k    ::Int64,
             h    ::Int64,
             ny   ::Int64,
             model::NTuple{3,Bool})

Return exponential expression for one time evaluation of covariates.
"""
function exp_expr(k    ::Int64,
                  h    ::Int64,
                  ny   ::Int64,
                  model::NTuple{3,Bool})

  # last non β parameters
  bbase = h*(k+1) + 2*k*h + k*(k-1)*h + h*(h-1)

  # ncov
  ncov = model[1]*k + model[2]*k + model[3]*k*(k-1)

  # y per parameter
  yppar = isone(ny) ? 1 : div(ny,ncov)

  # starting indices for models 2 and 3
  m2s = model[1]*k*h
  m3s = m2s + model[2]*k*h

  y2s = model[1]*yppar*k
  y3s = y2s + model[2]*yppar*k

  # start expression
  ex = quote end
  pop!(ex.args)

  isone(ny) && push!(ex.args, :(r1 = r[1]))

  # speciation
  if model[1]
    for j = Base.OneTo(h), i = Base.OneTo(k)
      coex = Expr(:call, :+)
      for yi in Base.OneTo(yppar)
        rex = isone(ny) ? :(r1) : :(r[$(yi+yppar*(i-1))])
        push!(coex.args, :(p[$(bbase + yi + yppar*((i-1) + k*(j-1)))] * $rex))
      end
      push!(ex.args, :(eaft[$(i + k*(j-1))] = p[$(i + (k+1)*(j-1))]*exp($coex)))
    end
  end
  # extinction
  if model[2]
    for j = Base.OneTo(h), i = Base.OneTo(k)
      coex = Expr(:call, :+)
      for yi in Base.OneTo(yppar)
        rex = isone(ny) ? :(r1) : :(r[$(y2s + yi+yppar*(i-1))])
        push!(coex.args,
          :(p[$(bbase + yi + m2s*yppar + yppar*((i-1) + k*(k-1)*(j-1)))] * $rex))
      end
      push!(ex.args, :(eaft[$(i + k*(j-1) + m2s)] =
        p[$((k+1)*h + i + k*(j-1))]*exp($coex)))
    end
  end
  # transition
  if model[3]
    for j = Base.OneTo(h), i = Base.OneTo(k*(k-1))
      coex = Expr(:call, :+)
      for yi in Base.OneTo(yppar)
        rex = isone(ny) ? :(r1) : :(r[$(y3s + yi+yppar*(i-1))])
        push!(coex.args,
          :(p[$(bbase + yi + m3s*yppar + yppar*((i-1) + k*(k-1)*(j-1)))] * $rex))
      end
      push!(ex.args, :(eaft[$(i + k*(j-1) + m3s)] =
        p[$(h*(2k+1) + i + k*(k-1)*(j-1))]*exp($coex)))
    end
  end

  return ex
end

"""
    noevents_expr(si   ::Int64,
                  s    ::Sgh,
                  ls   ::Int64,
                  oa   ::Array{Int64,1},
                  k    ::Int64,
                  h    ::Int64,
                  ns   ::Int64,
                  model::NTuple{3,Bool},
                  ext  ::Bool)

Return expression for no events in Vector form.
"""
function noevents_expr(si   ::Int64,
                       s    ::Sgh,
                       ls   ::Int64,
                       oa   ::Array{Int64,1},
                       k    ::Int64,
                       h    ::Int64,
                       ns   ::Int64,
                       model::NTuple{3,Bool},
                       ext  ::Bool)

  # starting indices for models 2 and 3
  m2s = model[1]*k*h
  m3s = m2s + model[2]*k*h

  # if extinction ODE
  wu = ext ? ns : 0

  ts = isone(ls) ? 0 : k*h*k

  # between-region speciation
  ex = Expr(:call, :+)
  ls > 1 && push!(ex.args, :($(2.0^(ls-1) - 1.0) * p[$(k+1+s.h*(k+1))]))

  for v in s.g

    #speciation
    if model[1]
      push!(ex.args, :(eaft[$(v + s.h*k)]))
    else
      push!(ex.args, :(p[$(v + s.h*(k+1))]))
    end

    # extinction
    if model[2]
      push!(ex.args, :(eaft[$(m2s + v + s.h*k)]))
    else
      push!(ex.args, :(p[$(v + (k+1)*h + s.h*k + ts)]))
    end

    # dispersal
    if model[3]
      for j = oa
        j -= v <= j ? 1 : 0
        push!(ex.args,
          :(eaft[$(m3s + s.h*k*(k-1) + (k-1)*(v-1) + j)]))
      end
    else
      for j = oa
        j -= v <= j ? 1 : 0
        push!(ex.args,
          :(p[$(2h*k + h + s.h*k*(k-1) + (k-1)*(v-1) + j)]))
      end
    end
  end

  # add hidden state shifts
  for hi in setdiff(0:(h-1), s.h)
    hi -= s.h <= hi ? 0 : -1
    push!(ex.args, :(p[$(h*(3k + 1 + k*(k-1)) + s.h*(h-1) + hi)]))
  end

  # multiply by u
  ex = :(-1.0 * $ex * u[$(si + wu)])

  return ex
end

"""
    noevents_expr(si     ::Int64,
                  s      ::Sgh,
                  ls     ::Int64,
                  oa     ::Array{Int64,1},
                  k      ::Int64,
                  h      ::Int64,
                  ns     ::Int64,
                  i      ::Int64,
                  model  ::NTuple{3,Bool},
                  ext    ::Bool)

Return expression for no events in Matrix form.
"""
function noevents_expr(si     ::Int64,
                       s      ::Sgh,
                       ls     ::Int64,
                       oa     ::Array{Int64,1},
                       k      ::Int64,
                       h      ::Int64,
                       ns     ::Int64,
                       i      ::Int64,
                       model  ::NTuple{3,Bool})

  # starting indices for models 2 and 3
  m2s = model[1]*k*h
  m3s = m2s + model[2]*k*h

  ts = isone(ls) ? 0 : k*h*k

  # between-region speciation
  ex = Expr(:call, :+)
  ls > 1 && push!(ex.args, :($(2.0^(ls-1) - 1.0) * p[$(k+1+s.h*(k+1))]))

  for v in s.g

    #speciation
    if model[1]
      push!(ex.args, :(eaft[$(v + s.h*k)]))
    else
      push!(ex.args, :(p[$(v + s.h*(k+1))]))
    end

    # extinction
    if model[2]
      push!(ex.args, :(eaft[$(m2s + v + s.h*k)]))
    else
      push!(ex.args, :(p[$(v + (k+1)*h + s.h*k + ts)]))
    end

    # dispersal
    if model[3]
      for j = oa
        j -= v <= j ? 1 : 0
        push!(ex.args,
          :(eaft[$(m3s + s.h*k*(k-1) + (k-1)*(v-1) + j)]))
      end
    else
      for j = oa
        j -= v <= j ? 1 : 0
        push!(ex.args,
          :(p[$(2h*k + h + s.h*k*(k-1) + (k-1)*(v-1) + j)]))
      end
    end
  end

  # add hidden state shifts
  for hi in setdiff(0:(h-1), s.h)
    hi -= s.h <= hi ? 0 : -1
    push!(ex.args, :(p[$(h*(3k + 1 + k*(k-1)) + s.h*(h-1) + hi)]))
  end

  # multiply by u
  ex = :(-1.0 * $ex * u[$si, $i])

  return ex
end

"""
    localext_expr(s ::Sgh,
                  S ::Array{Sgh,1},
                  k ::Int64,
                  h ::Int64,
                  model::NTuple{3,Bool})

Return expression for local extinction in Vector form.
"""
function localext_expr(s ::Sgh,
                       S ::Array{Sgh,1},
                       k ::Int64,
                       h ::Int64,
                       model::NTuple{3,Bool})

  ex = Expr(:call, :+)

  if model[2]
    # starting indices for models 2 and 3
    m2s = model[1]*k*h

    for i = s.g
      push!(ex.args, :(eaft[$(m2s + i + s.h*k)] *
                       u[$(findfirst(x -> isSghequal(x,
                                          Sgh(setdiff(s.g, i),s.h)), S))]))
    end
  else
    for i = s.g
      push!(ex.args, :(p[$(i + (k+1)*h + s.h*k + k*h*k)] *
                       u[$(findfirst(x -> isSghequal(x,
                                          Sgh(setdiff(s.g, i),s.h)), S))]))
    end
  end

  return ex
end

"""
    localext_expr(s     ::Sgh,
                  S     ::Array{Sgh,1},
                  k     ::Int64,
                  h     ::Int64,
                  i     ::Int64
                  model ::NTuple{3,Bool})

Return expression for local extinction in Matrix form.
"""
function localext_expr(s     ::Sgh,
                       S     ::Array{Sgh,1},
                       k     ::Int64,
                       h     ::Int64,
                       i     ::Int64,
                       model ::NTuple{3,Bool})

  ex = Expr(:call, :+)

  if model[2]
    # starting indices for models 2 and 3
    m2s = model[1]*k*h

    for j = s.g
      push!(ex.args, :(eaft[$(m2s + j + s.h*k)] *
                       u[$(findfirst(x -> isSghequal(x,
                                          Sgh(setdiff(s.g, j),s.h)), S)),$i]))
    end
  else
    for j = s.g
      push!(ex.args, :(p[$(j + (k+1)*h + s.h*k + k*h*k)] *
                       u[$(findfirst(x -> isSghequal(x,
                                          Sgh(setdiff(s.g, j),s.h)), S)),$i]))
    end
  end

  return ex
end

"""
    dispersal_expr(s    ::Sgh,
                   oa   ::Array{Int64,1},
                   S    ::Array{Sgh,1},
                   ns   ::Int64,
                   k    ::Int64,
                   h    ::Int64,
                   model::NTuple{3,Bool},
                   ext  ::Bool)

Return expression for dispersal in Vector form.
"""
function dispersal_expr(s    ::Sgh,
                        oa   ::Array{Int64,1},
                        S    ::Array{Sgh,1},
                        ns   ::Int64,
                        k    ::Int64,
                        h    ::Int64,
                        model::NTuple{3,Bool},
                        ext  ::Bool)

  m3s = model[1]*k*h + model[2]*k*h

  wu = ext ? ns : 0

  ida = findall(x -> length(union(s.g, x.g))     == length(s.g)+1 &&
                     length(intersect(s.g, x.g)) == length(s.g)   &&
                     s.h == x.h, S)

  ex = Expr(:call, :+)

  # if dispersal covariate
  if model[3]
    for a = s.g, (i, j) = enumerate(oa)
      j -= a <= j ? 1 : 0
      push!(ex.args, :(eaft[$(m3s + s.h*k*(k-1) + (k-1)*(a-1) + j)] *
                       u[$(ida[i] + wu)]))
    end
  else
    for a = s.g, (i, j) = enumerate(oa)
      j -= a <= j ? 1 : 0
      push!(ex.args, :(p[$(2h*k + h + s.h*k*(k-1) + (k-1)*(a-1) + j)] *
                       u[$(ida[i] + wu)]))
    end
  end

  return ex
end

"""
    dispersal_expr(s    ::Sgh,
                   oa   ::Array{Int64,1},
                   S    ::Array{Sgh,1},
                   ns   ::Int64,
                   k    ::Int64,
                   h    ::Int64,
                   i    ::Int64,
                   model::NTuple{3,Bool})

Return expression for dispersal in Matrix form.
"""
function dispersal_expr(s    ::Sgh,
                        oa   ::Array{Int64,1},
                        S    ::Array{Sgh,1},
                        ns   ::Int64,
                        k    ::Int64,
                        h    ::Int64,
                        i    ::Int64,
                        model::NTuple{3,Bool})

  m3s = model[1]*k*h + model[2]*k*h

  ida = findall(x -> length(union(s.g, x.g))     == length(s.g)+1 &&
                     length(intersect(s.g, x.g)) == length(s.g)   &&
                     s.h == x.h, S)

  ex = Expr(:call, :+)

  # if dispersal covariate
  if model[3]
    for a = s.g, (l, j) = enumerate(oa)
      j -= a <= j ? 1 : 0
      push!(ex.args, :(eaft[$(m3s + s.h*k*(k-1) + (k-1)*(a-1) + j)] *
                       u[$(ida[l]), $i]))
    end
  else
    for a = s.g, (l, j) = enumerate(oa)
      j -= a <= j ? 1 : 0
      push!(ex.args, :(p[$(2h*k + h + s.h*k*(k-1) + (k-1)*(a-1) + j)] *
                       u[$(ida[l]), $i]))
    end
  end

  return ex
end

"""
    hidtran_expr(s  ::Sgh,
                 S  ::Array{Sgh,1},
                 ns ::Int64,
                 k  ::Int64,
                 h  ::Int64,
                 i  ::Int64)

Return expression for hidden states transitions in Matrix form.
"""
function hidtran_expr(s  ::Sgh,
                      S  ::Array{Sgh,1},
                      ns ::Int64,
                      k  ::Int64,
                      h  ::Int64,
                      i  ::Int64)

  hs = findall(x -> isequal(s.g, x.g) &&
                   !isequal(s.h, x.h), S)

  ex = Expr(:call, :+)
  for j in hs
    hi = S[j].h
    hi -= s.h <= hi ? 0 : -1
    push!(ex.args, :(p[$(h*(3k + 1 + k*(k-1)) + s.h*(h-1) + hi)] *
                     u[$j, $i]))
  end
  return ex
end

"""
    hidtran_expr(s  ::Sgh,
                 S  ::Array{Sgh,1},
                 ns ::Int64,
                 k  ::Int64,
                 h  ::Int64,
                 ext::Bool)

Return expression for hidden states transitions in Vector form.
"""
function hidtran_expr(s  ::Sgh,
                      S  ::Array{Sgh,1},
                      ns ::Int64,
                      k  ::Int64,
                      h  ::Int64,
                      ext::Bool)

  wu = ext ? ns : 0

  hs = findall(x -> isequal(s.g, x.g) &&
                   !isequal(s.h, x.h), S)

  ex = Expr(:call, :+)
  for i in hs
    hi = S[i].h
    hi -= s.h <= hi ? 0 : -1
    push!(ex.args, :(p[$(h*(3k + 1 + k*(k-1)) + s.h*(h-1) + hi)] *
                     u[$(i + wu)]))
  end
  return ex
end

"""
    wrspec_expr(si   ::Int64,
                s    ::Sgh,
                ns   ::Int64,
                k    ::Int64,
                i    ::Int64,
                model::NTuple{3,Bool})

Return expression for within-region speciation in Matrix form.
"""
function wrspec_expr(si   ::Int64,
                     s    ::Sgh,
                     ns   ::Int64,
                     k    ::Int64,
                     i    ::Int64,
                     model::NTuple{3,Bool})

  # if model speciation
  if model[1]
    if isone(length(s.g))
      ex = Expr(:call, :*, 2.0)
      for j = s.g
        push!(ex.args, :(eaft[$(j + s.h*k)] *
               rE[$(j + s.h*(2^k-1))] * u[$(j + s.h*(2^k-1)), $i]))
      end
    else
      ex = Expr(:call, :+)
      for j = s.g
        push!(ex.args, :(eaft[$(j + s.h*k)] *
          (rE[$(j + s.h*(2^k-1))] * u[$(si), $i] +
           rE[$(si)] * u[$(j + s.h*(2^k-1)), $i])))
      end
    end

  # if model *NOT* speciation
  else
    if isone(length(s.g))
      ex = :(2.0 * p[$si] * rE[$(si)] * u[$(si)])
    else
      ex = Expr(:call, :+)
      for j = s.g
        push!(ex.args, :(p[$(j + s.h*(k+1))] *
          (rE[$(j + s.h*(2^k-1))] * u[$(si), $i] +
           rE[$(si)] * u[$(j + s.h*(2^k-1)), $i])))
      end
    end
  end

  return ex
end

"""
    wrspec_expr(si   ::Int64,
                s    ::Sgh,
                ns   ::Int64,
                k    ::Int64,
                model::NTuple{3,Bool})

Return expression for within-region speciation in Vector form.
"""
function wrspec_expr(si   ::Int64,
                     s    ::Sgh,
                     ns   ::Int64,
                     k    ::Int64,
                     model::NTuple{3,Bool})

  # if model speciation
  if model[1]
    if isone(length(s.g))
      ex = Expr(:call, :*, 2.0)
      for i = s.g
        push!(ex.args, :(eaft[$(i + s.h*k)] *
               u[$(i + s.h*(2^k-1) + ns)] * u[$(i + s.h*(2^k-1))]))
      end
    else
      ex = Expr(:call, :+)
      for i = s.g
        push!(ex.args, :(eaft[$(i + s.h*k)] *
          (u[$(i + s.h*(2^k-1) + ns)] * u[$(si)] +
           u[$(si + ns)] * u[$(i + s.h*(2^k-1))])))
      end
    end

  # if model *NOT* speciation
  else
    if isone(length(s.g))
      ex = :(2.0 * p[$si] * u[$(si + ns)] * u[$(si)])
    else
      ex = Expr(:call, :+)
      for i = s.g
        push!(ex.args, :(p[$(i + s.h*(k+1))] *
          (u[$(i + s.h*(2^k-1) + ns)] * u[$(si)] +
           u[$(si + ns)] * u[$(i + s.h*(2^k-1))])))
      end
    end
  end

  return ex
end

"""
    brspec_expr(s   ::Sgh,
                S   ::Array{Sgh,1},
                ns  ::Int64,
                k   ::Int64,
                i   ::Int64,
                ext ::Bool,
                flow::Bool)

Return expression for between-region speciation in Vector form.
"""
function brspec_expr(s   ::Sgh,
                     S   ::Array{Sgh,1},
                     ns  ::Int64,
                     k   ::Int64,
                     ext ::Bool,
                     flow::Bool)

  if ext
    va = vicsubsets(s.g)[1:div(end,2)]
    wu1 = flow ? 0 : ns
    wu2 = flow ? 0 : ns
  else
    va = vicsubsets(s.g)
    wu1 = flow ? 0 : ns
    wu2 = 0
  end

  ex = Expr(:call, :+)
  for (la, ra) = va
  push!(ex.args,
    :(u[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h + wu1)] *
      u[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h + wu2)]))
  end

  ex = :(p[$(k+1+s.h*(k+1))] * $ex)

  return ex
end

"""
    brspec_expr(s   ::Sgh,
                S   ::Array{Sgh,1},
                ns  ::Int64,
                k   ::Int64,
                i   ::Int64)

Return expression for between-region speciation in Matrix form.
"""
function brspec_expr(s   ::Sgh,
                     S   ::Array{Sgh,1},
                     ns  ::Int64,
                     k   ::Int64,
                     i   ::Int64)

  va = vicsubsets(s.g)

  ex = Expr(:call, :+)

  for (la, ra) = va
  push!(ex.args,
    :(rE[$(findfirst(x -> isequal(ra,x.g), S) + (2^k-1)*s.h)] *
       u[$(findfirst(x -> isequal(la,x.g), S) + (2^k-1)*s.h), $i]))
  end

  ex = :(p[$(k+1+s.h*(k+1))] * $ex)

  return ex
end

"""
    ext_expr(s    ::Sgh,
             S    ::Array{Sgh,1},
             ns   ::Int64,
             k    ::Int64,
             h    ::Int64,
             model::NTuple{3,Bool})

Return expression for extinction.
"""
function ext_expr(s    ::Sgh,
                  S    ::Array{Sgh,1},
                  ns   ::Int64,
                  k    ::Int64,
                  h    ::Int64,
                  model::NTuple{3,Bool},
                  flow ::Bool)

  wu = flow ? 0 : ns

  m2s = model[1]*k*h

  if isone(length(s.g))
    if model[2]
      ex = Expr(:call, :+)
      for i in s.g push!(ex.args, :(eaft[$(m2s + k*s.h + i)])) end
    else
      ex = Expr(:call, :+)
      for i in s.g push!(ex.args, :(p[$((k+1)*h + k*s.h + i)])) end
    end
  else
    if model[2]
      ex = Expr(:call, :+)
      for i = s.g
        push!(ex.args, :(eaft[$(m2s + k*s.h + i)] *
                         u[$(findfirst(x -> isSghequal(x,
                                            Sgh(setdiff(s.g, i),s.h)), S) + wu)]))
      end
    else
      ex = Expr(:call, :+)
      for i = s.g
        push!(ex.args, :(p[$(i + (k+1)*h + s.h*k + k*h*k)] *
                         u[$(findfirst(x -> isSghequal(x,
                                            Sgh(setdiff(s.g, i),s.h)), S) + wu)]))
      end
    end
  end

  return ex
end

"""
    wrsext_expr(si   ::Int64,
                s    ::Sgh,
                ns   ::Int64,
                k    ::Int64,
                model::NTuple{3,Bool},
                flow ::Bool)

Return expression of extinction for within-region speciation.
"""
function wrsext_expr(si   ::Int64,
                     s    ::Sgh,
                     ns   ::Int64,
                     k    ::Int64,
                     model::NTuple{3,Bool},
                     flow ::Bool)

  wu = flow ? 0 : ns

  # if speciation model
  if model[1]
    if isone(length(s.g))
      for i = s.g
        ex = :(eaft[$(i + s.h*k)] * u[$(i + s.h*(2^k-1) + wu)]^2)
      end
    else
      ex = Expr(:call, :+)
      for i = s.g
        push!(ex.args, :(eaft[$(i + s.h*k)] *
                         u[$(i + s.h*(2^k-1) + wu)] * u[$(si + wu)]))
      end
    end
  # if *NOT* speciation model
  else
    if isone(length(s.g))
      ex = :(p[$si] * u[$(si + wu)]^2)
    else
      ex = Expr(:call, :+)
      for i = s.g
        push!(ex.args, :(p[$(i + s.h*(k+1))] *
                         u[$(i + s.h*(2^k-1) + wu)] * u[$(si + wu)]))
      end
    end
  end

  return ex
end

"""
    geohisse_full(du  ::Array{Float64,1},
                  u   ::Array{Float64,1},
                  p   ::Array{Float64,1},
                  t   ::Float64,
                  r   ::Array{Float64,1},
                  eaft::Array{Float64,1},
                  af! ::Function,
                  ::Val{k},
                  ::Val{h},
                  ::Val{ny},
                  ::Val{model}) where {k, h, ny, model}

Creates Covariate GeoHiSSE ODE equation function for `k` areas,
`h` hidden states and `ny` covariates.
"""
@generated function geohisse_full(du  ::Array{Float64,1},
                                  u   ::Array{Float64,1},
                                  p   ::Array{Float64,1},
                                  t   ::Float64,
                                  r   ::Array{Float64,1},
                                  eaft::Array{Float64,1},
                                  af! ::Function,
                                  ::Val{k},
                                  ::Val{h},
                                  ::Val{ny},
                                  ::Val{model}) where {k, h, ny, model}

  # n states
  ns = (2^k - 1)*h

  # create individual areas subsets
  S = create_states(k, h)

  # start Expression for ODE equations
  eqs = quote end
  popfirst!(eqs.args)

  # add environmental function
  push!(eqs.args, :(af!(t, r)))

  # compute exponential for each β*covariable
  expex = exp_expr(k, h, ny, model)

  for ex in expex.args
    push!(eqs.args, ex)
  end

  for si = Base.OneTo(ns)

    # state range
    s = S[si]

    # length of geographic range
    ls = length(s.g)

    # which single areas do not occur in r
    oa = setdiff(1:k, s.g)

    #=
    likelihoods
    =#

    # no events
    nev = noevents_expr(si, s, ls, oa, k, h, ns, model, false)

    # remove if !isone(lr)
    if ls > 1
      # local extinction
      lex = localext_expr(s, S, k, h, model)

      # between-region speciation
      brs = brspec_expr(s, S, ns, k, false, false)
    end

    # dispersal
    # remove if lr == k
    if ls != k
      dis = dispersal_expr(s, oa, S, ns, k, h, model, false)
    end

    # hidden states transitions
    hid = h > 1 ? hidtran_expr(s, S, ns,k, h, false) : :0.0

    # within-region speciation
    wrs = wrspec_expr(si, s, ns, k, model)

    # push `D` equation to to eqs
    if isone(ls)
      push!(eqs.args,
        :(du[$si] = $nev + $dis + $hid + $wrs))
    elseif ls == k
      push!(eqs.args,
        :(du[$si] = $nev + $lex + $hid + $wrs + $brs))
    else
      push!(eqs.args,
        :(du[$si] = $nev + $lex + $dis + $hid + $wrs + $brs))
    end

    #=
    extinctions
    =#

    # no events
    nev = noevents_expr(si, s, ls, oa, k, h, ns, model, true)

    # extinction
    ext = ext_expr(s, S, ns, k, h, model, false)

    # dispersal and extinction
    if ls != k
      dis = dispersal_expr(s, oa, S, ns, k, h, model, true)
    end

    # hidden states transitions
    hid = h > 1 ? hidtran_expr(s, S, ns,k, h, true) : :0.0

    # within-region extinction
    wrs = wrsext_expr(si, s, ns, k, model, false)

    # between-region extinction
    if ls > 1
      brs = brspec_expr(s, S, ns, k, true, false)
    end

    # push `E` equation to to eqs
    if isone(ls)
      push!(eqs.args,
        :(du[$(si + ns)] = $nev + $ext + $dis + $hid + $wrs))
    elseif ls == k
      push!(eqs.args,
        :(du[$(si + ns)] = $nev + $ext + $hid + $wrs + $brs))
    else
      push!(eqs.args,
        :(du[$(si + ns)] = $nev + $ext + $dis + $hid + $wrs + $brs))
    end

  end

  ## aesthetic touches
  # sort `du`s
  eqs.args[[(end-(ns*2)+1):end...]] =
    eqs.args[append!([(end-(ns*2)+1):2:end...],
            [(end-(ns*2)+2):2:end...])]

  #Core.println(eqs)

  return quote
    @inbounds begin
      $eqs
    end
    return nothing
  end

end

"""
    make_egeohisse(::Val{k},
                   ::Val{h},
                   ::Val{ny},
                   ::Val{model},
                   af!::Function) where {k, h, ny, model}

Make closure for EGeoHiSSE.
"""
function make_egeohisse(::Val{k},
                        ::Val{h},
                        ::Val{ny},
                        ::Val{model},
                        af!::Function) where {k, h, ny, model}

  r    = Array{Float64,1}(undef, ny)
  eaft = Array{Float64,1}(undef,
    model[1]*k*h + model[2]*k*h + model[3]*k*(k-1)*h)

  # make ode function with closure
  ode_fun = (du::Array{Float64,1},
             u::Array{Float64,1},
             p::Array{Float64,1},
             t::Float64) ->
    begin
      geohisse_full(du::Array{Float64,1},
                    u ::Array{Float64,1},
                    p ::Array{Float64,1},
                    t ::Float64,
                    r,eaft,af!,Val(k),Val(h),Val(ny),Val(model))
      return nothing
    end

  return ode_fun
end

"""
    geohisse_E(du  ::Array{Float64,1},
               u   ::Array{Float64,1},
               p   ::Array{Float64,1},
               t   ::Float64,
               r   ::Array{Float64,1},
               eaft::Array{Float64,1},
               af! ::Function,
               ::Val{k},
               ::Val{h},
               ::Val{ny},
               ::Val{model}) where {k, h, ny, model}

Creates Covariate GeoHiSSE Extinction ODE equation for `k` areas,
`h` hidden states and `ny` covariates for the specific `model`.
"""
@generated function geohisse_E(du  ::Array{Float64,1},
                               u   ::Array{Float64,1},
                               p   ::Array{Float64,1},
                               t   ::Float64,
                               r   ::Array{Float64,1},
                               eaft::Array{Float64,1},
                               af! ::Function,
                               ::Val{k},
                               ::Val{h},
                               ::Val{ny},
                               ::Val{model}) where {k, h, ny, model}

  # n states
  ns = (2^k - 1)*h

  # create individual areas subsets
  S = create_states(k, h)

  # start Expression for ODE equations
  eqs = quote end
  popfirst!(eqs.args)

  # add environmental function
  push!(eqs.args, :(af!(t, r)))

  # compute exponential for each β*covariable
  expex = exp_expr(k, h, ny, model)

  for ex in expex.args
    push!(eqs.args, ex)
  end

  for si = Base.OneTo(ns)

    # state range
    s = S[si]

    # length of geographic range
    ls = length(s.g)

    # which single areas do not occur in r
    oa = setdiff(1:k, s.g)

    # no events
    nev = noevents_expr(si, s, ls, oa, k, h, ns, model, false)

    # extinction
    ext = ext_expr(s, S, ns, k, h, model, true)

    # dispersal and extinction
    if ls != k
      dis = dispersal_expr(s, oa, S, ns, k, h, model, false)
    end

    # hidden states transitions
    hid = h > 1 ? hidtran_expr(s, S, ns, k, h, false) : :0.0

    # within-region extinction
    wrs = wrsext_expr(si, s, ns, k, model, true)

    # between-region extinction
    if ls > 1
      brs = brspec_expr(s, S, ns, k, true, true)
    end

    # push `E` equation to to eqs
    if isone(ls)
      push!(eqs.args,
        :(du[$(si)] = $nev + $ext + $dis + $hid + $wrs))
    elseif ls == k
      push!(eqs.args,
        :(du[$(si)] = $nev + $ext + $hid + $wrs + $brs))
    else
      push!(eqs.args,
        :(du[$(si)] = $nev + $ext + $dis + $hid + $wrs + $brs))
    end

  end

  # print for checking
  #Core.println(eqs)

  return quote
    @inbounds begin
      $eqs
    end
    return nothing
  end

end

"""
    make_egeohisse_E(::Val{k},
                     ::Val{h},
                     ::Val{ny},
                     ::Val{model},
                     af!::Function) where {k, h, ny, model}

Make closure for Extinction function of EGeoHiSSE.
"""
function make_egeohisse_E(::Val{k},
                        ::Val{h},
                        ::Val{ny},
                        ::Val{model},
                        af!::Function) where {k, h, ny, model}

  r    = Array{Float64,1}(undef, ny)
  eaft = Array{Float64,1}(undef,
    model[1]*k*h + model[2]*k*h + model[3]*k*(k-1)*h)

  # make ode function with closure
  ode_fun = (du::Array{Float64,1},
             u::Array{Float64,1},
             p::Array{Float64,1},
             t::Float64) ->
    begin
      geohisse_E(du::Array{Float64,1},
                 u ::Array{Float64,1},
                 p ::Array{Float64,1},
                 t ::Float64,
                 r, eaft, af!, Val(k), Val(h), Val(ny), Val(model))
      return nothing
    end

  return ode_fun
end

"""
    geohisse_M(du  ::Array{Float64,1},
               u   ::Array{Float64,1},
               p   ::Array{Float64,1},
               t   ::Float64,
               r   ::Array{Float64,1},
               eaft::Array{Float64,1},
               af! ::Function,
               ::Val{k},
               ::Val{h},
               ::Val{ny},
               ::Val{model}) where {k, h, ny, model}

Creates Covariate GeoHiSSE Likelihood ODE equation function for `k` areas,
`h` hidden states and `ny` covariates in Matrix form.
"""
@generated function geohisse_M(du  ::Array{Float64,2},
                               u   ::Array{Float64,2},
                               pp  ::Array{Array{Float64,1},1},
                               t   ::Float64,
                               r   ::Array{Float64,1},
                               eaft::Array{Float64,1},
                               rE  ::Array{Float64,1},
                               af! ::Function,
                               afE!::Function,
                               idxl::Int64,
                               ::Val{k},
                               ::Val{h},
                               ::Val{ny},
                               ::Val{model}) where {k, h, ny, model}

  # n states
  ns = (2^k - 1)*h

  # create individual areas subsets
  S = create_states(k, h)

  # start Expression for ODE equations
  eqs = quote end
  popfirst!(eqs.args)

  # add parameter subset
  push!(eqs.args, :(p = pp[idxl]))

  # add extinction function
  push!(eqs.args, :(afE!(t, rE, pp)))

  # add environmental function
  push!(eqs.args, :(af!(t, r)))

  # compute exponential for each β*covariable
  expex = exp_expr(k, h, ny, model)

  for ex in expex.args
    push!(eqs.args, ex)
  end

  for i in Base.OneTo(ns), si = Base.OneTo(ns)

    # state range
    s = S[si]

    # length of geographic range
    ls = length(s.g)

    # which single areas do not occur in r
    oa = setdiff(1:k, s.g)

    # no events
    nev = noevents_expr(si, s, ls, oa, k, h, ns, i, model)

    # remove if !isone(lr)
    if ls > 1
      # local extinction
      lex = localext_expr(s, S, k, h, i, model)

      # between-region speciation
      brs = brspec_expr(s, S, ns, k, i)
    end

    # dispersal
    # remove if lr == k
    if ls != k
      dis = dispersal_expr(s, oa, S, ns, k, h, i, model)
    end

    # hidden states transitions
    hid = h > 1 ? hidtran_expr(s, S, ns,k, h, i) : :0.0

    # within-region speciation
    wrs = wrspec_expr(si, s, ns, k, i, model)

    # push `D` equation to to eqs
    if isone(ls)
      push!(eqs.args,
        :(du[$si, $i] = $nev + $dis + $hid + $wrs))
    elseif ls == k
      push!(eqs.args,
        :(du[$si, $i] = $nev + $lex + $hid + $wrs + $brs))
    else
      push!(eqs.args,
        :(du[$si, $i] = $nev + $lex + $dis + $hid + $wrs + $brs))
    end

  end

  ## aesthetic touches
  # sort `du`s

  #Core.println(eqs)

  return quote
    @inbounds begin
      $eqs
    end
    return nothing
  end

end

"""
    make_egeohisse_M(::Val{k},
                     ::Val{h},
                     ::Val{ny},
                     ::Val{model},
                     af!::Function) where {k, h, ny, model}

Make closure for EGeoHiSSE given `E(t)` and in matrix form.
"""
function make_egeohisse_M(::Val{k},
                          ::Val{h},
                          ::Val{ny},
                          ::Val{model},
                          af! ::Function,
                          afE!::Function,
                          idxl::Int64) where {k, h, ny, model}

  ns = (2^k - 1)*h

  r    = Array{Float64,1}(undef, ny)
  eaft = Array{Float64,1}(undef,
    model[1]*k*h + model[2]*k*h + model[3]*k*(k-1)*h)
  rE  = Array{Float64,1}(undef,ns)

  # make ode function with closure
  ode_fun = (du::Array{Float64,2},
             u ::Array{Float64,2},
             pp::Array{Array{Float64,1},1},
             t ::Float64) ->
    begin
      geohisse_M(du::Array{Float64,2},
                 u ::Array{Float64,2},
                 pp::Array{Array{Float64,1},1},
                 t ::Float64,
                 r, eaft, rE, af!, afE!, idxl,
                 Val(k), Val(h), Val(ny), Val(model))
      return nothing
    end

  return ode_fun
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
"""
```julia
    loop( iter )
```

Loop over `iter` until it finishes.
"""
function loop( iter )
    x = nothing

    for y in iter
        x = y
    end

    return x
end

ext: jl


text:
# Fast version with BasicCounting algorithm

include("./LinkedList.jl")
#using ListItem

mutable struct AdwinObject
    DELTA                   :: Float64
    mintMinimLongitudWindow :: Int
    mdbldelta               :: Float64
    mintTime                :: Int
    mintClock               :: Int
    mdblWidth               :: Float64  # Mean of Width = mdblWidth/Number of items
    #BUCKET
    MAXBUCKETS              :: Int
    lastBucketRow           :: Int
    TOTAL                   :: Float64
    VARIANCE                :: Float64
    WIDTH                   :: Int
    BucketNumber            :: Int
    #
    Detect                  :: Int
    numberDetections        :: Int
    DetectTwice             :: Int
    blnBucketDeleted        :: Bool
    BucketNumberMAX         :: Int
    mintMinWinLength        :: Int
    listRowBuckets          :: List
    AdwinObject() = new()
end

function initAdwinObject!(AO)
    AO.DELTA = 0.002
    AO.mintMinimLongitudWindow = 10
    AO.mdbldelta = 0.002
    AO.mintTime = 0
    AO.mintClock = 32
    AO.mdblWidth = 0 # Mean of Width = mdblWidth/Number of items
    #BUCKET
    AO.MAXBUCKETS = 5
    AO.lastBucketRow = 0
    AO.TOTAL = 0
    AO.VARIANCE = 0
    AO.WIDTH = 0
    AO.BucketNumber = 0
    AO.Detect = 0
    AO.numberDetections = 0
    AO.DetectTwice  = 0
    AO.blnBucketDeleted = false
    AO.BucketNumberMAX = 0
    AO.mintMinWinLength = 5
    #AO.listRowBuckets = List()
end

function initBuckets!(ao::AdwinObject)
    ao.listRowBuckets = List()
    ao.lastBucketRow = 0
    ao.TOTAL = 0
    ao.VARIANCE = 0
    ao.WIDTH = 0
    ao.BucketNumber = 0
end

function setInput!(ao::AdwinObject, intEntrada::Float64, delta::Float64)
    blnChange = false
    blnExit = false
    cursor :: ListIt = ListIt()
    ao.mintTime += 1
    # 		//1,2)Increment window in one element
    insertElement!(ao, intEntrada) # -- 1
    bitBucketDeleted = false
    #		//3)Reduce  window
    if ao.mintTime % ao.mintClock == 0 && ao.WIDTH > ao.mintMinimLongitudWindow #line 401
      # println("... setInput!: Enter first If")
      blnReduceWidth= true
      while blnReduceWidth
        # println("... setInput!: Enter While")
        blnReduceWidth = false
        blnExit = false # line 406
        n0=0
        n1=ao.WIDTH
        u0=0
        u1= ao.TOTAL # getTotal();
        v0=0
        v1=ao.VARIANCE
        n2=0
        u2=0

        cursor = ao.listRowBuckets.tail
        i = ao.lastBucketRow
        while !blnExit && !(cursor == cursor.next) # line 418
          println("... setInput!: Enter second While")
          for k = 0 : cursor.bucketSizeRow  # line 419
            println(k)
          end
        end
        # println(myCounter(ao.listRowBuckets.head))
      end
    end
end

function setInput!(ao::AdwinObject, intEntrada::Float64)
    setInput!(ao, intEntrada, ao.mdbldelta)
end

function insertElementBucket!(ao::AdwinObject, value::Float64, variance::Float64, node::Node)
    insertBucket!(node, value, variance) # todo -- 3
    ao.BucketNumber += 1
    if ao.BucketNumber > ao.BucketNumberMAX
      ao.BucketNumberMAX = ao.BucketNumber
    end
end

# function compressBuckets!(ao::AdwinObject)
#     # Traverse the list of buckets in increasing order
#     n1::Int
#     n2::Int
#     u1::Float64
#     u2::Float64
#     incVariance::Float64
#     cursor::ListIt
#     nextNode::ListIt
#     cursor = ao.listRowBuckets.head
#     i::Int = 0
#     while !isnull(cursor)
#         # Find the number of buckets in a row
#         k = cursor.bucketSizeRow
#         # If the row is full, merge buckets
#         if k == ao.MAXBUCKETS+1
#             nextNode = cursor.next
#             if isnull(nextNode)
#                 ao.listRowBuckets # line 365
#             end
#         end
#     end
# end

function compressBuckets!(ao::AdwinObject) # --2
  # 		//Traverse the list of buckets in increasing order
  n1::Int=0
  n2::Int=0
  u1::Float64=0
  u2::Float64=0
  cursor :: Node = ao.listRowBuckets.head
  nextNode:: Node = EmptyNode()
  i=0
  while ! typeof(cursor) <: EmptyNode # --3
    # //Find the number of buckets in a row
    #k = cursor.bucketSizeRow
    # 			//If the row is full, merge buckets
    #if k == ao.MAXBUCKETS
    #  nextNode = cursor.next
    #else
    #  break
    #end
    cursor = cursor.next
  end
end

function insertElement!(ao::AdwinObject, value::Float64)
    ao.WIDTH += 1
    insertElementBucket!(ao, 0.0, value, ao.listRowBuckets.head)
    incVariance::Float64 = 0.0
    if ao.WIDTH > 1
        incVariance = (ao.WIDTH-1)*(value-ao.TOTAL/(ao.WIDTH-1))*(value-ao.TOTAL/(ao.WIDTH-1))/ao.WIDTH
    end
    ao.VARIANCE += incVariance
    ao.TOTAL += value
    # println("w: $value")
    # println(myCounter(ao.listRowBuckets.head))
    #compressBuckets!(ao)
end

function run_adwin()
    AO = AdwinObject()
    initAdwinObject!(AO)
    initBuckets!(AO)
    for p in 0:999
        x = p < 500 ? 1000.0 : 500.0
        setInput!(AO, x)
    end
    println("... run_adwin(): OK")
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
module BeastUtils

include("MatrixUtils.jl")
# include("TreeUtils.jl")
# include("Simulation.jl")
# include("BeastNames.jl")
# include("XMLConstructor.jl")
include("Logs.jl")
include("DataStorage.jl")
include("XMLTools.jl")
include("ESS.jl")
include("RunBeast.jl")
include("PosteriorSummary.jl")
include("Nexus.jl")
# include("SlowLikelihoods.jl")

end


text:
"""
Centered content in a card with :content, :caption, :width
"""
centered_content_tpl = """

  {{{:content}}}

      {{{:caption}}}

"""

"""

Take an image file and encode it

Caption uses LaTeX markup, not markdown.

"""
mutable struct ImageFile
    f
    caption
    alt
    width
    content
end
ImageFile(f,caption=""; alt="A Figure", width=nothing) = ImageFile(f, caption, alt, width)
function ImageFile(f, caption, alt, width)
    imgfile = tempname() * ".gif"
    io = open(imgfile, "w")
    show(io, "image/png", f)
    close(io)
    ImageFile(imgfile, caption, alt, width)
end
function ImageFile(f::AbstractString, caption, alt, width)
    data = base64encode(read(f, String))
    content = Mustache.render(gif_to_img_tpl, data=data, alt=alt)
    caption = Markdown.parse(caption)
    ImageFile(f, caption, alt, width, content)
end

" template for an base64 encoded image"
gif_to_img_tpl = """

"""

# """

function gif_to_data(imgfile::AbstractString, caption="", width=480)
    data = base64encode(read(imgfile, String))
    content = Mustache.render(gif_to_img_tpl, data=data, alt="A Figure")
    Mustache.render(centered_content_tpl, content=content, caption=md(caption), width=width)
end

# take an image and process
function gif_to_data(img, caption="")
    imgfile = tempname() * ".gif"
    io = open(imgfile, "w")
    show(io, "image/png", img)
    close(io)
    gif_to_data(imgfile, caption)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "sinewave"
version = v"0.3.0"

# Collection of sources required to complete build
sources = [
    "https://github.com/jakubwro/sinewave.git" =>
    "193b8b776d26cac02cabd264e03b917a928769be",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd sinewave
export CPPFLAGS="-I${prefix}/include"
export LDFLAGS="-L${libdir}"
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libsinewave", :libsinewave)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "FFTW_jll"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Distances
using Clustering
using MatrixNetworks
using ImageSegmentation
using SparseArrays
export
    edge_weight,
    boundary_map,
    Precision,
    FMeasure,
    FBoundary,
    SegmentationCovering,
    VariationOfInformation,
    RandIndex,
    FMeasureRegions,
    PRObjectsAndParts,
    BoundaryDisplacementError
using MLBase

function evaluate(f, cl_seg::SegmentedImage, gt_seg::SegmentedImage)
    return evaluate(f, labels_map(cl_seg), labels_map(gt_seg))
end
function evaluate(f, cl_seg::Matrix{T1}, gt_seg::Matrix{T2}) where T1<:Integer where T2<:Integer
    return evaluate(f, convert.(Integer, cl_seg), convert.(Integer, gt_seg))
end

"""
```
struct FBoundary
```

Learning to Detect Natural Image Boundaries Using Local Brightness, Color, and Texture Cues
David R. Martin, Member, IEEE, Charless C. Fowlkes, and Jitendra Malik, Member, IEEE

# Members
- `dmax::Float64`
"""
struct FBoundary
    dmax::Float64
    FBoundary(;dmax::Float64=0.0075) = new(dmax)
end

function pixel_mapping(cl::AbstractMatrix, gt::AbstractMatrix, dmax)
    pixToIdxCL = Dict{Tuple{Integer,Integer},Integer}()
    pixToIdxGT = Dict{Tuple{Integer,Integer},Integer}()
    idxToPixCL = Vector{Tuple{Integer,Integer}}()
    idxToPixGT = Vector{Tuple{Integer,Integer}}()
    ee = Vector()
    matchableGT = falses(size(cl))
    matchableCL = falses(size(cl))
    nCL = 0
    nGT = 0
    r = round(Integer,ceil(dmax))
    for j in findall(cl)
        (rc,cc) = Tuple(CartesianIndices(size(cl))[j])
        row_range = max(rc-r,1):min(rc+r, size(cl,1))
        col_range = max(cc-r,1):min(cc+r,size(cl,2))
        for CI in CartesianIndices((row_range, col_range))
            if gt[CI]
                k =  LinearIndices(size(gt))[CI[1],CI[2]]
                dist = sqrt((CI[1]-rc)^2 +  (CI[2]-cc)^2)
                if dist <= dmax
                    matchableCL[rc,cc] = true
                    matchableGT[CI] = true
                    push!(ee, ( ( (rc,cc),  (CI[1], CI[2])), round(Integer,dist*100)))
                end
            end
        end
    end
    for J in CartesianIndices(size(cl))
        if (matchableCL[J])
            nCL = nCL +1
            pixToIdxCL[(J[1],J[2])] = nCL
            push!(idxToPixCL, (J[1],J[2]))
        end
        if matchableGT[J]
            nGT = nGT +1
            pixToIdxGT[(J[1],J[2])] = nGT
            push!(idxToPixGT, (J[1], J[2]))
        end
    end
    return (pixToIdxCL, idxToPixCL, nCL, matchableCL, pixToIdxGT, idxToPixGT, nGT, matchableGT, ee)
end
function evaluate(f::FBoundary, cl_seg::Matrix{T1}, gt_seg::Matrix{T2}) where T1<:Integer where T2<:Integer

     cl = thinning(boundary_map(BoundaryGradient(), cl_seg))
     gt = thinning(boundary_map(BoundaryGradient(), gt_seg))
     scale_cost =  sqrt( size(cl,1)^2 + size(cl,2)^2)
     dmax = f.dmax * scale_cost
     noutliers = 6
     cl_detected = findall(cl)
     gt_detected = findall(gt)
     outlier_weight =-100*f.dmax*scale_cost

    (pixToIdxCL, idxToPixCL, nCL, matchableCL, pixToIdxGT, idxToPixGT, nGT, matchableGT, ee) = pixel_mapping(cl, gt, dmax)

     I = sizehint!(Vector{Int64}(),(nCL+nGT)*6)
     J = sizehint!(Vector{Int64}(),(nCL+nGT)*6)
     w     = sizehint!(Vector{Float64}(),(nCL+nGT)*6)

    Aidx(x) = Integer(x)
    Aoutlier(x) = Integer(nCL+x)
    Bidx(x) = Integer(nCL+nGT+x)
    Boutlier(x) = Integer(nCL+nGT+nGT+x)

     n = 2*(nCL+nGT)
    for (edge, weight) in ee
        push!(I, Aidx(pixToIdxCL[edge[1]]))
        push!(J, Bidx(pixToIdxGT[edge[2]]))
        push!(w    , weight)

    end
    for i=1:nCL
        for j in rand(1:nCL-1, noutliers)
            if (i==j)
                j=j+1
            end

            push!(I,Aidx(i))
            push!(J, Boutlier(j))
            push!(w, outlier_weight)
        end
    end

    for i=1:nGT
        for j in rand(1:nGT-1, noutliers)
            if (i==j)
                j=j+1
            end

            push!(I,Aoutlier(j))
            push!(J,Bidx(i))
            push!(w, outlier_weight)
        end
    end
    for i=1:max(nGT,nCL)
        for j in rand(1:min(nGT,nCL), noutliers)
            if (nGT<nCL)

                push!(I, Aoutlier(j))
                push!(J,Boutlier(i))
            else
                push!(I, Aoutlier(i))
                push!(J, Boutlier(j))
            end
            push!(w, outlier_weight)
        end
    end

    for i=1:nCL
        push!(I,Aidx(i))
        push!(J,Boutlier(i))
        push!(w, outlier_weight)
    end

    for i=1:nGT
        push!(I, Aoutlier(i))
        push!(J,Bidx(i))
        push!(w, outlier_weight)
    end

     match = bipartite_matching(sparse(I,J,w,n,n))
    (cl_edge, gt_edge) = edge_list(match)
     matchCL = zeros(Integer,size(cl))
     matchG = zeros(Integer,size(cl))
     ff = 0
    for j=1:length(cl_edge)
         v1 = cl_edge[j]
         v2 = gt_edge[j] - (nCL+nGT)
        if (v1 > nCL+nGT)
             t = v1
            v1 = v2
            v2 = t
        end
        if (v1<= nCL) &&  (v2<=nGT)
             pix1 = idxToPixCL[v1]
            if cl[pix1[1],pix1[2]]
                matchCL[pix1[1], pix1[2]] = v2
            end

             pix2 = idxToPixGT[v2]
            if gt[pix2[1],pix2[2]]

                matchG[pix2[1],pix2[2]] = v1
            end
        end
    end

    cntR = sum(matchG.>0)
    sumR = length(gt_detected)
    cntP = sum(matchCL.>0)
    sumP = length(cl_detected)
    rec = cntR/(sumR +  eps());
    prec = cntP/(sumP + eps());
    if sumR == 0
        @warn("Infinit recall")
    end
    if sumP == 0
        @warn("Infinit precision")
    end
    return ((2*prec*rec)/(prec+rec + eps()), prec, rec)
end

function evaluate(d, cl::Matrix{T}, gt::Matrix{T}) where T<:Integer
     N = prod(size(cl))
     clusters_cl = [(cl.==j) for j in unique(cl)]
     clusters_gt = [(gt.==j) for j in unique(gt)]

     sum_tot = 0
    for c in clusters_cl
         length = sum(c)
        sum_tot = sum_tot + maximum([length*edge_weight(d, c,g) for g in clusters_gt])
    end
    return sum_tot / N

end
"""
Precision
"""
struct Precision
end
function edge_weight(d::Precision, cl::BitArray{2}, gt::BitArray{2})
    return precision(roc(vec(gt),vec(cl)))
end
"""
FMeasure
"""
struct FMeasure
end
function edge_weight(d::FMeasure, cl::BitArray{2}, gt::BitArray{2})
    return f1score(roc(vec(gt),vec(cl)))
end
"""
Segmentation Covering
"""
struct SegmentationCovering end
function edge_weight(d::SegmentationCovering, cl::BitArray{2}, gt::BitArray{2})
    return 1-jaccard(Integer.(vec(gt)),Integer.(vec(cl)))
end
"""
Variation of Information
"""
struct VariationOfInformation
    normalize::Bool
end
function VariationOfInformation()
    return VariationOfInformation(true)
end

function H(clusters::Vector{BitArray{2}})
     n = prod(size(clusters[1]))
     tot_sum = 0
    for c in clusters
         R = sum(c)
        tot_sum = tot_sum + R*log(R/n)
    end
    return (-1/n)*tot_sum
end
function I(clusters1::Vector{BitArray{2}}, clusters2::Vector{BitArray{2}})
     tot_sum = 0
     n = prod(size(clusters1[1]))
    for c1 in clusters1
         R1 = sum(c1)
        for c2 in clusters2
             R2 = sum(c2)
             intersection = sum(c1 .& c2)
            if (intersection>0)
                tot_sum = tot_sum + intersection *log(((n*intersection)/(R1*R2)))
            end
        end
    end
    return (1/n)*tot_sum

end
function evaluate(d::VariationOfInformation, cl::Matrix{T}, gt::Matrix{T}) where T<:Integer
    ff = varinfo(maximum(cl), vec(cl), maximum(gt),vec(gt))
    clusters_cl = [(cl.==j) for j in unique(cl)]
    clusters_gt = [(gt.==j) for j in unique(gt)]
    voi =  H(clusters_cl) + H(clusters_gt) - 2*I(clusters_cl,clusters_gt)
    if (d.normalize)
        return voi / 2*log(max(maximum(cl),maximum(gt)))
    else
        return voi
    end
end
"""
RandIndex
"""
struct RandIndex end
function evaluate(d::RandIndex, c1::Matrix{T}, gt::Matrix{T}) where T<:Integer
     n = prod(size(c1))
     M = zeros(2,2)
    (a,b,c,d) = randindex(vec(c1), vec(gt))
    return b
end
"""
FMeasure Regions
"""
struct FMeasureRegions end

""""
```julia
function relabel(c1::Matrix}, part_bimap=Dict{Integer,Integer}()) where T<:Integer
```
It relabels a partition in scanning order. Bimaps are the look up tables of the relabeling.
  - author Jordi Pont Tuset <jordi.pont@upc.edu>
"""
function relabel(c1::Matrix{T}, bimap=Dict{Integer,Integer}()) where T<:Integer
     partition_out = zeros(Integer, size(c1))
     max_region = 1;
    for I in CartesianIndices(size(c1))
        if !haskey(bimap, c1[I])
            bimap[c1[I]] = max_region;
            partition_out[I] = max_region;
            max_region=max_region+1;
        else
            partition_out[I] = bimap[c1[I]]
        end
    end
    return (partition_out, max_region-1)
end
"""
```julia
function evaluate(d::FMeasureRegions, c1::Matrix{T}, gt::Matrix{T}) where T<:Integer
```

# References
- https://github.com/jponttuset/seism/blob/c0f10b4509d364e39e5dba28fe57c2d157c0808e/src/misc/intersection_matrix.hpp
- https://github.com/JuliaStats/Clustering.jl/blob/master/src/randindex.jl
"""

function evaluate(d::FMeasureRegions, c1::Matrix{T}, gt::Matrix{T}) where T<:Integer

    @assert(size(c1)==size(gt))

    (partition1_relab, num_reg_1) = relabel(c1)
    (partition2_relab, num_reg_2) = relabel(gt)
    c = Clustering.counts(partition1_relab,
                          partition2_relab,
                          (1:maximum(partition1_relab),
                           1:maximum(partition2_relab))) # form contingency matrix
    n = round(Int,sum(c))
    nis = sum(sum(c, dims=2).^2)        # sum of squares of sums of rows
    njs = sum(sum(c, dims=1).^2)        # sum of squares of sums of columns
    t2 = sum(c.^2)                # sum over rows & columnns of nij^2
    t3 = .5*(nis+njs)
    n11 = (t2-n)/2
    n00 = (n^2 - nis -njs +t2)/2
    n10 = (nis-t2)/2
    n01 = (njs-t2)/2

    precision = n11/(n11+n10);
    recall    = n11/(n11+n01);
    if (precision+recall>0)
        return (2*precision*recall/(precision+recall), precision, recall)
    else
        return (0,0,0)
    end
end

"""

http://www.cv-foundation.org/openaccess/content_cvpr_2013/papers/Pont-Tuset_Measures_and_Meta-Measures_2013_CVPR_paper.pdf

Measures and Meta-Measures for the Supervised Evaluation of Image Segmentation Jordi Pont-Tuset and Ferran Marques. Universitat Politecnica de Catalunya BarcelonaTech

"""
struct PRObjectsAndParts
    object_threshold::Float64
    part_threshold::Float64
    B::Float64
    PRObjectsAndParts(;
                      object_threshold::Float64=0.95,
                      part_threshold::Float64=0.25,
                      B::Float64=0.1) = new(object_threshold,
                                            part_threshold,
                                            B)
end

function calculate_regions(num_reg_gt::Integer, num_reg_part::Integer,  intersect_matrix::Matrix)
     image_area        = 0;
     region_areas_gt   = zeros(Integer,num_reg_gt)
     region_areas_part = zeros(Integer, num_reg_part);
    for ii=1:num_reg_gt
        for jj=1:num_reg_part
            image_area += intersect_matrix[jj,ii];
            region_areas_part[jj] = region_areas_part[jj] + intersect_matrix[jj,ii];
            region_areas_gt[ii]   = region_areas_gt[ii]   + intersect_matrix[jj,ii];
        end
    end
    return (region_areas_gt, region_areas_part, image_area)
end
function get_candidates(image_area::Integer, region_areas::Vector)

     area_percentile = 0.99;
     area_map  = Vector{Tuple{Float64, Integer}}() # Mapping between each region area and its id
     candidates = Dict{Integer,Bool}()
    #Get candidates in the partition (remove percentile of small area)
    for ii=1:length(region_areas)
        push!(area_map,((region_areas[ii])/float(image_area),ii));
    end
    sort!(area_map, by=x->x[1],rev=true)
     curr_pct = 0;
    for (area_prop, area) in area_map
        if (curr_pct < area_percentile)
            candidates[area] = true
        else
            candidates[area] = false
        end
        curr_pct += area_prop
    end
    return candidates
end
function evaluate(c::PRObjectsAndParts,
                  assignments_seg::Matrix{T},
                  assignments_gt::Matrix{T}) where T<:Integer

    NOT_CLASSIFIED = 1
    OBJECT         = 2
    PART           = 3

    assignments_seg = relabel(assignments_seg)[1]
    assignments_gt  = relabel(assignments_gt)[1]
     intersect_matrix = Clustering.counts(assignments_seg, assignments_gt,(1:maximum(assignments_seg),1:maximum(assignments_gt))) # form contingency matrix

     num_reg_part        = size(intersect_matrix,1)
     classification_part = ones(Integer,num_reg_part)
     prec_part           = zeros(num_reg_part);
     mapping_gt          = Dict{Integer,Integer}()
     mapping_part        = Dict{Integer,Integer}()
     num_reg_gt          = size(intersect_matrix,2)
     classification_gt   = ones(Integer, num_reg_gt)
     recall_gt           = zeros(num_reg_gt);

    (region_areas_gt, region_areas_part, image_area) =  calculate_regions(num_reg_gt, num_reg_part, intersect_matrix)

     candidate_part = get_candidates(image_area, region_areas_part)
     candidate_gt   = get_candidates(image_area, region_areas_gt)

    # Scan through table and find all OBJECT mappings */
    for ii=1:num_reg_gt
        for jj=1:num_reg_part
             recall    = intersect_matrix[jj,ii]/float(region_areas_gt[ii]);
             precision = intersect_matrix[jj,ii]/float(region_areas_part[jj]);
            # Ignore those regions with tiny area */
            if(candidate_gt[ii]==true && candidate_part[jj]==true)
                # Is it an object candidate? */
                if(recall >= c.object_threshold)  &&  (precision >= c.object_threshold)
                    classification_gt[ii]  = OBJECT;
                    classification_part[jj] = OBJECT;
                    mapping_gt[ii] = jj;
                    mapping_part[jj]        = ii;
                elseif (recall >= c.part_threshold)  &&  (precision >= c.object_threshold)
                    if (classification_part[jj] == NOT_CLASSIFIED)
                        classification_part[jj] = PART;
                        mapping_part[jj]        = ii;
                    end
                elseif(recall >= c.object_threshold)  &&  (precision >= c.part_threshold)
                        # Cannot have a classification already */
                        classification_gt[ii] = PART;
                        mapping_gt[ii] = jj;
                end
            end

            #Get _recall_gt and _prec_part (no matter if candidates or not), discarding objects
            if(precision >= c.object_threshold) && (recall < c.object_threshold)
                recall_gt[ii] += recall;
            elseif (recall >= c.object_threshold) && (precision < c.object_threshold)
                prec_part[jj] += precision;
            end
        end
    end

    #Count everything
     num_objects_part = 0;
     num_objects_gt = 0;
     num_parts_part = 0;
     num_parts_gt = 0;
     num_underseg_part = 0;
     num_overseg_gt = 0;
     num_candidates_part = 0;
     num_candidates_gt = 0;

    for jj=1:num_reg_part
        num_candidates_part += candidate_part[jj];

        if (classification_part[jj]==PART)
            num_parts_part = num_parts_part+1;
        elseif(classification_part[jj]==OBJECT)
            num_objects_part = num_objects_part +1 ;
        elseif (candidate_part[jj]) # Compute degree of undersegmentation
            num_underseg_part += prec_part[jj];
        end
    end
    for ii=1:num_reg_gt
        num_candidates_gt += candidate_gt[ii];
        if (classification_gt[ii]==PART)
                num_parts_gt = num_parts_gt+1;
        elseif(classification_gt[ii]==OBJECT)
                num_objects_gt = num_objects_gt+1;
        elseif(candidate_gt[ii])
                num_overseg_gt =  num_overseg_gt + recall_gt[ii];
        end
    end
    # Precision and recall
     precision = (num_objects_part + num_underseg_part + c.B*num_parts_part)/float(num_candidates_part);
     recall    = (num_objects_gt   + num_overseg_gt    + c.B*num_parts_gt  )/float(num_candidates_gt);

    # F-measure for Region Detection
    if(precision==0 && recall==0)
        f_measure = (0,0,0) ;
    else
        f_measure = (2*precision*recall/(precision+recall), precision, recall)
    end
end

"""
Boundary Displacement Error
"""
struct BoundaryDisplacementError
end
using Images
function evaluate(cfg::BoundaryDisplacementError, cl::Matrix{T}, gt::Matrix{T})  where T<:Integer
    @assert(size(cl)==size(gt))
    # Generate boundary maps
    boundary1 = boundary_map(BoundaryGradient(), cl)
    boundary2 = boundary_map(BoundaryGradient(), gt)

    # boundary1 and boundary2 are now binary boundary masks. compute their distance transforms:
    D1 = distance_transform(feature_transform(boundary1));
    D2 = distance_transform(feature_transform(boundary2));

    # compute the distance of the pixels in boundary1 to the nearest pixel in% boundary2:
    dist_12 = sum(boundary1 .* D2 );
    dist_21 = sum(boundary2 .* D1 );
    avgError_12 = dist_12 / sum(boundary1);
    avgError_21 = dist_21 / sum(boundary2);
    return (avgError_12 + avgError_21) / 2;
end


text:
# =====================================================================
function update_ff!(ff::Array{Complex{Float32},1},
                    zp::Union{PZResp, PZResp64},
                    f::Array{Float32,1})
  Z = zp.z
  P = zp.p

  i = 0
  γ = 0.0f0
  g = 0.0f0
  d = zero(Complex{Float32})
  ϵ = eps(Float32)^2

  @inbounds while i < length(f)
    i = i + 1
    cf = im*f[i]
    n = 1.0f0
    for z in Z
      n *= cf-z
    end
    for p in P
      d = conj(cf-p)
      n *= d
      n /= max(ϵ, abs2(d))
    end
    g = abs2(n)
    if g > γ
      γ = g
    end
    setindex!(ff, n, i)
  end

  rmul!(ff, zp.a0)

  return sqrt(γ)*zp.a0
end

# frequencies in radians/s from -pi*fs to pi*fs
function fill_f!(f::Array{Float32,1}, fs::Float32, N2::Int64)
  N = div(N2, 2)
  fn = Float32(pi*fs)
  df = fn/Float32(N)
  f[:] .= 0:N2-1
  for i = N+2:N2
    f[i] -= N2
  end
  rmul!(f, df)
  return f
end

function update_resp!(f::AbstractArray,
                      ff_old::AbstractArray,
                      ff_new::AbstractArray,
                      N2::Int64,
                      fs::Float32,
                      resp_old::Union{PZResp, PZResp64},
                      resp_new::Union{PZResp, PZResp64},
                      wl::Float32)

  fill_f!(f, fs, N2)
  update_ff!(ff_new, resp_new, f)
  γ = update_ff!(ff_old, resp_old, f)

  # allows manual watermarking for ill-behaved translations with scaling isssues
  wm = γ * wl
  for i = 1:N2
    setindex!(ff_new, (ff_new[i]*conj(ff_old[i])) / (abs2(ff_old[i]) + wm), i)
  end

  return nothing
end

function update_resp_vecs!( Xw::Array{Complex{Float32},1},
                            f::Array{Float32,1},
                            ff_old::Array{Complex{Float32},1},
                            ff_new::Array{Complex{Float32},1},
                            N2::Int64 )
  resize!(Xw, N2)
  resize!(f, N2)
  resize!(ff_old, N2)
  resize!(ff_new, N2)
  xfl = reinterpret(Float32, Xw)
  xre = view(xfl, 1:2:2*N2-1)
  return xfl, xre
end


text:
function [Grid] = FacesInEdges(Grid)
NumFacesPerEdge=zeros(Grid.NumEdges,1);
for iF=1:Grid.NumFaces
  Face=Grid.Faces(iF);
  for iE=1:size(Face.E,2)
    NumFacesPerEdge(Face.E(iE))=NumFacesPerEdge(Face.E(iE))+1;
  end
end
FacesPerEdge=zeros(Grid.NumEdges,max(NumFacesPerEdge));
NumFacesPerEdge=zeros(Grid.NumEdges,1);
for iF=1:Grid.NumFaces
  Face=Grid.Faces(iF);
  for iE=1:size(Face.E,2)
    NumFacesPerEdge(Face.E(iE))=NumFacesPerEdge(Face.E(iE))+1;
    FacesPerEdge(Face.E(iE),NumFacesPerEdge(Face.E(iE)))=iF;
  end
end
for iE=1:Grid.NumEdges
  Grid.Edges(iE).F=FacesPerEdge(iE,1:NumFacesPerEdge(iE));
end
end


text:
#
# This file is part of module
#   GenerateMediumDict (ModiaMedia/dict/GenerateMediumDict.jl)
#

dict["ConstantPropertyLiquidWater"] = ModiaMedia.SimpleMedium(
                        mediumName = "ConstantPropertyLiquidWater",
                        fluidConstants = ModiaMedia.BasicFluidConstants(
                            chemicalFormula="H2O",
                            structureFormula="H2O",
                            casRegistryNumber="7732-18-5",
                            iupacName="oxidane",
                            molarMass=0.018015268),
                        data = ModiaMedia.SimpleMediumData(
                            cp_const=4184,
                            cv_const=4184,
                            d_const=995.586,
                            eta_const=1.e-3,
                            lambda_const=0.598,
                            a_const=1484,
                            T_min=ModiaMedia.from_degC(-1),
                            T_max=ModiaMedia.from_degC(130),
                            T0=273.15,
                            MM_const=0.018015268)
                    )


text:
@testset "FluxExtensions.triangularloss" begin
	d = rand(10)
	y = rand(1:3,10)
	g = Flux.Tracker.ngradient((d) -> FluxExtensions.triangularloss(d,y), d)[1];
	@test all( g - FluxExtensions.triangularloss_back(d, y, 1.0) .== 0 )
	d = Flux.param(d)
	Flux.back!(FluxExtensions.triangularloss(d,y))
	@test all(d.grad - g .== 0)
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
# Based on https://github.com/cloud-oak/Tqdm.jl by @cloud-oak under Mozilla Public License 2.0
# Modified for Knet by Deniz Yuret

using Printf
import Base: length, size, iterate, eltype, IteratorSize, IteratorEltype, haslength, @propagate_inbounds

"""
    progress(itr; width, alpha, interval)

Return an iterator which acts exactly like `itr`, but prints a progressbar as new values are
requested:

    2.70e-01  21.83%┣███▉              ┫ 13101/60000 [00:12/00:53, 1137.05i/s]

Here `2.70e-01` is the exponential average of values generated by itr (only displayed for
iterators with numeric values). `21.83%` is the percentage, `13101` is the number of
iterations completed, `60000` is the total number of iterations. `00:12` is elapsed seconds,
`00:53` is the estimated total seconds, `1137.05i/s` is the average number of iterations
completed per second. If the speed is less than 1, the average number of seconds per
iteration (s/i) is reported instead.  The percent, total iterations, and completion time are
omitted for iterators whose size is unknown.

`progress!(...)` is equivalent to `(for x in progress(...) end)`, i.e.  iterates over the
object created by `progress(...)` and returns `nothing`.

An integer `itr` is treated as `1:itr`, i.e. `progress(n::Integer)` is equivalent to
`progress(1:n)`

Keyword arguments:
* `width=max(64,displaysize()[2])`: controls display width. The default width can be
  controlled using `ENV["COLUMNS"]`.
* `interval=1.0`: minimum time interval in seconds between progressbar updates.
* `alpha=1.0`: controls the exponential average displayed for numeric iterators:
  `avg = alpha * val + (1-alpha) * avg`
"""
progress, progress!

mutable struct Progress{I}
    iter::I
    current::Int
    nprint::Int
    start_time::UInt
    print_time::UInt
    print_interval::UInt
    width::Int
    alpha::Float64
    avg::Float64
end

progress(iter::I; width=max(64,displaysize()[2]), alpha=1.0, interval=1.0) where {I} =
    Progress{I}(iter,0,0,time_ns(),0,Int(1e9*interval),width,alpha,Inf)

progress(i::Integer; o...)=progress(1:n; o...)
progress!(x...; o...)=(for x in progress(x...; o...) end)

length(p::Progress) = length(p.iter)
size(p::Progress) = size(p.iter)
eltype(p::Progress) = eltype(p.iter)
IteratorSize(::Type{Progress{I}}) where {I} = IteratorSize(I)
IteratorEltype(::Type{Progress{I}}) where {I} = IteratorEltype(I)

@propagate_inbounds function iterate(p::Progress, s...)
    next = iterate(p.iter, s...)
    if next !== nothing
        p.current += 1
        (x, s) = next
        if p.alpha > 0 && value(x) isa Number
            p.avg = (p.avg === Inf ? value(x) : p.alpha * value(x) + (1-p.alpha) * p.avg)
        end
    end
    display_progress(p, next === nothing)
    return next
end

function display_progress(p::Progress, last=false)
    curr_time = time_ns()
    if !last && (curr_time < p.print_time + p.print_interval)
        return
    end
    p.print_time = curr_time
    p.nprint += 1
    seconds   = (curr_time - p.start_time) * 1e-9
    speed     = p.current / seconds

    fval_string = isfinite(p.avg) ? @sprintf("%.2e  ", p.avg) : ""

    if haslength(p)
        ETA = (length(p) - p.current) / speed
        percentage_string = string(@sprintf("%.2f%%",p.current/length(p)*100))
        status_string = string(p.current, "/", length(p),
                               " [", format_time(seconds), "/", format_time(seconds+ETA),
                               ", ", format_speed(speed),"]")

    else
        ETA = Inf
        percentage_string = ""
        status_string = string(p.current, # "/", length(p),
                               " [", format_time(seconds), # "-", format_time(ETA),
                               ", ", format_speed(speed),"]")
    end

    print("\r")
    print(fval_string)
    print(percentage_string)
    print("┣")

    width = p.width - length(fval_string) - length(percentage_string) - length(status_string) - 2

    if (haslength(p))
        cellvalue = length(p) / width
        full_cells, remain = divrem(p.current, cellvalue)
        full_cells = round(Int, full_cells)
        print(repeat("█", full_cells))
        if (full_cells < width)
	    part = floor(Int, 8 * remain / cellvalue)
	    print(EIGHTS[part])
            print(repeat(" ", width - full_cells - 1))
        end
    else
        offset = p.nprint % 8
        print(repeat(" ", offset))
        print("/")
        segments, remain = divrem(width - offset - 1, 8)
        print(repeat("       /", Int(segments)))
        print(repeat(" ", Int(remain)))
    end
    print("┫ ")

    print(status_string)
    last && println()
end

function format_time(seconds)
    if seconds != Inf
        mins,s  = divrem(round(Int,seconds), 60)
        h, m    = divrem(mins, 60)
    else
        h=0;m=Inf;s=Inf
    end
    if h!=0
         return @sprintf("%02d:%02d:%02d",h,m,s)
    else
         return @sprintf("%02d:%02d",m,s)
    end
end

format_speed(s)=(s >= 1 ? @sprintf("%.2fi/s",s) : @sprintf("%.2fs/i",1/s))

EIGHTS = Dict(0 => ' ',
	      1 => '▏',
	      2 => '▎',
	      3 => '▍',
	      4 => '▌',
	      5 => '▋',
	      6 => '▊',
	      7 => '▉',
	      8 => '█')


text:
export lsqr, lsqr!

# Adapted from the BSD-licensed Matlab implementation at
#    http://www.stanford.edu/group/SOL/software/lsqr.html
#
# LSQR solves Ax = b or min ||b - Ax||_2 if damp = 0,
# or   min ||(b) - (  A   )x||   otherwise.
#          ||(0)   (damp*I) ||_2

#-----------------------------------------------------------------------
# LSQR uses an iterative (conjugate-gradient-like) method.
# For further information, see
# 1. C. C. Paige and M. A. Saunders (1982a).
#    LSQR: An algorithm for sparse linear equations and sparse least squares,
#    ACM TOMS 8(1), 43-71.
# 2. C. C. Paige and M. A. Saunders (1982b).
#    Algorithm 583.  LSQR: Sparse linear equations and least squares problems,
#    ACM TOMS 8(2), 195-209.
# 3. M. A. Saunders (1995).  Solution of sparse rectangular systems using
#    LSQR and CRAIG, BIT 35, 588-604.
#
# Input parameters:
# atol, btol  are stopping tolerances.  If both are 1.0e-9 (say),
#             the final residual norm should be accurate to about 9 digits.
#             (The final x will usually have fewer correct digits,
#             depending on cond(A) and the size of damp.)
# conlim      is also a stopping tolerance.  lsqr terminates if an estimate
#             of cond(A) exceeds conlim.  For compatible systems Ax = b,
#             conlim could be as large as 1.0e+12 (say).  For least-squares
#             problems, conlim should be less than 1.0e+8.
#             Maximum precision can be obtained by setting
#             atol = btol = conlim = zero, but the number of iterations
#             may then be excessive.
# maxiter      is an explicit limit on iterations (for safety).
#
#              Michael Saunders, Systems Optimization Laboratory,
#              Dept of MS&E, Stanford University.
#
# Adapted for Julia by Timothy E. Holy with the following changes:
#    - Allow an initial guess for x
#    - Eliminate printing
#-----------------------------------------------------------------------
function lsqr_method!(log::ConvergenceHistory, x, A, b; damp=0, atol=sqrt(eps(Adivtype(A,b))), btol=sqrt(eps(Adivtype(A,b))), conlim=real(one(Adivtype(A,b)))/sqrt(eps(Adivtype(A,b))), maxiter::Int=max(size(A,1), size(A,2)))
    # Sanity-checking
    m = size(A,1)
    n = size(A,2)
    length(x) == n || error("x should be of length ", n)
    length(b) == m || error("b should be of length ", m)
    for i = 1:n
        isfinite(x[i]) || error("Initial guess for x must be finite")
    end

    # Initialize
    T = Adivtype(A, b)
    Tr = real(T)
    itn = istop = 0
    ctol = conlim > 0 ? convert(Tr, 1/conlim) : zero(Tr)
    Anorm = Acond = ddnorm = res2 = xnorm = xxnorm = z = sn2 = zero(Tr)
    cs2 = -one(Tr)
    dampsq = abs2(damp)
    tmpm = similar(b, T, m)
    tmpn = similar(x, T, n)

    log[:atol] = atol
    log[:btol] = btol
    log[:ctol] = ctol

    # Set up the first vectors u and v for the bidiagonalization.
    # These satisfy  beta*u = b-A*x,  alpha*v = A'u.
    u = b - A*x
    v = copy(x)
    beta = norm(u)
    alpha = zero(Tr)
    if beta > 0
        log.mtvps=1
        scale!(u, inv(beta))
        Ac_mul_B!(v,A,u)
        alpha = norm(v)
    end
    if alpha > 0
        scale!(v, inv(alpha))
    end
    w = copy(v)
    wrho = similar(w)

    Arnorm = alpha*beta
    if Arnorm == 0
        return
    end

    rhobar = alpha
    phibar = bnorm = rnorm = r1norm = r2norm = beta

    #------------------------------------------------------------------
    #     Main iteration loop.
    #------------------------------------------------------------------
    while itn < maxiter
        nextiter!(log,mvps=1)
        itn += 1

        # Perform the next step of the bidiagonalization to obtain the
        # next beta, u, alpha, v.  These satisfy the relations
        #      beta*u  =  A*v  - alpha*u,
        #      alpha*v  =  A'*u - beta*v.

        # Note that the following three lines are a band aid for a GEMM: X: C := αAB + βC.
        # This is already supported in A_mul_B! for sparse and distributed matrices, but not yet dense
        A_mul_B!(tmpm, A, v)
        scale!(u, -alpha)
        LinAlg.axpy!(one(eltype(tmpm)), tmpm, u)
        beta = norm(u)
        if beta > 0
            log.mtvps+=1
            scale!(u, inv(beta))
            Anorm = sqrt(abs2(Anorm) + abs2(alpha) + abs2(beta) + dampsq)
            # Note that the following three lines are a band aid for a GEMM: X: C := αA'B + βC.
            # This is already supported in Ac_mul_B! for sparse and distributed matrices, but not yet dense
            Ac_mul_B!(tmpn, A, u)
            scale!(v, -beta)
            LinAlg.axpy!(one(eltype(tmpn)), tmpn, v)
            alpha  = norm(v)
            if alpha > 0
                scale!(v, inv(alpha))
            end
        end

        # Use a plane rotation to eliminate the damping parameter.
        # This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
        rhobar1 = sqrt(abs2(rhobar) + dampsq)
        cs1     = rhobar/rhobar1
        sn1     = damp  /rhobar1
        psi     = sn1*phibar
        phibar  = cs1*phibar

        # Use a plane rotation to eliminate the subdiagonal element (beta)
        # of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
        rho     =   sqrt(abs2(rhobar1) + abs2(beta))
        cs      =   rhobar1/rho
        sn      =   beta   /rho
        theta   =   sn*alpha
        rhobar  = - cs*alpha
        phi     =   cs*phibar
        phibar  =   sn*phibar
        tau     =   sn*phi

        # Update x and w
        t1      =   phi  /rho
        t2      = - theta/rho

        LinAlg.axpy!(t1, w, x)
        scale!(w, t2)
        LinAlg.axpy!(one(t2), v, w)
        copy!(wrho, w)
        scale!(wrho, inv(rho))
        ddnorm += norm(wrho)

        # Use a plane rotation on the right to eliminate the
        # super-diagonal element (theta) of the upper-bidiagonal matrix.
        # Then use the result to estimate  norm(x).
        delta   =   sn2*rho
        gambar  =  -cs2*rho
        rhs     =   phi - delta*z
        zbar    =   rhs/gambar
        xnorm   =   sqrt(xxnorm + abs2(zbar))
        gamma   =   sqrt(abs2(gambar) + abs2(theta))
        cs2     =   gambar/gamma
        sn2     =   theta /gamma
        z       =   rhs   /gamma
        xxnorm +=   abs2(z)

        # Test for convergence.
        # First, estimate the condition of the matrix  Abar,
        # and the norms of  rbar  and  Abar'rbar.
        Acond   =   Anorm*sqrt(ddnorm)
        res1    =   abs2(phibar)
        res2    =   res2 + abs2(psi)
        rnorm   =   sqrt(res1 + res2)
        Arnorm  =   alpha*abs(tau)

        # 07 Aug 2002:
        # Distinguish between
        #    r1norm = ||b - Ax|| and
        #    r2norm = rnorm in current code
        #           = sqrt(r1norm^2 + damp^2*||x||^2).
        #    Estimate r1norm from
        #    r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
        # Although there is cancellation, it might be accurate enough.
        r1sq    =   abs2(rnorm) - dampsq*xxnorm
        r1norm  =   sqrt(abs(r1sq));   if r1sq < 0 r1norm = - r1norm; end
        r2norm  =   rnorm
        push!(log, :resnorm, r1norm)

        # Now use these norms to estimate certain other quantities,
        # some of which will be small near a solution.
        test1   =   rnorm /bnorm
        test2   =   Arnorm/(Anorm*rnorm)
        test3   =   inv(Acond)
        t1      =   test1/(1 + Anorm*xnorm/bnorm)
        rtol    =   btol + atol*Anorm*xnorm/bnorm
        push!(log, :cnorm, test3)
        push!(log, :anorm, test2)
        push!(log, :rnorm, test1)

        # The following tests guard against extremely small values of
        # atol, btol  or  ctol.  (The user may have set any or all of
        # the parameters  atol, btol, conlim  to 0.)
        # The effect is equivalent to the normal tests using
        # atol = eps,  btol = eps,  conlim = 1/eps.
        if itn >= maxiter  istop = 7; end
        if 1 + test3  <= 1 istop = 6; end
        if 1 + test2  <= 1 istop = 5; end
        if 1 + t1     <= 1 istop = 4; end

        # Allow for tolerances set by the user
        if  test3 <= ctol  istop = 3; end
        if  test2 <= atol  istop = 2; end
        if  test1 <= rtol  istop = 1; end
    end
    shrink!(log)
    setconv(log, istop > 0)
    x
end

function lsqr!(x, A, b; maxiter::Int=max(size(A,1), size(A,2)), kwargs...)
    T = Adivtype(A, b)
    z = zero(T)
    history = ConvergenceHistory()
    reserve!(history,[:resnorm,:anorm,:rnorm,:cnorm],maxiter)
    lsqr_method!(history, x, A, b; maxiter=maxiter, kwargs...)
    x, history
end

lsqr(A, b; kwargs...) = lsqr!(zerox(A, b), A, b; kwargs...)


text:
using ComparativeAutograder

tests = Array{TestCase, 1}()

for _ in 1:10
    push!(
        tests,
        TestCase((rand(), ))
    )
end

testsuite = TestSuite(
    "square",
    tests,
)

runtestsuite(testsuite)


text:
@testset "Examples" begin
    @testset "Naive sin" begin
        function sin_naive!(res, x)
            s, t, u = zero(x), zero(x), zero(x)
            tol = one(x)
            Arblib.mul_2exp!(tol, tol, -precision(tol))
            # @show tol
            k = 0
            while true
                Arblib.pow!(t, x, UInt(2k + 1))
                Arblib.fac!(u, UInt(2k + 1))
                Arblib.div!(t, t, u)
                Arblib.abs!(u, t)

                if u ≤ tol
                    # @show u
                    Arblib.add_error!(s, u)
                    break
                end
                if iseven(k)
                    Arblib.add!(s, s, t)
                else
                    Arblib.sub!(s, s, t)
                end
                k += 1
            end
            Arblib.set!(res, s)
            return res
        end

        @test let
            str = ""
            prec = 64
            while true
                x = Arb("2016.1"; prec = prec)
                y = zero(x)
                y = sin_naive!(y, x)
                str *= sprint(print, "Using $(lpad(prec, 5)) bits, sin(x) = ")
                str *= sprint(println, Arblib.string_nice(y, 10))
                y < zero(y) && break
                prec *= 2
            end
            str
        end == """
        Using    64 bits, sin(x) = [+/- 2.67e+859]
        Using   128 bits, sin(x) = [+/- 1.30e+840]
        Using   256 bits, sin(x) = [+/- 3.60e+801]
        Using   512 bits, sin(x) = [+/- 3.01e+724]
        Using  1024 bits, sin(x) = [+/- 2.18e+570]
        Using  2048 bits, sin(x) = [+/- 1.22e+262]
        Using  4096 bits, sin(x) = [-0.7190842207 +/- 1.20e-11]
        """
    end

    @testset "Logistic" begin
        """
            logistic(n[; x0 = 0.5, r = 3.75, digits = 10])
        Compute the nth iterate of the logistic map `x_{n + 1} = r x_n (1 -
        x_n).
        """
        function logistic(n; x0 = Arb(0.5), r = Arb(3.75), digits = 10)
            goalprec = digits * 3.3219280948873623 + 3

            prec = 64
            if typeof(x0) == String
                x = Arb(x0, prec = prec)
            else
                x = Arb(prec = prec)
            end

            t = zero(x)
            while true
                if typeof(x0) == String
                    x = Arb(x0, prec = prec)
                else
                    Arblib.set!(x, x0)
                end

                for i = 1:n
                    Arblib.sub!(t, x, 1)
                    Arblib.neg!(t, t)
                    Arblib.mul!(x, x, t)
                    Arblib.mul!(x, x, r)

                    if Arblib.rel_accuracy_bits(x) < goalprec
                        break
                    end
                end

                if Arblib.rel_accuracy_bits(x) >= goalprec
                    break
                end

                prec *= 2
                x = setprecision(x, prec)
                r = setprecision(r, prec)
                t = setprecision(t, prec)
            end

            return x
        end

        @test string(logistic(10, digits = 10)) == "[0.645367290830930 +/- 3.45e-16]"
        @test string(logistic(10, digits = 20)) ==
              "[0.64536729083093027156146131423635101 +/- 5.18e-36]"
    end
end


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11
@doc """
  SIGMOID Compute sigmoid functoon
    J = SIGMOID(z) computes the sigmoid of z.
""" ->
function sigmoid(z)
  g = SIGMOID(z)
  return g
end


text:
"""
    to_linear_operator(lind::GraphLindbladian, [T=eltype(lind)])

Converts the Lindbladian to a `KLocalOperatorSum` with eltype `T` for the non
Hermitian Hamiltonian and a set of KLocalOperator(Sum) for the dissipators.

This structure has minimal memory cost, and is NEEDED when the need to represent
systems in very big lattices arises. Evidently, there is a small performance
price to pay when using those structures.
"""
to_linear_operator(lind::GraphLindbladian, T=nothing) =
    to_linear_operator(hamiltonian(lind), jump_operators_graph(lind), T)

function to_linear_operator(ham::GraphOperator, c_ops::Vector, T::Union{Nothing, Type{<:Number}}=nothing)

    ham_locs = ham.LocalOperators

    # default type
    T = isnothing(T) ?  eltype(first(ham_locs).data) : T

    op_loc = KLocalOperatorRow(T, [1], [length(basis(first(ham_locs)))],
                               first(ham_locs).data)

    H = KLocalOperatorSum(op_loc)

    # Add the local hamiltonian terms
    for (i, h_loc) = enumerate(ham.LocalOperators)
        i == 1 && continue

        dim = length(basis(h_loc))
        sum!(H, KLocalOperatorRow(T, [i], [dim], h_loc.data))
    end

    # Hoppings
    for edge=edges(graph(ham))
        connections = ham.EdgeOpList[edge]
        for (coeff, l_id, r_id)=connections
            hl = ham.HopOperatorsList[l_id]
            hr = ham.HopOperatorsList[r_id]

            dim_l = length(basis(hl))
            dim_r = length(basis(hr))
            sum!(H, KLocalOperatorRow(T, [edge.src, edge.dst],
                                      [dim_l, dim_r],
                                      coeff*((hl⊗hr).data)))
        end
    end

    # store the vector of local-acting, individual dissipators
    loss_ops   = []
    loss_ops_t = []

    # dissipations
    # TODO support nonlocal dissipators
    for (i, L) = enumerate(c_ops)

        # Find local basis where the operator is non-zero
        nz_sites = Int[]
        for (j, L_j) = enumerate(L.LocalOperators)
            length(L_j.data.nzval) == 0 && continue
            push!(nz_sites, j)
        end

        isempty(nz_sites) && continue

        hilb_dims = length.(basis.(L.LocalOperators[nz_sites]))

        L_nz      = tensor(L.LocalOperators[nz_sites]...)

        op        = KLocalOperatorRow(T, nz_sites, hilb_dims, L_nz.data)

        op_hnh    = KLocalOperatorRow(T, nz_sites, hilb_dims,
                                      -im/2*(L_nz'*L_nz).data)

        sum!(H, op_hnh)
        push!(loss_ops, op)
        push!(loss_ops_t, transpose(op))
    end

    T = isempty(loss_ops) ? Nothing : typeof(first(loss_ops))

    loss_ops   = Vector{T}(loss_ops)
    loss_ops_t = Vector{T}(loss_ops_t)

    return (H, loss_ops, loss_ops_t)
end

"""
    to_linear_operator(lind::GraphLindbladian, [T=eltype(lind)])

Converts the Hamiltonian to a `KLocalOperatorSum` with eltype `T`.

This structure has minimal memory cost, and is NEEDED when the need to represent
systems in very big lattices arises. Evidently, there is a small performance
price to pay when using those structures.
"""
function to_linear_operator(op::GraphOperator, T::Union{Nothing, Type{<:Number}}=nothing)

    op_locs = op.LocalOperators

    # default type
    T = isnothing(T) ?  eltype(first(op_locs).data) : T

    op_loc = KLocalOperatorRow(T, [1], [length(basis(first(op_locs)))],
                        first(op_locs).data)

    res_op = KLocalOperatorSum(op_loc)

    # Add the local hamiltonian terms
    for (i, op_loc) = enumerate(op.LocalOperators)
        i == 1 && continue # the first one was already added by the constructor.

        dim = length(basis(op_loc))
        sum!(res_op, KLocalOperatorRow(T, [i], [dim], op_loc.data))
    end

    # Hoppings
    for edge=edges(graph(op))
        edge ∉ keys(op.EdgeOpList) && continue

        connections = op.EdgeOpList[edge]
        for (coeff, l_id, r_id)=connections
            hl = op.HopOperatorsList[l_id]
            hr = op.HopOperatorsList[r_id]

            dim_l = length(basis(hl))
            dim_r = length(basis(hr))
            sum!(res_op, KLocalOperatorRow(T, [edge.src, edge.dst],
                                           [dim_l, dim_r],
                                           coeff*((hl⊗hr).data)))
        end
    end

    return res_op
end


text:
using SequGen
using Unitful, Test

@testset "Constant" begin
  t = 1:10
  for c in (5, 1.23, u"1s", u"1K")
    s = Constant(c)
    @test sample(s, t) == fill(c, 10)
  end
end

@testset "Trend" begin
  for Δt in (1, 1.0, 1.12u"s", 2.5u"hr"),
    t₀ in (zero(Δt), oneunit(Δt)),
    Δx in (3, 3.14, 3.4u"m", 7.5u"K")

    t = collect(range(t₀, t₀ + Δt, length = 101))
    s = Trend(t₀, Δx / Δt)
    x = collect(range(zero(Δx), Δx, length = 101))
    @test sample(s, t) ≈ x
  end
end


text:
"""
Struct for Lewis Integration Method

		bsProcess=LewisMethod(A,N)

Where:\n
		A	=	volatility of the process.
		N	=	volatility of the process.
"""
mutable struct LewisMethod{num <: Number, num_1 <: Integer} <: AbstractIntegralMethod
    A::num
    N::num_1
    function LewisMethod(A::num, N::num_1) where {num <: Number, num_1 <: Integer}
        @assert(A > 0.0, "A must be positive")
        @assert(N > 2, "N must be greater than 2")
        return new{num, num_1}(A, N)
    end
end

export LewisMethod;

"""
Documentation Lewis Method
"""
function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, method::LewisMethod, abstractPayoff::FinancialMonteCarlo.EuropeanOption, ::FinancialMonteCarlo.BaseMode = FinancialMonteCarlo.SerialMode())
    T = abstractPayoff.T
    K = abstractPayoff.K
    A = method.A
    N = method.N
    S0 = mcProcess.underlying.S0
    r = FinancialMonteCarlo.integral(zero_rate.r, T) / T
    d = FinancialMonteCarlo.integral(FinancialMonteCarlo.dividend(mcProcess), T) / T
    cf = FinancialFFT.CharactheristicFunction(mcProcess, T)
    corr = FinancialFFT.CharactheristicExponent(-1im, mcProcess, T)
    CharFunc(v) = cf(v) * exp(-v * 1im * corr)
    mod = log(S0 / K) + (r - d) * T
    func_(z) = real_mod(exp(-z * 1im * mod) * CharFunc(-z - 1im * 0.5) / (z^2 + 0.25))
    int_1 = midpoint_definite_integral(func_, -A, A, N)
    price = S0 * (1 - exp(-mod / 2) * int_1 / (2 * pi)) * exp(-d * T)
    return call_to_put(price, mcProcess.underlying, zero_rate, abstractPayoff)
end

function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, method::LewisMethod, abstractPayoffs::Array{U}, ::FinancialMonteCarlo.BaseMode = FinancialMonteCarlo.SerialMode()) where {U <: FinancialMonteCarlo.EuropeanOption}
    TT = unique([opt.T for opt in abstractPayoffs])
    zero_typed = FinancialMonteCarlo.predict_output_type_zero(mcProcess, zero_rate, abstractPayoffs)
    prices = Array{typeof(zero_typed)}(undef, length(abstractPayoffs))

    for T in TT
        index_same_t = findall(op -> (op.T == T), abstractPayoffs)
        payoffs = abstractPayoffs[index_same_t]
        prices[index_same_t] .= [pricer(mcProcess, zero_rate, method, payoff) for payoff in payoffs]
    end

    return prices
end


text:
module Remark

using Base64
using Dates
using Documenter.Utilities.DOM
using EnglishText: SemanticText
using FunctionalCollections: PersistentHashMap
using SchemeSyntax
using SExpressions.Parser
using SExpressions.Lists
using SExpressions.Keywords

include("RemarkStates.jl")
include("stdlib.jl")
using .RemarkStates

const ListOrArray = Union{List, Array}

function makeenv(ass=Dict(), modules=[])
    Env = Module(gensym(:Env))
    Core.eval(Env, quote
        using Base.Iterators
        using SExpressions.Lists
        using SExpressions.Keywords
        using SchemeSyntax
        using SchemeSyntax.RacketExtensions
        using SchemeSyntax.R5RS
        using Remarkable.Remark.StdLib
    end)
    for (k, v) in ass
        Core.eval(Env, :(const $k = $v))
    end
    for touse in modules
        Core.eval(Env, :(const $(nameof(touse)) = $touse))
        Core.eval(Env, :(using .$(nameof(touse))))
    end
    Env
end

setindex(x, y, z) = assoc(x, z, y)

"""
Return `true` if `α` is a list, and each element in `α` is a list.
"""
islisty(α::List) = all(β -> isa(β, List), α)
islisty(_) = false

"""
Monad-ish: thread second argument through successive calls of `f` to elements
of `xs`.
"""
function acc2(f, xs, acc)
    res = []
    for x in xs
        r, acc = f(x, acc)
        push!(res, r)
    end
    res, acc
end

function handleinclude(obj, kind::Keyword, state)
    if kind == Keyword("remark")
        url = evaluate!(state, obj)
        file = relativeto(state, url)
        α = Parser.parsefile(file)
        acc2(domify, α, state)
    else
        error("Unknown included object type $kind")
    end
end

function handleinclude(ρ, state)
    if length(ρ) ≠ 2
        error("include must take exactly two arguments")
    end
    handleinclude(car(ρ), cadr(ρ), state)
end

function handleremark(ρ, state)
    if isempty(ρ)
        error("remark requires a nonempty body expression")
    end
    StdLib.setstate!(state)
    domify(evaluateall!(state, ρ), state)
end

function handleremarks(ρ, state)
    if isempty(ρ)
        error("remarks requires a nonempty body expression")
    end
    StdLib.setstate!(state)
    result = evaluateall!(state, ρ)
    if result === nothing
        nothing, state
    else
        acc2(domify, result, state)
    end
end

flattentree(::Nothing) = DOM.Node[]
flattentree(xs::ListOrArray)::Vector{DOM.Node} =
    vcat((flattentree(x) for x in xs)...)
flattentree(x) = DOM.Node[x]

"""
Evaluate a single entry in an attribute list.

Supported syntax is:

 - [class "value"] ⇒ equivalent to `class="value"` in HTML
 - [disabled]      ⇒ equivalent to `disabled=""`, or just `disabled`, in HTML

The value, if specified, can be any expression.
"""
function attribute(β, state)
    if length(β) == 1
        car(β) => ""
    elseif length(β) == 2
        attr, expr = β
        attr => string(evaluate!(state, expr))
    else
        throw(ArgumentError("Expected [attr] or [attr \"value\"], got $β"))
    end
end

function domify(head, ρ, state)
    error("Invalid HTSX head: $head")
end

function domify(head::Symbol, ρ, state)
    if head == :include
        handleinclude(ρ, state)
    elseif head == :remark
        handleremark(ρ, state)
    elseif head == :remarks
        handleremarks(ρ, state)
    elseif isnil(ρ)
        DOM.Node(head, DOM.Attributes(undef, 0), DOM.Node[]), state
    else
        if islisty(car(ρ))  # is a list of attrs
            attrs = [attribute(β, state) for β in car(ρ)]
            content, state = acc2(domify, cdr(ρ), state)
        else  # is just another body element
            attrs = DOM.Attributes(undef, 0)
            content, state = acc2(domify, ρ, state)
        end
        children = flattentree(content)
        DOM.Node(head, attrs, collect(DOM.Node, children)), state
    end
end

quoted(x) = list(:quote, x)
function domify(head::Keyword, ρ, state)
    if head == Keyword("each")
        var, array, code = ρ
        doms = Core.eval(state.env, quote
            [$(tojulia(code)) for $var in $(tojulia(array))]
        end)
        objects = []
        for dom in doms
            res, state = acc2(domify, dom, state)
            push!(objects, res)
        end
        objects, state
    else
        error("Unsupported HTSX keyword $head")
    end
end

domify(::Nil, state) = error("Empty list not allowed here")
function domify(α::Cons, state)
    head = car(α)
    domify(head, cdr(α), state)
end

domify(s::String, state) = DOM.Node(s), state
domify(s::AbstractString, state) = domify(String(s), state)
domify(s::SemanticText, state) = domify(string(s), state)
domify(i::Number, state) = domify(string(i), state)
domify(::Nothing, state) = nothing, state

domify(x, state) = error("Can’t serialize $(repr(x))")

function show_html(io::IO, dom)
    join(io, (stringmime("text/html", p) for p in dom))
end

function tohtml(io::IO, α::List, tmpls=PersistentHashMap{Symbol,Any}();
                file=joinpath(pwd(), "_implicit.rem"),
                modules=[])
    println(io, "<!DOCTYPE html>")
    state = RemarkState(makeenv(tmpls, modules), file)
    dom, _ = acc2(domify, α, state)
    show_html(io, flattentree(dom))
end

function tohtml(io::IO, f::AbstractString,
                tmpls=PersistentHashMap{Symbol,Any}();
                modules=[])
    tohtml(io, Parser.parsefile(f), tmpls; file=abspath(f), modules=modules)
end

tohtml(α::Union{List,AbstractString}, tmpls=PersistentHashMap{Symbol,Any}()) =
    sprint(tohtml, α, tmpls)

end


text:
const DATASET_METADATA_FILES = ("_common_metadata", "_metadata")

"""
    Parquet.Dataset(path; kwargs...)

Returns the table contained in the parquet dataset in an Tables.jl compatible format.
A dataset comprises of multiple parquet files and optionally some metadata files.

These options if provided are passed along while reading each parquet file in the dataset:
- `filter`: Filter function that takes the path to partitioned file and returns boolean to indicate whether to include the partition while loading. All partitions are loaded by default.
- `batchsize`: Maximum number of rows to read in each batch (default: row count of first row group). Applied to each file in the partition.
- `use_threads`: Whether to use threads while reading the file; applicable only for Julia v1.3 and later and switched on by default if julia processes is started with multiple threads.
- `column_generator`: Function to generate a partitioned column when not found in the partitioned table. Parameters provided to the function: table, column index, length of column to generate. Default implementation determines column values from the table path.

One can easily convert the returned object to any Tables.jl compatible table e.g. DataFrames.DataFrame via

```
using DataFrames
df = DataFrame(read_parquet(path))
```
"""
struct Dataset <: Tables.AbstractColumns
    path::String
    filter::Function
    ncols::Int
    kwargs::NamedTuple{(:batchsize, :use_threads, :column_generator), Tuple{Union{Nothing,Signed}, Bool, Function}}
    schema::Tables.Schema
    lookup::Dict{Symbol, Int} # map column name => index
    columns::Vector{AbstractVector}
    tables::Vector{Table}

    function Dataset(path;
            filter::Function=(path)->true,
            batchsize::Union{Nothing,Signed}=nothing,
            column_generator::Function=column_generator,
            use_threads::Bool=(nthreads() > 1))

        isdir(path) || error("Invalid Dataset path. Not a directory - $path")
        sch = dataset_schema(string(path))
        ncols = length(sch.names)
        lookup = Dict{Symbol, Int}(nm => i for (i, nm) in enumerate(sch.names))
        kwargs = (batchsize=batchsize, use_threads=use_threads, column_generator=dataset_column_generator)
        new(path, filter, ncols, kwargs, sch, lookup, AbstractVector[], Table[])
    end
end

const PARTITION_DATE_FORMATS = [dateformat"Y-m-d", dateformat"Y-m-d HH:MI:SS", dateformat"Y-m-dTHH:MI:SS"]
const PARTITION_DATETIME_FORMATS = [dateformat"Y-m-d HH:MI:SS", dateformat"Y-m-dTHH:MI:SS", dateformat"Y-m-d"]
function parse_date(missingstrval)
    for format in PARTITION_DATE_FORMATS
        try
            return Date(missingstrval, format)
        catch ex
            (format == last(PARTITION_DATE_FORMATS)) && rethrow()
        end
    end
end
function parse_datetime(missingstrval)
    for format in PARTITION_DATETIME_FORMATS
        try
            return DateTime(missingstrval, format)
        catch ex
            (format == last(PARTITION_DATETIME_FORMATS)) && rethrow()
        end
    end
end

function dataset_column_generator(table::Table, colidx::Int, len::Int)
    table_path = getfield(table, :path)
    schema = getfield(table, :schema)
    colname = schema.names[colidx]
    coltype = schema.types[colidx]

    pattern = Regex("\\S+[\\/\\\\]?$(colname)=([a-zA-Z0-9 :\\-\\.]*)[\\/\\\\]?\\S+")
    matches = match(pattern, table_path)
    if (matches !== nothing) && (length(matches.captures) == 1)
        missingstrval = matches.captures[1]
        nm_coltype = nonmissingtype(coltype)
        if nm_coltype <: Real
            missingval = parse(nm_coltype, lowercase(missingstrval))
        elseif nm_coltype <: Date
            missingval = parse_date(missingstrval)
        elseif nm_coltype <: DateTime
            missingval = parse_datetime(missingstrval)
        elseif nm_coltype <: String
            missingval = string(missingstrval)
        else
            error("Unhandled dataset partitioned column type $nm_coltype for column $colname of table $table_path")
        end
    else
        missingval = nonmissingtype(coltype) === coltype ? undef : missing
    end

    if missingval === undef || missingval === missing
        Array{coltype}(missingval, len)
    else
        fill!(Array{coltype}(undef, len), missingval)
    end
end

function close(dataset::Dataset)
    empty!(getfield(dataset, :columns))
    for table in getfield(dataset, :tables)
        close(table)
    end
    empty!(getfield(dataset, :tables))
    nothing
end

function dataset_schema(path::String)
    schema = nothing

    # look for _common_metadata or _metadata file
    for name in DATASET_METADATA_FILES
        meta_file = joinpath(path, name)
        if isfile(meta_file)
            schema = tables_schema(Parquet.File(meta_file))
            break
        end
    end

    if schema === nothing
        # else extract schema from any file in the dataset
        for (root, dirs, files) in walkdir(path)
            for file in files
                full_filename = joinpath(root, file)
                if is_par_file(full_filename)
                    schema = tables_schema(Parquet.File(full_filename))
                    break
                end
            end
            (schema === nothing) || break
        end
    end

    schema
end

"""
Iterator to iterate over partitions of a parquet dataset, returned by the `Tables.partitions(dataset)` method.
Each partition is a Parquet.Table.
"""
struct DatasetPartitions
    dataset::Dataset
    ncols::Int
    filter::Function

    function DatasetPartitions(dataset::Dataset, filter::Function)
        new(dataset, getfield(dataset, :ncols), filter)
    end
end

function iterated_partition(partitions::DatasetPartitions, cursor)
    partition = nothing
    walker, root, files, fileidx, step = cursor
    schema = getfield(partitions.dataset, :schema)

    while partition === nothing
        if 0 < fileidx <= length(files)  # we are iterating on files in a directory
            file = files[fileidx]
            fileidx += 1
            if !(file in DATASET_METADATA_FILES)
                full_filename = joinpath(root, file)
                if partitions.filter(full_filename) && is_par_file(full_filename)
                    partition = Table(full_filename, schema; getfield(partitions.dataset, :kwargs)...)
                end
            end
        else # walk further into directory tree
            itervals = (fileidx == 0) ? iterate(walker) : iterate(walker, step)
            (itervals === nothing) && (return nothing)  # end of directory tree

            dirinfo, step = itervals
            root, dirs, files = dirinfo
            fileidx = 1
        end
    end
    partition, (walker, root, files, fileidx, step)
end
Base.iterate(partitions::DatasetPartitions, cursor) = iterated_partition(partitions, cursor)
Base.iterate(partitions::DatasetPartitions) = iterated_partition(partitions, (walkdir(getfield(partitions.dataset, :path)), "", [], 0, nothing))

loaded(dataset::Dataset) = !isempty(getfield(dataset, :columns))
function load(dataset::Dataset)
    tables = getfield(dataset, :tables)
    columns = getfield(dataset, :columns)
    ncols = getfield(dataset, :ncols)
    empty!(tables)
    empty!(columns)
    for table in Tables.partitions(dataset)
        push!(tables, table)
        if isempty(columns)
            for colidx in 1:ncols
                push!(columns, ChainedVector([Tables.getcolumn(table, colidx)]))
            end
        else
            for colidx in 1:ncols
                append!(columns[colidx], Tables.getcolumn(table, colidx))
            end
        end
    end
    nothing
end

Tables.istable(::Dataset) = true
Tables.columnaccess(::Dataset) = true
Tables.schema(d::Dataset) = getfield(d, :schema)
Tables.columnnames(d::Dataset) = getfield(d, :schema).names
Tables.columns(d::Dataset) = Tables.CopiedColumns(d)
Tables.getcolumn(d::Dataset, nm::Symbol) = Tables.getcolumn(d, getfield(d, :lookup)[nm])
function Tables.getcolumn(d::Dataset, i::Int)
    loaded(d) || load(d)
    getfield(d, :columns)[i]
end
Tables.partitions(d::Dataset) = DatasetPartitions(d, getfield(d, :filter))