text:
export maximal_order, pmaximal_order,poverorder, MaximalOrder, ring_of_integers

###############################################################################
#
#  Maximal Order interface
#
###############################################################################
@doc Markdown.doc"""
    maximal_order(O::NfAbsOrd; index_divisors::Vector{fmpz}, discriminant::fmpz, ramified_primes::Vector{fmpz}) -> NfAbsOrd

Returns the maximal order of the number field that contains $O$. Additional information can be supplied if they are already known, as the ramified primes,
the discriminant of the maximal order or a set of integers dividing the index of $O$ in the maximal order.
"""
function MaximalOrder(O::NfAbsOrd{S, T}; index_divisors::Vector{fmpz} = fmpz[], discriminant::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[]) where {S, T}
  K = nf(O)
  try
    # First check if the number field knows its maximal order
    M = _get_maximal_order(K)::typeof(O)
    return M
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    M = new_maximal_order(O, index_divisors = index_divisors, disc = discriminant, ramified_primes = ramified_primes)
    M.ismaximal = 1
    _set_maximal_order(K, M)
    return M
  end
end

@doc Markdown.doc"""
    maximal_order(K::Union{AnticNumberField, NfAbsNS}; discriminant::fmpz, ramified_primes::Vector{fmpz}) -> NfAbsOrd

Returns the maximal order of $K$. Additional information can be supplied if they are already known, as the ramified primes
or the discriminant of the maximal order.

# Example

```julia-repl
julia> Qx, xx = FlintQQ["x"];
julia> K, a = NumberField(x^3 + 2, "a");
julia> O = MaximalOrder(K);
```
"""
function MaximalOrder(K::AnticNumberField; discriminant::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[])
  try
    c = _get_maximal_order(K)::NfAbsOrd{AnticNumberField, nf_elem}
    return c
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    #O = MaximalOrder(K)::NfOrd
    O = new_maximal_order(any_order(K))
    O.ismaximal = 1
    _set_maximal_order(K, O)
    return O
  end
end

@doc Markdown.doc"""
    ring_of_integers(K::AnticNumberField) -> NfAbsOrd

This function returns the ring of integers of $K$.
"""
function ring_of_integers(x::T; kw...) where T
  return maximal_order(x; kw...)
end

function maximal_order(f::T) where T <: Union{fmpz_poly, fmpq_poly}
  K = number_field(f, cached = false)[1]
  return maximal_order(K)
end

################################################################################
#
#  function to get an order which is maximal at some primes
#
################################################################################
@doc Markdown.doc"""
    pmaximal_overorder_at(O::NfOrd, primes::Array{fmpz, 1}) - > NfOrd

Given a set of prime numbers, this function returns an overorder of $O$ which
is maximal at those primes.
"""
function pmaximal_overorder_at(O::NfOrd, primes::Array{fmpz, 1})
  if isempty(primes)
    return O
  end

  OO = O

  if !isdefining_polynomial_nice(nf(O)) || !isinteger(gen_index(O))
    for i in 1:length(primes)
      p = primes[i]
      @vprint :NfOrd 1 "Computing p-maximal overorder for $p ..."
      OO += pmaximal_overorder(OO, p)
      if !(p in OO.primesofmaximality)
        push!(OO.primesofmaximality, p)
      end
    end
    return OO
  end

  ind = index(O)
  EO = EquationOrder(nf(O))
  M = zero_matrix(FlintZZ, 2 * degree(O), degree(O))
  for i in 1:length(primes)
    p = primes[i]
    @vprint :NfOrd 1 "Computing p-maximal overorder for $p ..."
    if divisible(ind, p)
      OO = sum_as_Z_modules(OO, pmaximal_overorder(O, p), M)
    else
      O1 = pmaximal_overorder(EO, p)
      if divisible(index(O1), p)
        OO = sum_as_Z_modules(OO, O1, M)
      end
    end
    if !(p in OO.primesofmaximality)
      push!(OO.primesofmaximality, p)
    end
    @vprint :NfOrd 1 "done\n"
  end
  return OO
end
################################################################################
#
#  Buchmann Lenstra heuristic
#
################################################################################

#  Buchmann-Lenstra for simple absolute number fields.
function new_maximal_order(O::NfOrd; index_divisors::Vector{fmpz} = fmpz[], disc::fmpz = fmpz(-1), ramified_primes::Vector{fmpz} = fmpz[])

  K = nf(O)
  if degree(K) == 1
    O.ismaximal = 1
    return O
  end

  if isdefining_polynomial_nice(K) && (isequation_order(O) || isinteger(gen_index(O)))
    Zx, x = PolynomialRing(FlintZZ, "x", cached = false)
    f1 = Zx(K.pol)
    ds = rres(f1, derivative(f1))
  else
    ds = discriminant(O)
  end

  #First, factorization of the discriminant given by the snf of the trace matrix
  M = trace_matrix(O)
  l = divisors(M, ds)
  if !isempty(index_divisors)
    append!(l, index_divisors)
  end
  if !isempty(ramified_primes)
    append!(l, ramified_primes)
  end
  if disc != -1
    push!(l, disc)
  end
  l = coprime_base(l)
  @vprint :NfOrd 1 "Factors of the discriminant: $l\n "
  l1 = fmpz[]
  OO = O
  @vprint :NfOrd 1 "Trial division of the discriminant\n "
  auxmat = zero_matrix(FlintZZ, 2*degree(K), degree(K))
  first = true
  for d in l
    if disc != -1
      u = divexact(discriminant(OO), disc)
      if isone(gcd(u, d))
        continue
      end
    end
    fac = factor_trial_range(d)[1]
    rem = d
    for (p,v) in fac
      rem = divexact(rem, p^v)
    end
    @vprint :NfOrd 1 "Computing the maximal order at $(collect(keys(fac)))\n "
    O1 = pmaximal_overorder_at(O, collect(keys(fac)))
    if first
      OO = O1
      first = false
    else
      @vtime :NfOrd 3 OO = sum_as_Z_modules(OO, O1, auxmat)
    end
    rem = abs(rem)
    if !isone(rem)
      if disc != -1
        u = divexact(discriminant(OO), disc)
        if isone(gcd(u, rem))
          continue
        end
      end
      push!(l1, rem)
    end
  end
  if isempty(l1) || discriminant(OO) == disc
    OO.ismaximal = 1
    return OO
  end
  for i=1:length(l1)
    a, b = ispower(l1[i])
    if a>1
      if isprime(b)
        O1 = pmaximal_overorder(O, b)
        OO = sum_as_Z_modules(OO, O1, auxmat)
        l1[i] = 0
      else
        l1[i]=b
      end
    end
  end
  ll1 = fmpz[x for x in l1 if !iszero(x)]
  if isempty(ll1)
    OO.ismaximal = 1
    return OO
  end
  O1, Q = _TameOverorderBL(OO, ll1)
  if !isempty(Q) && discriminant(O1) != disc
    @vprint :NfOrd 1 "I have to factor $Q\n "
    for el in Q
      d = factor(el).fac
      O1 = pmaximal_overorder_at(O, collect(keys(d)))
      OO = sum_as_Z_modules(OO, O1, auxmat)
    end
  end
  O1.ismaximal = 1
  return O1

end

function _TameOverorderBL(O::NfOrd, lp::Array{fmpz,1})

  OO = O
  M = coprime_base(lp)
  Q = fmpz[]
  while !isempty(M)
    @vprint :NfOrd 1 M
    q = pop!(M)
    if isprime(q)
      OO1 = pmaximal_overorder(O, q)
      if valuation(discriminant(OO1), q) < valuation(discriminant(OO), q)
        OO += OO1
      end
    else
      OO, q1 = _cycleBL(OO, q)
      if q1 == q
        push!(Q, q)
      elseif !isone(q1)
        push!(M, q1)
        push!(M, divexact(q, q1))
        M = coprime_base(M)
      end
    end
  end
  if isempty(Q)
    OO.ismaximal = 1
  end
  return OO, Q

end

function _qradical(O::NfOrd, q::fmpz)

  d = degree(O)
  R = ResidueRing(FlintZZ, q, cached=false)
  #First, we compute the q-radical as the kernel of the trace matrix mod q.
  #By theory, this is free if q is prime; if I get a non free module, I have found a factor of q.
  @vprint :NfOrd 1 "radical computation\n "
  rk, M = kernel(trace_matrix(O), R, side = :left)
  if iszero(rk)
    @vprint :NfOrd 1 "The radical is equal to the ideal generated by q"
    return fmpz(1), ideal(O, q)
  end
  M = howell_form(M)
  for i = 1:nrows(M)
    if iszero_row(M, i)
      break
    end
    j = i
    while iszero(M[i,j])
      j += 1
    end
    if !isone(M[i,j])
      @vprint :NfOrd 1 "Split: $(M[i,j])"
      return lift(M[i,j]), NfOrdIdl()
    end
  end
  # Now, we have the radical.
  # We want to compute the ring of multipliers.
  # So we need to lift the matrix.
  @vprint :NfOrd 1 "Computing hnf of basis matrix \n "
  MatIdeal = zero_matrix(FlintZZ, d, d)
  for i=1:nrows(M)
    for j=1:degree(O)
      MatIdeal[i, j] = M[i,j].data
    end
  end
  gens = NfOrdElem[O(q)]
  for i=1:nrows(M)
    if !iszero_row(MatIdeal, i)
      push!(gens, elem_from_mat_row(O, MatIdeal, i))
    end
  end
  M2 = _hnf_modular_eldiv(MatIdeal, fmpz(q),  :lowerleft)
  I = NfOrdIdl(O, M2)
  I.gens = gens
  return fmpz(1), I
end

function _cycleBL(O::NfOrd, q::fmpz)

  q1, I = _qradical(O, q)
  if !isone(q1)
    return O, q1
  elseif isdefined(I, :princ_gens) && q == I.princ_gens
    return O, fmpz(1)
  end
  @vprint :NfOrd 1 "ring of multipliers\n"
  O1 = ring_of_multipliers(I)
  @vprint :NfOrd 1 "ring of multipliers computed\n"
  while discriminant(O1) != discriminant(O)
    if isone(gcd(discriminant(O1), q))
      return O1, fmpz(1)
    end
    O = O1
    q1, I = _qradical(O, q)
    if !isone(q1)
      return O, q1
    elseif isdefined(I, :princ_gens) && q == I.princ_gens
      return O, fmpz(1)
    end
    @vprint :NfOrd 1 "ring of multipliers\n"
    O1 = ring_of_multipliers(I)
  end
  @vprint :NfOrd 1 "couldn't enlarge\n"
  # (I:I)=OO. Now, we want to prove tameness (or get a factor of q)
  # We have to test that (OO:a)/B is a free Z/qZ module.
  inva = colon(ideal(O, 1), I, true)
  M1 = basis_mat_inv(inva)
  @assert isone(M1.den)
  G1 = divisors(M1.num, q)
  for i = 1:length(G1)
    q1 = gcd(q, G1[i])
    if q1 != q && !isone(q1)
      @vprint :NfOrd 1 "Found the factor $q1"
      return O, q1
    end
  end
  @vprint :NfOrd 1 "...is free\n"
  return _cycleBL2(O, q, I)

end

function _cycleBL2(O::NfOrd, q::fmpz, I::NfOrdIdl)

  h = 2
  ideals = Array{NfOrdIdl,1}(undef, 3)
  ideals[1] = I
  ideals[2] = I*I
  ideals[3] = ideals[2] * I
  while true
    if h > degree(O)
      error("Not found!")
    end
    I1 = (ideals[1] + ideal(O, q))*(ideals[3] + ideal(O, q))
    I2 = (ideals[2] + ideal(O, q))^2
    M2 = basis_mat(I2, copy = false)*basis_mat_inv(I1, copy = false)
    @assert isone(M2.den)
    G2 = divisors(M2.num, q)
    if isempty(G2)
      h += 1
      ideals[1] = ideals[2]
      ideals[2] = ideals[3]
      ideals[3] = ideals[2]*I
      continue
    end
    for i = 1:length(G2)
      q1 = gcd(q, G2[i])
      if q1 != q && !isone(q1)
        return O, q1
      end
    end
    break
  end
  f, r = ispower(q, h)
  if f
    return O, r
  else
    return O, q
  end
end

function TameOverorderBL(O::NfOrd, lp::Array{fmpz,1}=fmpz[])

  # First, we hope that we can get a factorization of the discriminant by computing
  # the structure of the group OO^*/OO
  OO=O
  list = append!(elementary_divisors(trace_matrix(OO)), primes_up_to(degree(O)))
  l=coprime_base(list)
  #Some trivial things, maybe useless
  for i=1:length(l)
    a,b=ispower(l[i])
    if a>1
      l[i]=b
    end
    if isprime(l[i])
      @vprint :NfOrd 1 "pmaximal order at $(l[i])\n"
      OO1=pmaximal_overorder(O, l[i])
      if valuation(discriminant(OO1), l[i])<valuation(discriminant(OO), l[i])
        OO+=OO1
      end
      l[i]=0
    end
  end
  push!(l, discriminant(OO))
  append!(l,lp)
  filter!(x-> !iszero(x), l)
  for i=1:length(l)
    l[i]=abs(l[i])
  end
  M=coprime_base(l)
  Q=fmpz[]
  while !isempty(M)
    @vprint :NfOrd 1 M
    q=M[1]
    if isprime(q)
      OO1=pmaximal_overorder(O, q)
      if valuation(discriminant(OO1), q)< valuation(discriminant(OO), q)
        OO+=OO1
      end
      filter!(x-> x!=q, M)
      continue
    end
    OO, q1 = _cycleBL(OO,q)
    if isone(q1)
      filter!(x->x!=q, M)
    elseif q1 == q
      push!(Q, q)
      filter!(x-> x != q, M)
    else
      push!(M, q1)
      push!(M, divexact(q,q1))
      M = coprime_base(M)
    end
  end
  if isempty(Q)
    OO.ismaximal=1
  end
  return OO, Q

end

################################################################################
#
#  p-overorder
#
################################################################################

function _poverorder(O::NfAbsOrd, p::fmpz)
  @vtime :NfOrd 3 I = pradical(O, p)
  if isdefined(I, :princ_gen) && I.princ_gen == p
    return O
  end
  @vtime :NfOrd 3 R = ring_of_multipliers(I)
  return R
end

@doc Markdown.doc"""
    poverorder(O::NfOrd, p::fmpz) -> NfOrd
    poverorder(O::NfOrd, p::Integer) -> NfOrd

This function tries to find an order that is locally larger than $\mathcal O$
at the prime $p$: If $p$ divides the index $[ \mathcal O_K : \mathcal O]$,
this function will return an order $R$ such that
$v_p([ \mathcal O_K : R]) < v_p([ \mathcal O_K : \mathcal O])$. Otherwise
$\mathcal O$ is returned.
"""
function poverorder(O::NfAbsOrd, p::fmpz)
  if isequation_order(O) && issimple(nf(O))
    #return dedekind_poverorder(O, p)
    return polygons_overorder(O, p)
  else
    return _poverorder(O, p)
  end
end

function poverorder(O::NfAbsOrd, p::Integer)
  return poverorder(O, fmpz(p))
end

################################################################################
#
#  p-maximal overorder
#
################################################################################

@doc Markdown.doc"""
    pmaximal_overorder(O::NfOrd, p::fmpz) -> NfOrd
    pmaximal_overorder(O::NfOrd, p::Integer) -> NfOrd

This function finds a $p$-maximal order $R$ containing $\mathcal O$. That is,
the index $[ \mathcal O_K : R]$ is not dividible by $p$.
"""
function pmaximal_overorder(O::NfAbsOrd, p::fmpz)
  @vprint :NfOrd 1 "computing p-maximal overorder for $p ... \n"
  if p in O.primesofmaximality
    return O
  end
  d = discriminant(O)
  if !iszero(rem(d, p^2))
    push!(O.primesofmaximality, p)
    return O
  end

  @vprint :NfOrd 1 "extending the order at $p for the first time ... \n"
  i = 1
  OO = poverorder(O, p)
  dd = discriminant(OO)
  while d != dd
    if !iszero(rem(dd, p^2))
      break
    end
    i += 1
    @vprint :NfOrd 1 "extending the order at $p for the $(i)th time ... \n"
    d = dd
    OO = poverorder(OO, p)
    dd = discriminant(OO)
  end
  push!(OO.primesofmaximality, p)
  return OO
end

function pmaximal_overorder(O::NfAbsOrd, p::Integer)
  return pmaximal_overorder(O, fmpz(p))
end

################################################################################
#
#  Ring of multipliers
#
################################################################################

@doc Markdown.doc"""
    ring_of_multipliers(I::NfAbsOrdIdl) -> NfAbsOrd
Computes the order $(I : I)$, which is the set of all $x \in K$
with $xI \subseteq I$.
"""
function ring_of_multipliers(a::NfAbsOrdIdl)
  O = order(a)
  n = degree(O)
  bmatinv = basis_mat_inv(a, copy = false)
  if isdefined(a, :gens) && length(a.gens) < n
    B = a.gens
  else
    B = basis(a, copy = false)
  end
  @assert length(B) > 0
  id_gen = zero_matrix(FlintZZ, 2*n, n)
  m = zero_matrix(FlintZZ, n*length(B), n)
  ind = 1
  for i = 1:length(B)
    if i != 1
      c = matrix(FlintZZ, 1, n, coordinates(B[i]))
      reduce_mod_hnf_ll!(c, id_gen)
      if iszero(c)
        continue
      end
    end
    M = representation_matrix_mod(B[i], minimum(a)*bmatinv.den)
    _copy_matrix_into_matrix(id_gen, 1, 1, M)
    hnf_modular_eldiv!(id_gen, minimum(a), :lowerleft)
    mod!(M, minimum(a)*bmatinv.den)
    mul!(M, M, bmatinv.num)
    M = transpose(M)
    _copy_matrix_into_matrix(m, n*(ind-1)+1, 1, M)
    if view(id_gen, n+1:2*n, 1:n) == basis_mat(a, copy = false)
      m = view(m, 1:n*ind, 1:n)
      break
    end
    ind += 1
  end
  if !isone(bmatinv.den)
    divexact!(m, m, bmatinv.den)
  end
  mhnf = hnf_modular_eldiv!(m, minimum(a))
  s = prod(mhnf[i,i] for i = 1:n)
  if isone(s)
    return O
  end
  # mhnf is upper right HNF
  mhnf = transpose(view(mhnf, 1:n, 1:n))
  b = FakeFmpqMat(pseudo_inv(mhnf))
  mul!(b, b, basis_mat(O, copy = false))
  @hassert :NfOrd 1 defines_order(nf(O), b)[1]
  O1 = NfAbsOrd(nf(O), b)
  if isdefined(O, :disc)
    O1.disc = divexact(O.disc, s^2)
  end
  if isdefined(O, :index)
    O1.index = s*O.index
    O1.gen_index = fmpq(O1.index)
  end
  if isdefined(O, :basis_mat_inv)
    O1.basis_mat_inv = O.basis_mat_inv * mhnf
  end
  return O1
end

################################################################################
#
#  p-radical
#
################################################################################

function pradical_trace(O::NfAbsOrd, p::Union{Integer, fmpz})
  d = degree(O)
  M = trace_matrix(O)
  F = GF(p, cached = false)
  M1 = change_base_ring(M, F)
  k, B = nullspace(M1)
  if iszero(k)
    return ideal(O, p)
  end
  M2 = zero_matrix(FlintZZ, d, d)
  for i = 1:ncols(B)
    for j = 1:d
      M2[i, j] = FlintZZ(B[j, i].data)
    end
  end
  gens = elem_type(O)[O(p)]
  for i=1:ncols(B)
    if !iszero_row(M2,i)
      push!(gens, elem_from_mat_row(O, M2, i))
    end
  end
  M2 = hnf_modular_eldiv!(M2, fmpz(p), :lowerleft)
  I = ideal(O, M2)
  I.minimum = p
  I.gens = gens
  return I
end

@doc Markdown.doc"""
    factor_shape_refined(f::gfp_poly)

Given a polynomial f over a finite field, it returns an array having one
entry for every irreducible factor giving its degree and its multiplicity.
"""
function factor_shape_refined(x::gfp_poly) where {T <: RingElem}
  res = Tuple{Int, Int}[]
  square_fac = factor_squarefree(x)
  for (f, i) in square_fac
    discdeg = factor_distinct_deg(f)
    for (j, g) in discdeg
      num = divexact(degree(g), j)
      for l = 1:num
        push!(res, (j, i))
      end
    end
  end
  return res
end

function pradical_frobenius1(O::NfOrd, p::Union{Integer, fmpz})
  R = GF(p, cached = false)
  d = degree(O)
  K = nf(O)
  Rx = PolynomialRing(R, "x", cached = false)[1]
  res = factor_shape_refined(Rx(K.pol))
  md = 1
  for i = 1:length(res)
    md = max(md, res[i][2])
  end
  j = clog(fmpz(md), p)
  sqf = factor_squarefree(Rx(K.pol))
  p1 = one(Rx)
  for (x, v) in sqf
    if v > 1
      p1 = mul!(p1, p1, x)
    end
  end
  gen2 = O(lift(K, p1))
  M1 = representation_matrix_mod(gen2, fmpz(p))
  hnf_modular_eldiv!(M1, fmpz(p), :lowerleft)
  nr = 0
  indices = Int[]
  for i = 1:d
    if !isone(M1[i, i])
      push!(indices, i)
      nr += 1
    end
  end
  A = zero_matrix(R, d, nr + d)
  B = basis(O, copy = false)
  ind = 0
  for i in 1:d
    if !(i in indices)
      continue
    end
    t = powermod(B[i], p^j, p)
    ind += 1
    if iszero(t)
      continue
    end
    ar = coordinates(t, copy = false)
    for k in 1:d
      A[k, ind] = ar[k]
    end
  end
  for s = 1:d
    for i = 1:s
      A[i, s+nr] = R(M1[s, i])
    end
  end
  X = right_kernel_basis(A)
  gens = elem_type(O)[O(p), gen2]
  if isempty(X)
    I = ideal(O, p)
    I.gens = gens
    return I
  end
  #First, find the generators
  for i = 1:length(X)
    coords = zeros(FlintZZ, d)
    for j=1:nr
      coords[indices[j]] = lift(X[i][j])
    end
    if !iszero(coords)
      push!(gens, O(coords))
    end
  end
  #Then, construct the basis matrix of the ideal
  m = zero_matrix(FlintZZ, nr + d, d)
  for i = 1:length(X)
    for j = 1:nr
      m[i, indices[j]] = lift(X[i][j])
    end
  end
  for i = 1:d
    for s = 1:d
      m[i+nr, s] = M1[i, s]
    end
  end
  hnf_modular_eldiv!(m, fmpz(p), :lowerleft)
  I = NfAbsOrdIdl(O, view(m, nr+1:nr+d, 1:d))
  I.minimum = p
  I.gens = gens
  return I
end

function pradical_frobenius(O::NfAbsOrd, p::Union{Integer, fmpz})

  #First, I need an exponent for the maximum of the nilpotency indices.
  R = GF(p, cached = false)
  d = degree(O)
  K = nf(O)
  if issimple(K) && isdefining_polynomial_nice(K) && contains_equation_order(O)
    return pradical_frobenius1(O, p)
  end
  j = clog(fmpz(d), p)
  @assert p^(j-1) < d
  @assert d <= p^j

  A = zero_matrix(R, degree(O), degree(O))
  B = basis(O, copy = false)
  for i in 1:d
    t = powermod(B[i], p^j, p)
    if iszero(t)
      continue
    end
    ar = coordinates(t)
    for k in 1:d
      A[k, i] = ar[k]
    end
  end
  X = right_kernel_basis(A)
  gens = elem_type(O)[O(p)]
  if isempty(X)
    I = ideal(O, p)
    I.gens = gens
    return I
  end
  #First, find the generators
  for i = 1:length(X)
    coords = Array{fmpz,1}(undef, d)
    for j=1:d
      coords[j] = lift(X[i][j])
    end
    push!(gens, O(coords))
  end
  #Then, construct the basis matrix of the ideal
  m = zero_matrix(FlintZZ, d, d)
  for i = 1:length(X)
    for j = 1:d
      m[i, j] = lift(X[i][j])
    end
  end
  mm = hnf_modular_eldiv!(m, fmpz(p), :lowerleft)
  I = NfAbsOrdIdl(O, mm)
  I.minimum = p
  I.gens = gens
  return I

end

@doc Markdown.doc"""
    pradical(O::NfOrd, p::{fmpz|Integer}) -> NfAbsOrdIdl

Given a prime number $p$, this function returns the $p$-radical
$\sqrt{p\mathcal O}$ of $\mathcal O$, which is
just $\{ x \in \mathcal O \mid \exists k \in \mathbf Z_{\geq 0} \colon x^k
\in p\mathcal O \}$. It is not checked that $p$ is prime.
"""
function pradical(O::NfAbsOrd, p::Union{Integer, fmpz})
  if p isa fmpz
    if nbits(p) < 64
      return pradical(O, Int(p))
    end
  end
  d = degree(O)

  #Trace method if the prime is large enough
  if p > d
    return pradical_trace(O, p)
  else
    res = pradical_frobenius(O, p)
    return res
  end
end


text:
correct = ( pi/4 + atan(3.0) ) / 2
for endpt in [neither, both, left, right]
    x, w = legendre(40, endpt)
    Q = sum(w .* 1 ./ ( 1 .+ (2x .- 1).^2))
    @test Q ≈ correct
end

ext: jl


text:
#
#Function to calculate the distance between two points considering PBC
#
function distpbc(p1,p2,side)
           dx = p1[1] - p2[1]
            x = mod(dx,side)
             if x > side/2
              x = x - side
             end
             if x < -side/2
              x = x + side
             end
           dy = p1[2] - p2[2]
            y = mod(dy,side)
             if y > side/2
              y = y - side
             end
             if y < -side/2
              y = y + side
             end
             d12 = sqrt(x^2 + y^2)
           end


text:
function MDP_fixture()
    dmodel = NoCrashIDMMOBILModel(3, PhysicalParam(3, lane_length=100.0))
    rmodel = NoCrashRewardModel()
    return NoCrashMDP{typeof(rmodel)}(dmodel, rmodel, 1.0)
end
function state_fixture()
    # XXX are we going to git rid of v_min, etc?
    lane_length=100.0
    pp = PhysicalParam(3, lane_length=lane_length)
    return MLState(0.0,0.0, CarState[CarState(lane_length/2.0, 2.0, pp.v_med, 0.0, Multilane.NORMAL,0),
                                CarState(lane_length*3/4, 1.0, pp.v_min, 0.0,
                                         IDMMOBILBehavior("normal", pp.v_min, 0.0, 1),0),
                                CarState(lane_length*1/4, 3.0, pp.v_max, 0.0,
                                         IDMMOBILBehavior("normal", pp.v_max, 0.0, 2),0)])
end

function test_occupation()
    println("\t\tTestng occupation_overlap")
    @test Multilane.occupation_overlap(1.0, 2.0) == false
    @test Multilane.occupation_overlap(1.0, 1.99) == true
    @test Multilane.occupation_overlap(0.9, 2.0) == false
    @test Multilane.occupation_overlap(0.9, 1.1) == true
    @test Multilane.occupation_overlap(0.9, 1.99) == true
end

# not done
# function test_e_brake()
#     mdp = MDP_fixture()
#     pp = mdp.dmodel.phys_param
#     s = MLState(false,0.0,0.0, CarState[CarState(0.0, 2.0, pp.v_med, 0.0, Multilane.NORMAL),
#                                 CarState(pp.l_car+0.1, 2.0, pp.v_min, 0.0,
#                                          BehaviorModel("normal"), pp.v_min, 0.0, 1)])
# end

"""
Test that cars are regenerated when there are some missing
"""
function test_encounter()
    println("\t\tTesting car regeneration for encounters")
    mdp = MDP_fixture()
    mdp.dmodel.p_appear = 1.0
    s = state_fixture()
    deleteat!(s.cars, 2)
    rng = MersenneTwister(2)
    sp, r = generate_sr(mdp, s, MLAction(0.0,0.0), rng)
    @test length(sp.cars) == 3
end

"""
Return a CarState for the ego vehicle with the specified y.
"""
ego_state_spec_y(y) = CarState(0.0, y, 30.0, 0.0, Multilane.NORMAL,1)

function test_snapback()
    println("\t\tTesting snapping back after lane change too far")
    mdp = MDP_fixture()
    rng = MersenneTwister(2)
    mdp.dmodel.phys_param.dt = 0.75
    s = state_fixture()
    s.cars[1] = ego_state_spec_y(2.0)
    a1 = MLAction(0,1.0)
    sp,r = generate_sr(mdp, s, a1, rng)
    @test sp.cars[1].y == 2.0 + 1.0*0.75
    a2 = MLAction(0,2.0)
    sp,r = generate_sr(mdp, s, a2, rng)
    @test sp.cars[1].y == 3.0
    s.cars[1] = ego_state_spec_y(1.99)
    a3 = MLAction(0,1.0)
    sp,r = generate_sr(mdp, s, a3, rng)
    @test sp.cars[1].y == 2.0
    s.cars[1] = ego_state_spec_y(2.01)
    a4 = MLAction(0,-1.0)
    sp,r = generate_sr(mdp, s, a4, rng)
    @test sp.cars[1].y == 2.0
    s.cars[1] = ego_state_spec_y(1.01)
    a5 = MLAction(0.,0.1)
    sp,r = generate_sr(mdp, s, a5, rng)
    @test sp.cars[1].y == 1.01 + 0.1*0.75
end

function test_model()
    println("\tTesting Model...")
    test_occupation()
    test_encounter()
    test_snapback()
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
push!(LOAD_PATH, "/Users/thomasriddick/Documents/workspace/Dynamic_HD_Code/Dynamic_HD_Julia_Code/")
using Profile
using ArgParse
using HDDriverModule: drive_hd_model,drive_hd_and_lake_model
using FieldModule: UnstructuredField,Field,repeat
using IOModule: load_river_parameters, load_river_initial_values,get_ncells
using IOModule: load_lake_initial_values,load_lake_parameters
using IOModule: get_additional_grid_information
using GridModule: UnstructuredGrid

function pass_arguments()
  settings = ArgParseSettings()
  @add_arg_table settings begin
  "--hd-para-file", "-p"
    help = "Filepath to HD parameters file"
    required = true
  "--hd-init-file", "-i"
    help = "Filepath to HD initialisation file (optional)"
  "--lake-para-file", "-l"
    help = "Filepath to lake model parameters file (optional)"
  "--lake-init-file", "-n"
    help = "Filepath to lake model initialisation file (optional)"
  "--drainage-file", "-d"
    help = "Filepath to file containing drainage values"
  "--runoff-file", "-r"
    help = "Filepath to file containing runoff values"
  "--lake-evaporation-file", "-e"
    help = "Filepath to file containing values for evaporation from lakes"
  "--timesteps", "-t"
    help = "Number of timesteps to run"
    arg_type = Int
    default=100
  "--step-length","-s"
    help = "Length of each timestep in seconds"
    arg_type = Float64
    default = 86400.0
  end
  return parse_args(settings)
end

function main()
  args = pass_arguments()
  timesteps = args["timesteps"]
  step_length = args["step-length"]
  ncells = get_ncells(args["hd-para-file"])
  local clat::Array{Float64,1}
  local clon::Array{Float64,1}
  local clat_bounds::Array{Float64,2}
  local clon_bounds::Array{Float64,2}
  clat,clon,clat_bounds,clon_bounds = get_additional_grid_information(args["hd-para-file"])
  grid = UnstructuredGrid(ncells,clat,clon,clat_bounds,clon_bounds)
  lake_grid = UnstructuredGrid(ncells,clat,clon,clat_bounds,clon_bounds)
  surface_model_grid = UnstructuredGrid(ncells,clat,clon,clat_bounds,clon_bounds)
  river_parameters = load_river_parameters(args["hd-para-file"],grid;
                                           day_length=86400.0,step_length=step_length)
  local drainages::Array{Field{Float64},1}
  local runoffs::Array{Field{Float64},1}
  local lake_evaporations::Array{Field{Float64},1}
  if args["drainage-file"] != nothing
    drainages = load_drainage_fields(args["drainage-file"],grid,
                                     last_timestep=12)
    drainages = [divide(x,30.0) for x in drainages]
  else
    drainage::Field{Float64} = UnstructuredField{Float64}(river_parameters.grid,
                                                          100*0.0000000227*step_length*2.6*10000000000)
    #drainages = repeat(drainage,Int(round(timesteps/30)+1))
    drainages = [drainage]
  end
  if args["runoff-file"] != nothing
    runoffs = load_runoff_fields(args["runoff-file"],grid,
                                 last_timestep=12)
    runoffs = [divide(x,30.0) for x in runoffs]
  else
    runoffs = deepcopy(drainages)
  end
  if args["lake-evaporation-file"] != nothing
    lake_evaporations = load_lake_evaporation_fields(args["lake-evaporation-file"],grid,
                                                     last_timestep=12)
    lake_evaporations = [divide(x,30.0) for x in lake_evaporations]
  else
    lake_evaporation_none ::Field{Float64} = UnstructuredField{Float64}(river_parameters.grid,0.0)
    lake_evaporation_high ::Field{Float64} =
      UnstructuredField{Float64}(river_parameters.grid,
                                 100000*0.0000000227*step_length*2.6*10000000000)
    lake_evaporations_none = repeat(lake_evaporation_none,Int(round(timesteps/60)+1))
    lake_evaporations_high = repeat(lake_evaporation_high,Int(round(timesteps/60)+1))
    lake_evaporations = vcat(lake_evaporations_none,lake_evaporations_high)
  end
  if args["lake-para-file"] != nothing
    lake_parameters = load_lake_parameters(args["lake-para-file"],lake_grid,grid,surface_model_grid)
    drainages_copy = deepcopy(drainages)
    runoffs_copy = deepcopy(runoffs)
    if args["hd-init-file"] != nothing
      river_fields = load_river_initial_values(args["hd-init-file"],grid,river_parameters)
    end
    if args["lake-init-file"] != nothing
      initial_water_to_lake_centers::UnstructuredField{Float64},
      initial_spillover_to_rivers::UnstructuredField{Float64} =
        load_lake_initial_values(args["lake-init-file"],lake_grid,grid)
      if args["hd-init-file"] != nothing
        drive_hd_and_lake_model(river_parameters,river_fields,
                                lake_parameters,drainages,runoffs,
                                lake_evaporations,timesteps,true,
                                initial_water_to_lake_centers,
                                initial_spillover_to_rivers;print_timestep_results=false)
        # Profile.clear()
        # Profile.init(delay=0.01)
        # @time drive_hd_and_lake_model(river_parameters,river_fields,
        #                               lake_parameters,drainages_copy,runoffs_copy,
        #                               timesteps,true,initial_water_to_lake_centers,
        #                               initial_spillover_to_rivers;print_timestep_results=true)
        # Profile.print()
        # r = Profile.retrieve();
        # f = open("/Users/thomasriddick/Downloads/profile.bin", "w")
        # Serialization.serialize(f, r)
        # close(f)
      else
        drive_hd_and_lake_model(river_parameters,lake_parameters,
                                drainages,runoffs,
                                lake_evaporations,timesteps,true,
                                initial_water_to_lake_centers,
                                initial_spillover_to_rivers;
                                print_timestep_results=true)
      end
    end
  else
    if args["hd-init-file"] != nothing
      river_fields = load_river_initial_values(args["hd-init-file"],grid,river_parameters)
      drive_hd_model(river_parameters,river_fields,
                     drainages,runoffs,timesteps;
                     print_timestep_results=false,
                     output_timestep=160)
      # @time drive_hd_model(river_parameters,river_fields,
      #                      drainages_copy,runoffs_copy,timesteps;
      #                      print_timestep_results=false)
    else
      # drainages_copy = deepcopy(drainages)
      # runoffs_copy = deepcopy(runoffs)
      drive_hd_model(river_parameters,drainages,
                     runoffs,timesteps;print_timestep_results=false)
      # Profile.clear()
      # Profile.init(delay=0.0001)
      # @time drive_hd_model(river_parameters,drainages_copy,
      #                      runoffs_copy,timesteps;print_timestep_results=false)
      # Profile.print()
      # r = Profile.retrieve();
      # f = open("/Users/thomasriddick/Downloads/profile.bin", "w")
      # Serialization.serialize(f, r)
      # close(f)
    end
  end
end

empty!(ARGS)
# push!(ARGS,"-p/Users/thomasriddick/Documents/data/ICONHDdata/hdparafiles/hd_para_icon_r2b9_30_20191452_retuned_v3.nc")
# push!(ARGS,"-i/Users/thomasriddick/Documents/data/ICONHDdata/hdstartfiles/hdrestart_R02B09_015_G_241019_1337_v2.nc")
# push!(ARGS,"-t1920")
# push!(ARGS,"-s45")
push!(ARGS,"-p/Users/thomasriddick/Documents/data/temp/icon_lake_model_test/hdpara_icon.nc")
push!(ARGS,"-i/Users/thomasriddick/Documents/data/temp/icon_lake_model_test/hdrestart_R02B04_013_G_231019_1242_v2.nc")
push!(ARGS,"-l/Users/thomasriddick/Documents/data/temp/icon_lake_model_test/lakeparams.nc")
push!(ARGS,"-n/Users/thomasriddick/Documents/data/temp/icon_lake_model_test/lakestart.nc")
push!(ARGS,"-t3600")
push!(ARGS,"-s86400")
main()

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
mutable struct Der <: Component
	x
	dy
end
function print(io::IO,d::Der)
	if isa(d.dy,X) && isa(d.x,X)
		print(io,"∂$(d.x)/∂$(d.dy)")
	else
		print(io,"∂($(d.x))/∂($(d.dy))")
	end
end
getargs(d::Der)=[d.x,d.dy]
mutable struct DerOp <: NonAbelian
	dy
end
#*(dop::DerOp,ex::EX)=Der(ex,dop.dy)
function print(io::IO,dop::DerOp)
	if isa(dop.dy,X)
		print(io,"∂/∂$(dop.dy)")
	else
		print(io,"∂/∂($(dop.dy))")
	end
end
relations["Der"]=simplify(Equation[Der(:a,:x)≖0,Der(Oneable(:a)*:x,:x)≖:a,Der(Oneable(:a)*Pow(:x,:n),:x)≖:n*:a*Pow(:x,:n-1),Der(Oneable(:a)*Sqrt(:x),:x)≖0.5*:a*Pow(:x,-0.5)]) #Der(Pow(:a,:x),:x)≖log(:a)*Pow(:a,:x)
function matches(d::Der,pat::Der)
	mdx=matches(d.x,pat.x)
	mddy=matches(d.dy,pat.dy)
	validated=Dict[]
	for md in mdx
		for mdy in mddy
			md,mdy=simplify(md),simplify(mdy)
			con=false
			for key in keys(md)
				if key!=d.dy&&has(md[key],d.dy)
					con=true
				end
			end
			if con;continue;end
			if !clash(md,mdy)
				push!(validated,combine(md,mdy))
			end
		end
	end
	return validfilter(d,pat,validated)
end
matches(::N,::Der)=[]
matches(::Term,::Der)=[]
function matches(d::Der)
	ap=addparse(d.x)
	nap=Term[]
	for term in ap
		push!(nap,Factor[Der(Expression(term),d.dy)])
	end
	return Expression(nap)
end
function simplify(d::Der)
	return Der(simplify(d.x),simplify(d.dy))
end


text:
import StatsBase

"""
	mean(uv::AbstractUncertainValue, n::Int)

Compute the mean of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.mean(x::AbstractUncertainValue, n::Int) = resample(StatsBase.mean, x, n)

"""
	mode(uv::AbstractUncertainValue, n::Int)

Compute the mode of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.mode(x::AbstractUncertainValue, n::Int) = resample(StatsBase.mode, x, n)

"""
	quantile(uv::AbstractUncertainValue, q, n::Int)

Compute the quantile(s) `q` of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.quantile(x::AbstractUncertainValue, q, n::Int) = resample(StatsBase.quantile, x, n, q)

"""
	iqr(uv::AbstractUncertainValue, n::Int)

Compute the interquartile range (IQR), i.e. the 75th percentile minus the 25th percentile,
over an `n`-draw sample of an uncertain value.
"""
StatsBase.iqr(x::AbstractUncertainValue, n::Int) = resample(StatsBase.iqr, x, n, q)

"""
	median(uv::AbstractUncertainValue, n::Int)

Compute the median of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.median(x::AbstractUncertainValue, n::Int) = resample(StatsBase.median, x, n)

"""
	middle(uv::AbstractUncertainValue, n::Int)

Compute the middle of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.middle(x::AbstractUncertainValue, n::Int) = resample(StatsBase.middle, x, n)

"""
	std(uv::AbstractUncertainValue, n::Int)

Compute the standard deviation of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.std(x::AbstractUncertainValue, n::Int) = resample(StatsBase.std, x, n)

"""
	variance(uv::AbstractUncertainValue, n::Int)

Compute the variance of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.var(x::AbstractUncertainValue, n::Int) = resample(StatsBase.var, x, n)

"""
	genmean(uv::AbstractUncertainValue, p, n::Int)

Compute the generalized/power mean with exponent `p` of an uncertain value over an
`n`-draw sample of it.
"""
StatsBase.genmean(x::AbstractUncertainValue, p, n::Int) = resample(StatsBase.genmean, x, n, p)

"""
	genvar(uv::AbstractUncertainValue, n::Int)

Compute the generalized sample variance of an uncertain value over an
`n`-draw sample of it.
"""
StatsBase.genvar(x::AbstractUncertainValue, n::Int) = resample(StatsBase.genvar, x, n)

"""
	harmmean(uv::AbstractUncertainValue, n::Int)

Compute the harmonic mean of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.harmmean(x::AbstractUncertainValue, n::Int) = resample(StatsBase.harmmean, x, n)

"""
	geomean(uv::AbstractUncertainValue, n::Int)

Compute the geometric mean of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.geomean(x::AbstractUncertainValue, n::Int) = resample(StatsBase.geomean, x, n)

"""
	kurtosis(uv::AbstractUncertainValue, n::Int, m = mean(uv, n))

Compute the excess kurtosis of an uncertain value over an `n`-draw sample of it,
optionally specifying a center `m`).
"""
StatsBase.kurtosis(x::AbstractUncertainValue, n::Int, m = mean(x, n)) = resample(StatsBase.kurtosis, x, n, m)

"""
	moment(x::AbstractUncertainValue, k, n::Int, m = mean(x, n))

Compute the `k`-th order central moment of an uncertain value over an
`n`-draw sample of it, optionally specifying a center `m`.
"""
StatsBase.moment(x::AbstractUncertainValue, k, n::Int, m = mean(x, n)) = resample(StatsBase.moment, x, n, k, m)

"""
	percentile(x::AbstractUncertainValue, p, n::Int)

Compute the percentile(s) `p` of an uncertain value over an `n`-draw sample of it.
"""
StatsBase.percentile(x::AbstractUncertainValue, p, n::Int) = resample(StatsBase.percentile, x, n, p)

"""
	renyientropy(uv::AbstractUncertainValue, α, n::Int)

Compute the Rényi (generalized) entropy of order `α` of an uncertain value over an
`n`-draw sample of it.
"""
StatsBase.renyientropy(x::AbstractUncertainValue, α, n::Int) = resample(StatsBase.quantile, x, n, α)

"""
	rle(x::AbstractUncertainValue, n::Int)

Compute the run-length encoding of an uncertain value over a `n`-draw
sample of it as a tuple. The first element of the tuple is a vector of
values of the input and the second is the number of consecutive occurrences of
each element.
"""
StatsBase.rle(x::AbstractUncertainValue, n::Int) = resample(StatsBase.rle, x, n)

"""
	sem(x::AbstractUncertainValue, n::Int)

Compute the standard error of the mean of an uncertain value over a `n`-draw
sample of it, optionally specifying a center `m`, i.e.
`sqrt(var(x_draw, corrected = true) / length(x_draw))`.
"""
StatsBase.sem(x::AbstractUncertainValue, n::Int) = resample(StatsBase.sem, x, n)

"""
	skewness(x::AbstractUncertainValue, n::Int, m = mean(x, n))

Compute the standardized skewness of an uncertain value over an `n`-draw sample of
it, optionally specifying a center `m`.
"""
StatsBase.skewness(x::AbstractUncertainValue, n::Int; m = mean(x, n)) = resample(StatsBase.skewness, x, n, m)

"""
	span(x::AbstractUncertainValue, n::Int)

Compute the span of a collection, i.e. the range `minimum(x):maximum(x)`, of an
uncertain value over an `n`-draw sample of it.  The minimum and
maximum of the draws of `x` are computed in one pass using extrema.
"""
StatsBase.span(x::AbstractUncertainValue, n::Int) = resample(StatsBase.span, x, n)

"""
	summarystats(uv::AbstractUncertainValue, n::Int)

Compute summary statistics of an uncertain value over an `n`-draw sample of it. Returns
a `SummaryStats` object containing the mean, minimum, 25th percentile, median,
75th percentile, and maximum.
"""
StatsBase.summarystats(x::AbstractUncertainValue, n::Int) = resample(StatsBase.summarystats, x, n)

"""
	totalvar(uv::AbstractUncertainValue, n::Int)

Compute the total sample variance of an uncertain value over an
`n`-draw sample of it. For a single uncertain value, this is
equivalent to the sample variance.
"""
StatsBase.totalvar(x::AbstractUncertainValue, n::Int) = resample(StatsBase.totalvar, x, n)


text:
let
    ID_LIST = [ActionID(),ObjectID(),BotID{DeliveryBot}(),BotID{CleanUpBot}()]
    for id in ID_LIST
        @test read_abstract_id(TOML.parse(id)) == id
    end
end


text:
struct Closure{F, Free}
    frees :: Free
end

function (closure::Closure{F, X})(args...; kwargs...) where {F, X}
    F(closure.frees..., args...; kwargs...)
end


text:
#################### NaV ###############################

struct Na{F,D<:Real} <: FlowChannel(Sodium)
    gNa::D
    mNa::F
    hNa::F
end

Na(x) = Na(x, 0.0, 0.0)
m∞(::Na, V) = 1.0 / (1.0 + exp((V + 25.5) / -5.29))
h∞(::Na, V) = 1.0 / (1.0 + exp((V + 48.9) / 5.18))
τm(::Na, V) = 1.32 - 1.26 / (1 + exp((V + 120.0) / -25.0))
τh(::Na, V) = (0.67 / (1.0 + exp((V + 62.9) / -10.0))) * (1.5 + 1.0 / (1.0 + exp((V + 34.9) / 3.6)))

function (ch::Na)(n::Neuron; name=get_name(ch))

    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)
    @variables mNa(t) hNa(t)

    states, params = vardivide(V, mNa, hNa, I, g, E)
    eqs = [D(mNa) ~ (1 / τm(ch, V)) * (m∞(ch, V) - mNa),
        D(hNa) ~ (1 / τh(ch, V)) * (h∞(ch, V) - hNa),
        I ~ g * mNa^3 * hNa * (E - V)]

    defaultmap = [mNa => m∞(ch, n.somatic_parameters[Voltage]), hNa => h∞(ch, n.somatic_parameters[Voltage]), g => ch.gNa]
    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

#################### Slow calcium current #############################
#doesn't have to sense calcium. only ECa needs ca sensed but
struct CaS{F,D<:Real} <: FlowChannel(Calcium)
    gCaS::D
    mCaS::F
    hCaS::F
end

CaS(x) = CaS(x, 20.0, 0.0)
m∞(::CaS, V) = 1.0 / (1.0 + exp((V + 33.0) / -8.1))
h∞(::CaS, V) = 1.0 / (1.0 + exp((V + 60.0) / 6.2))
τm(::CaS, V) = 1.4 + 7.0 / (exp((V + 27.0) / 10.0) + exp((V + 70.0) / -13.0))
τh(::CaS, V) = 60.0 + 150.0 / (exp((V + 55.0) / 9.0) + exp((V + 65.0) / -16.0))

function (ch::CaS)(n::Neuron; name=get_name(ch))
    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)
    @variables mCaS(t) hCaS(t) V(t)

    eqs = [D(mCaS) ~ (1 / τm(ch, V)) * (m∞(ch, V) - mCaS),
        D(hCaS) ~ (1 / τh(ch, V)) * (h∞(ch, V) - hCaS),
        I ~ g * mCaS^3 * hCaS * (E - V)]

    states, params = vardivide(V, E, mCaS, hCaS, I, g)
    defaultmap = [mCaS => m∞(ch, n.somatic_parameters[Voltage]), hCaS => h∞(ch, n.somatic_parameters[Voltage]), g => ch.gCaS]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

#################### Transient calcium current ######################

struct CaT{F,D<:Real} <: FlowChannel(Calcium)
    gCaT::D
    mCaT::F
    hCaT::F
end

CaT(x) = CaT(x, 0.0, 0.0)
m∞(::CaT, V) = 1.0 / (1.0 + exp((V + 27.1) / -7.2))
h∞(::CaT, V) = 1.0 / (1.0 + exp((V + 32.1) / 5.5))
τm(::CaT, V) = 21.7 - 21.3 / (1.0 + exp((V + 68.1) / -20.5));
τh(::CaT, V) = 105.0 - 89.8 / (1.0 + exp((V + 55.0) / -16.9));

function (ch::CaT)(n::Neuron; name=get_name(ch))
    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)

    @variables mCaT(t) hCaT(t)
    eqs = [D(mCaT) ~ (1 / τm(ch, V)) * (m∞(ch, V) - mCaT),
        D(hCaT) ~ (1 / τh(ch, V)) * (h∞(ch, V) - hCaT),
        I ~ g * mCaT^3 * hCaT * (E - V)]

    states, params = vardivide(V, E, mCaT, hCaT, I, g)
    defaultmap = [mCaT => m∞(ch, n.somatic_parameters[Voltage]), hCaT => h∞(ch, n.somatic_parameters[Voltage]), g => ch.gCaT]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

#################### A-type potassium current #########################

struct Ka{F,D<:Real} <: FlowChannel(Potassium)
    gKa::D
    mKa::F
    hKa::F
end

Ka(x) = Ka(x, 0.0, 0.0)
m∞(::Ka, V) = 1.0 / (1.0 + exp((V + 27.2) / -8.7))
h∞(::Ka, V) = 1.0 / (1.0 + exp((V + 56.9) / 4.9))
τm(::Ka, V) = 11.6 - 10.4 / (1.0 + exp((V + 32.9) / -15.2))
τh(::Ka, V) = 38.6 - 29.2 / (1.0 + exp((V + 38.9) / -26.5))

function (ch::Ka)(n::Neuron; name=get_name(ch))
    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)

    @variables mKa(t) hKa(t)
    eqs = [D(mKa) ~ (1 / τm(ch, V)) * (m∞(ch, V) - mKa),
        D(hKa) ~ (1 / τh(ch, V)) * (h∞(ch, V) - hKa),
        I ~ g * mKa^3 * hKa * (E - V)]

    states, params = vardivide(V, E, mKa, hKa, I, g)
    defaultmap = [mKa => m∞(ch, n.somatic_parameters[Voltage]), hKa => h∞(ch, n.somatic_parameters[Voltage]), g => ch.gKa]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

################### Calcium-activated potassium current ########

struct KCa{F,D<:Real} <: FlowChannel(Calcium, Potassium)
    gKCa::D
    mKCa::F
end

KCa(x) = KCa(x, 0.0)
m∞(::KCa, V, Ca) = (Ca / (Ca + 3.0)) / (1.0 + exp((V + 28.3) / -12.6));
τm(::KCa, V) = 90.3 - 75.1 / (1.0 + exp((V + 46.0) / -22.7));

function (ch::KCa)(n::Neuron; name=get_name(ch))
    I, Ca, E, V = instantiate_hooks(n, ch, currents, sensed_ions, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)

    @variables mKCa(t)
    eqs = [D(mKCa) ~ (1 / τm(ch, V)) * (m∞(ch, V, Ca) - mKCa),
        I ~ (g * mKCa^4) * (E - V)]

    states, params = vardivide(V, Ca, E, mKCa, I, g)
    defaultmap = [mKCa => m∞(ch, n.somatic_parameters[Voltage], n.somatic_parameters[Calcium]), g => ch.gKCa]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

#################### Delayed rectifier potassium current ######################
struct Kdr{F,D<:Real} <: FlowChannel(Potassium)
    gKdr::D
    mKdr::F
end

Kdr(x) = Kdr(x, 0.0)
m∞(::Kdr, V) = 1.0 / (1.0 + exp((V + 12.3) / -11.8));
τm(::Kdr, V) = 7.2 - 6.4 / (1.0 + exp((V + 28.3) / -19.2));

function (ch::Kdr)(n::Neuron; name=get_name(ch))
    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)
    @variables mKdr(t)

    eqs = [D(mKdr) ~ (1 / τm(ch, V)) * (m∞(ch, V) - mKdr),
        I ~ (g * mKdr^4) * (E - V)]

    states, params = vardivide(V, E, mKdr, I, g)
    defaultmap = [mKdr => m∞(ch, n.somatic_parameters[Voltage]), g => ch.gKdr]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

#################### H current ####################

struct H{F,D<:Real} <: FlowChannel(Proton)
    gH::D
    mH::F
end

H(x) = H(x, 0.0)
m∞(::H, V) = 1.0 / (1.0 + exp((V + 70.0) / 6.0))
τm(::H, V) = (272.0 + 1499.0 / (1.0 + exp((V + 42.2) / -8.73)))

function (ch::H)(n::Neuron; name=get_name(ch))
    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)
    @variables mH(t)

    eqs = [D(mH) ~ (1 / τm(ch, V)) * (m∞(ch, V) - mH),
        I ~ g * mH * (E - V)]

    states, params = vardivide(V, E, mH, I, g)
    defaultmap = [mH => m∞(ch, n.somatic_parameters[Voltage]), g => ch.gH]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end

#################### Leak current #########################

struct leak{D<:Real} <: FlowChannel(Leak)
    gLeak::D
end

function (ch::leak)(n::Neuron; name=get_name(ch))
    I, E, V = instantiate_hooks(n, ch, currents, reversals, voltage)
    g, = instantiate_parameters(ch, conductances)

    eqs = [I ~ g * (E - V)]

    states, params = vardivide(V, E, I, g)
    defaultmap = [g => ch.gLeak]

    return ODESystem(eqs, t, states, params; defaults=defaultmap, name=name)
end


text:
module BinaryBandit

# precision constants: The rule of thumb is this: if you aggregate 10^N values, you will lose N digits of precision (http://www.ilikebigbits.com/2017_06_01_float_or_double.html); in our case N=3
# to be used relatively, i.e. if ( a - b ) > BB_numerical_precision_XX * min( abs( a ) + abs( b ) , floatmax( FloatXX ) (see also https://www.floating-point-gui.de/errors/comparison/)
const BB_numerical_precision_64 = 1.0e-13 # instead of 16
const BB_numerical_precision_32 = 1.0e-4 # instead of 7

include( "BB_DP_2.jl" )
include( "BB_evaluation_DP_2.jl" )

end # module

max_stars_repo_name: PeterJacko/BinaryBandit

max_issues_repo_name: PeterJacko/BinaryBandit

max_forks_repo_name: PeterJacko/BinaryBandit


text:
"""
    SpectralNorm(layer; ps=Flux.params(layer), power_iterations::Int=1)

Applies spectral normalization to `ps` parameters in the given layer.

!!!note:
This assumes that `layer` cannot contain non-trainable states.

Ref: [`https://arxiv.org/abs/1802.05957`](@ref)
"""
struct SpectralNorm{Re,P,S,PUV}
    layer_re::Re
    parameters::P
    states::S
    parameter_uvs::PUV
    power_iterations::Int
end

Flux.@functor SpectralNorm (parameters, states)

function Base.show(io::IO, sn::SpectralNorm)
    ps = Flux.params(sn)
    p = update_parameters(sn)
    l = wn.layer_re(p, sn.states)
    print(io, "SpectralNorm(")
    print(io, l)
    return print(") ", string(sum(length.(ps))), " Trainable Parameters")
end

l2normalize(x::AbstractArray{T}) where {T} = x ./ (norm(x, 2) + eps(T))

function construct_uv(x::AbstractArray)
    H, W = Flux.flatten(x)
    u = l2normalize(randn(eltype(x), H))
    v = l2normalize(randn(eltype(x), W))
    return (u, v)
end

function construct_uv(x::CuArray)
    H, W = Flux.flatten(x)
    u = l2normalize(CUDA.randn(eltype(x), H))
    v = l2normalize(CUDA.randn(eltype(x), W))
    return (u, v)
end

function spectral_normalization(x, u, v, power_iterations::Int=1)
    x_mat = Flux.flatten(x)
    for _ in 1:power_iterations
        v .= l2normalize(x_mat' * u)
        u .= l2normalize(x_mat * v)
    end
    Σ = sum(u .* (x_mat * v))
    return Σ
end

Zygote.@nograd spectral_normalization

function update_parameters(sn::SpectralNorm)
    return vcat(ntuple(i -> vec(sn.parameters[i] ./
                                spectral_normalization(sn.parameters[i], sn.parameter_uvs[i][1], sn.parameter_uvs[i][2],
                                                       sn.power_iterations)), length(sn.dims))...)
end

function SpectralNorm(layer; ps=Flux.params(layer), power_iterations::Int=1)
    if ps isa Array
        ps = Flux.params(ps)
    end

    _, s_, layer_re = destructure_parameters_states(layer, ps.params)

    parameters = []
    parameter_uvs = []
    foreach(p -> begin
                push!(parameter_uvs, construct_uv(p))
                push!(parameter, p)
            end, ps)

    return SpectralNorm(layer_re, tuple(parameters...), s_, tuple(parameter_uvs...), power_iterations)
end

(sn::SpectralNorm)(args...; kwargs...) = sn.layer_re(update_parameters(sn), sn.states)(args...; kwargs...)


text:
module Common

include("bytes.jl")
include("errors.jl")
include("Utils.jl")
include("BufferParser.jl")
include("CDObject.jl")
include("LibCrypto.jl")

end


text:
@testset "Vector: $TRef" for (TVec, T, TRef) in [
    (ArbVector, Arb, Arb),
    (AcbVector, Acb, Acb),
    (ArbRefVector, Arb, ArbRef),
    (AcbRefVector, Acb, AcbRef),
]
    V = TVec(4, prec = 128)
    @test size(V) == (4,)
    @test precision(V) == 128

    x = T(1.5)
    V[3] = x
    @test V[3] == x
    @test V[3] isa TRef
    @test !isempty(sprint(show, V[3]))
    @test precision(V[3]) == 128

    Arblib.add!(V, V, V)
    V2 = TVec(4, prec = 128)
    V2[3] = T(3)
    @test V == V2

    A = TVec([T(i + 1) for i = 2:5])
    @test A[end] == T(6)

    A = TVec([i + j for i = 1:4 for j = 1:4]; prec = 96)
    @test A[16] == T(8)
    @test precision(A) == 96
    @test precision(A[14]) == 96

    @test ref(A, 16) isa Union{ArbRef,AcbRef}

    # arithmetic
    AInt = [1, 2, 3, 4]
    BInt = [5, 6, 7, 8]
    A = TVec(AInt; prec = 96)
    B = TVec(BInt; prec = 96)
    @test A - B == AInt - BInt
    @test precision(A - B) == 96
    @test (A - B) isa TVec
    @test -B == -BInt
    @test -B isa TVec
    @test -B + A == A - B
    @test precision(-B + A) == 96

    C = TVec(length(A); prec = 96)
    Arblib.add!(C, A, B)
    @test C == AInt + BInt
    # test original signature
    D = TVec(length(A); prec = 96)
    Arblib.add!(D, A, B, length(AInt), 96)
    @test D == AInt + BInt

    @testset "similar" begin
        A = TVec(5; prec = 96)

        a = similar(A)
        @test a isa TVec
        @test precision(a) == precision(A)

        a = similar(A, TRef, 3)
        @test a isa TVec
        @test precision(a) == precision(A)

        for (ElT, VT, MT) in (
            (Arb, ArbVector, ArbMatrix),
            (Acb, AcbVector, AcbMatrix),
            (ArbRef, ArbRefVector, ArbRefMatrix),
            (AcbRef, AcbRefVector, AcbRefMatrix),
        )
            a = similar(A, ElT, 3)
            @test a isa VT
            @test precision(a) == precision(A)

            a = similar(A, ElT, (3, 2))
            @test a isa MT
            @test precision(a) == precision(A)
        end
    end

    @testset "set Rational" begin
        A = TVec(2)
        A[1, 1] = 5 // 8
        @test !iszero(A[1, 1])
    end
end

@testset "VectorRef: $T" for (T, TRef) in
                             [(ArbVector, ArbRefVector), (AcbVector, AcbRefVector)]
    A = T(5; prec = 96)
    A[4] = 3

    B = TRef(A)
    @test B isa TRef
    @test precision(B) == 96
    B[4] = 4
    @test A[4] == 4
    @test B[4] == 4
    C = T(B)
    @test C == A
    @test C[4] == 4
end


text:
n = 4
print("Testing NCSubSeq for ", n, " items:\n   ")
for a in NCSubSeq(n)
    print(" ", a)
end
println()

s = "Rosetta"
cs = split(s, "")
m = 10
n = length(NCSubSeq(length(s))) - m
println()
println("The first and last ", m, " NC sub-sequences of \"", s, "\":")
for (i,a) in enumerate(NCSubSeq(length(cs)))
    i <= m || n < i || continue
    println(@sprintf "%6d %s" i join(cs[a], ""))
    i == m || continue
    println("    .. ......")
end

t = {}
append!(t, collect(1:10))
append!(t, collect(20:10:40))
append!(t, big(50):50:200)
println()
println("Numbers of NC sub-sequences of a given length:")
for i in t
    println(@sprintf("%7d => ", i), length(NCSubSeq(i)))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
const bars = "-------------------------------------------------------------------------------"

"""

Print some information about the run.

"""
function title(R::Result,
               solute::Selection,
               solvent::Selection,
               nspawn::Int)

  println(bars)
  println("Starting MDDF calculation in parallel:")
  println("  $(R.nframes_read) frames will be considered")
  println("  Number of calculation threads: ",nspawn)
  println("  Solute: $(solute.natoms) atoms belonging to $(solute.nmols) molecules.")
  println("  Solvent: $(solvent.natoms) atom belonging to $(solvent.nmols) molecules.")

end

function title(R::Result,
               solute::Selection,
               solvent::Selection)

  println(bars)
  println("Starting MDDF calculation:")
  println("  $(R.nframes_read) frames will be considered")
  println("  Solute: $(solute.natoms) atoms belonging to $(solute.nmols) molecules.")
  println("  Solvent: $(solvent.natoms) atom belonging to $(solvent.nmols) molecules.")

end


text:
type View{A<:AbstractArray, B<:AbstractArray} <: Layer{A,B}
    output :: B
    gradInput :: A
end

function View(T::Type, input_size, output_size)
    @assert(prod(input_size) == prod(output_size),
            "View input / output size mismatch")
    return View(Array(T, output_size), Array(T, input_size))
end

function updateOutput!{A, B}(self :: View{A, B}, input :: A)
    self.output[:] = reshape(input, size(self.output))
    return self.output
end

function updateGradInput!{A, B}(self :: View{A, B},
                                input :: A,
                                gradOutput :: B)
    self.gradInput[:] = reshape(gradOutput, size(self.gradInput))
    return self.gradInput
end

export View


text:
module SFRDSP

using FFTW
using DSP
using Statistics
using LinearAlgebra
using SampledSignals

export stft2, istft2
export cubic
export xcorr2, circ_xcorr2, xc_align
export circ_smooth, circ_resample, circ_conv
export noise_psd
export fit_ir, IRFitException
export resample_filter2, resample2, pad

include("stft.jl")
include("interp.jl")
include("xcorr.jl")
include("circ.jl")
include("noise_psd.jl")
include("ir.jl")
include("resample.jl")

end # module

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
module Polytopes

using LightGraphs
using StaticArrays

include("polytope.jl")
include("show.jl")

export
        Polytope,
        Polygon,
        Polyhedron,
        rank,
        nverts,
        faces,
        facets,
        GeometricPolytope,
        Vertex,
        embeddim

end


text:
using gars

push!(ARGS, "arg")
gars.julia_main()


text:
mutable struct Bus
    bus_i::Int
    bustype::Int
    Pd::Float64
    Qd::Float64
    Gs::Float64
    Bs::Float64
    area::Int
    Vm::Float64
    Va::Float64
    baseKV::Float64
    zone::Int
    Vmax::Float64
    Vmin::Float64
end

mutable struct Line
    from::Int
    to::Int
    r::Float64
    x::Float64
    b::Float64
    rateA::Float64
    rateB::Float64
    rateC::Float64
    ratio::Float64 #TAP
    angle::Float64 #SHIFT
    status::Int
    angmin::Float64
    angmax::Float64
end
Line() = Line(0,0,0.,0.,0.,0.,0.,0.,0.,0.,0,0.,0.)

mutable struct Gener
    # .gen fields
    bus::Int
    Pg::Float64
    Qg::Float64
    Qmax::Float64
    Qmin::Float64
    Vg::Float64
    mBase::Float64
    status::Int
    Pmax::Float64
    Pmin::Float64
    Pc1::Float64
    Pc2::Float64
    Qc1min::Float64
    Qc1max::Float64
    Qc2min::Float64
    Qc2max::Float64
    ramp_agc::Float64
    # .gencost fields
    gentype::Int
    startup::Float64
    shutdown::Float64
    n::Int
    coeff::Array
end

mutable struct GenerVec{VI, VD}
    # .gen fields
    bus::VI
    Pg::VD
    Qg::VD
    Qmax::VD
    Qmin::VD
    Vg::VD
    mBase::VD
    status::VI
    Pmax::VD
    Pmin::VD
    Pc1::VD
    Pc2::VD
    Qc1min::VD
    Qc1max::VD
    Qc2min::VD
    Qc2max::VD
    ramp_agc::VD
    # .gencost fields
    gentype::VI
    startup::VD
    shutdown::VD
    n::VI
    #  coeff::Array
    coeff2::VD
    coeff1::VD
    coeff0::VD

    function GenerVec{VI, VD}(num_on) where {VI, VD}
        genvec = new{VI, VD}()
        genvec.bus = VI(undef, num_on)
        genvec.Pg = VD(undef, num_on)
        genvec.Qg = VD(undef, num_on)
        genvec.Qmax = VD(undef, num_on)
        genvec.Qmin = VD(undef, num_on)
        genvec.Vg = VD(undef, num_on)
        genvec.mBase = VD(undef, num_on)
        genvec.status = VI(undef, num_on)
        genvec.Pmax = VD(undef, num_on)
        genvec.Pmin = VD(undef, num_on)
        genvec.Pc1 = VD(undef, num_on)
        genvec.Pc2 = VD(undef, num_on)
        genvec.Qc1min = VD(undef, num_on)
        genvec.Qc1max = VD(undef, num_on)
        genvec.Qc2min = VD(undef, num_on)
        genvec.Qc2max = VD(undef, num_on)
        genvec.ramp_agc = VD(undef, num_on)
        genvec.gentype = VI(undef, num_on)
        genvec.startup = VD(undef, num_on)
        genvec.shutdown = VD(undef, num_on)
        genvec.n = VI(undef, num_on)
        genvec.coeff0 = VD(undef, num_on)
        genvec.coeff1 = VD(undef, num_on)
        genvec.coeff2 = VD(undef, num_on)
        return genvec
    end
end

struct OPFData
    buses::Array{Bus}
    lines::Array{Line}
    generators::Array{Gener}
    bus_ref::Int
    baseMVA::Float64
    BusIdx::Dict{Int,Int}    #map from bus ID to bus index
    FromLines::Array         #From lines for each bus (Array of Array)
    ToLines::Array           #To lines for each bus (Array of Array)
    BusGenerators::Array     #list of generators for each bus (Array of Array)
end

mutable struct Ybus{VD}
    YffR::VD
    YffI::VD
    YttR::VD
    YttI::VD
    YftR::VD
    YftI::VD
    YtfR::VD
    YtfI::VD
    YshR::Array{Float64}
    YshI::Array{Float64}

    Ybus{VD}(YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI) where {VD} = new(YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI)
end

function _sort_keys(data, k)
    bus2idx = Dict{Int,Int}()
    orig_index = Int[]
    for (_, e) in data
        push!(orig_index, e[k])
    end
    sorted_index = sort(orig_index)
    # Reverse mapping
    for (i, idx) in enumerate(sorted_index)
        bus2idx[idx] = i
    end
    return sorted_index, bus2idx
end

function _opf_loaddata_mat(case_name, lineOff=Line())
    # import data as Dict(...)
    baseMVA = 100
    data = PowerModels.parse_file(case_name)
    PowerModels.standardize_cost_terms!(data, order=2)

    num_buses = length(data["bus"])

    # Build mapping
    bus_index, bus2idx = _sort_keys(data["bus"], "bus_i")

    buses = Array{Bus}(undef, num_buses)
    bus_ref=-1
    for idx in bus_index
        i = bus2idx[idx]
        bus = data["bus"][string(idx)]
        @assert bus["bus_i"] > 0  #don't support nonpositive bus ids
        @assert bus["bus_i"] == idx
        buses[i] = Bus(
            bus["bus_i"],
            bus["bus_type"],
            0.0, 0.0, 0.0, 0.0,
            bus["area"],
            bus["vm"],
            bus["va"],
            bus["base_kv"],
            bus["zone"],
            bus["vmax"],
            bus["vmin"],
        )
        buses[i].Va *= pi/180
        if buses[i].bustype==3
            if bus_ref>0
                error("More than one reference bus present in the data")
            else
                bus_ref=i
            end
        end
    end

    for (_, load) in data["load"]
        load_bus = load["load_bus"]
        idx = bus2idx[load_bus]
        buses[idx].Pd = load["pd"] * baseMVA
        buses[idx].Qd = load["qd"] * baseMVA
    end

    for (_, shunt) in data["shunt"]
        shunt_bus = shunt["shunt_bus"]
        idx = bus2idx[shunt_bus]
        buses[idx].Gs = shunt["gs"]
        buses[idx].Bs = shunt["bs"]
    end

    # Get active branches
    active_branch = Int[]
    for (_, branch) in data["branch"]
        if branch["br_status"] == 1
            push!(active_branch, branch["index"])
        end
    end
    num_lines = length(active_branch)
    lines = Array{Line}(undef, num_lines)

    for (lit, idx) in enumerate(active_branch)
        branch = data["branch"][string(idx)]
        @assert branch["br_status"] == 1  # should be on since we discarded all other
        lines[lit] = Line(
            branch["f_bus"],
            branch["t_bus"],
            branch["br_r"],
            branch["br_x"],
            branch["b_fr"] + branch["b_to"],
            get(branch, "rate_a", 0.0) * baseMVA,
            get(branch, "rate_b", 0.0) * baseMVA,
            get(branch, "rate_c", 0.0) * baseMVA,
            branch["tap"],
            branch["shift"],
            branch["br_status"],
            branch["angmin"],
            branch["angmax"],
        )
    end

    # Get active generators
    active_gens = Int[]
    for (_, gen) in data["gen"]
        if gen["gen_status"] == 1
            push!(active_gens, gen["index"])
        end
    end
    num_gens = length(active_gens)
    generators = Array{Gener}(undef, num_gens)
    for (i, idx) in enumerate(active_gens)
        gen = data["gen"][string(idx)]
        @assert gen["gen_status"] == 1
        generators[i] = Gener(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, Array{Int}(undef, 0))

        generators[i].bus      = gen["gen_bus"]
        generators[i].Pg       = gen["pg"]
        generators[i].Qg       = gen["qg"]
        generators[i].Qmax     = gen["qmax"]
        generators[i].Qmin     = gen["qmin"]
        generators[i].Vg       = gen["vg"]
        generators[i].mBase    = gen["mbase"]
        generators[i].status   = gen["gen_status"]
        generators[i].Pmax     = gen["pmax"]
        generators[i].Pmin     = gen["pmin"]
        generators[i].Pc1      = get(gen, "pc1", 0.0)
        generators[i].Pc2      = get(gen, "pc2", 0.0)
        generators[i].Qc1min   = get(gen, "qc1min", 0.0)
        generators[i].Qc1max   = get(gen, "qc1max", 0.0)
        generators[i].Qc2min   = get(gen, "qc2min", 0.0)
        generators[i].Qc2max   = get(gen, "qc2max", 0.0)
        # generators[i].gentype  = costgen_arr[git,1]
        generators[i].startup  = gen["startup"]
        generators[i].shutdown = gen["shutdown"]
        generators[i].n        = gen["ncost"]
        generators[i].coeff    = gen["cost"] ./ [baseMVA^2, baseMVA, 1.0]
    end

    # build a dictionary between buses ids and their indexes
    busIdx = mapBusIdToIdx(buses)

    # set up the FromLines and ToLines for each bus
    FromLines,ToLines = mapLinesToBuses(buses, lines, busIdx)

    # generators at each bus
    BusGeners = mapGenersToBuses(buses, generators, busIdx)

    #println(generators)
    #println(bus_ref)
    return OPFData(buses, lines, generators, bus_ref, baseMVA, busIdx, FromLines, ToLines, BusGeners)
end

function _opf_loaddata_legacy(case_name, lineOff=Line())
    #
    # load buses
    #
    # bus_arr = readdlm("data/" * case_name * ".bus")
    bus_arr = readdlm(case_name * ".bus")
    num_buses = size(bus_arr,1)
    buses = Array{Bus}(undef, num_buses)
    bus_ref=-1
    for i in 1:num_buses
        @assert bus_arr[i,1]>0  #don't support nonpositive bus ids
        buses[i] = Bus(bus_arr[i,1:13]...)
        buses[i].Va *= pi/180
        if buses[i].bustype==3
            if bus_ref>0
                error("More than one reference bus present in the data")
            else
                bus_ref=i
            end
        end
        #println("bus ", i, " ", buses[i].Vmin, "      ", buses[i].Vmax)
    end

    #
    # load branches/lines
    #
    # branch_arr = readdlm("data/" * case_name * ".branch")
    branch_arr = readdlm(case_name * ".branch")
    num_lines = size(branch_arr,1)
    lines_on = findall((branch_arr[:,11].>0) .& ((branch_arr[:,1].!=lineOff.from) .| (branch_arr[:,2].!=lineOff.to)) )
    num_on   = length(lines_on)

    if lineOff.from>0 && lineOff.to>0
        println("opf_loaddata: was asked to remove line from,to=", lineOff.from, ",", lineOff.to)
        #println(lines_on, branch_arr[:,1].!=lineOff.from, branch_arr[:,2].!=lineOff.to)
    end
    if length(findall(branch_arr[:,11].==0))>0
        println("opf_loaddata: ", num_lines-length(findall(branch_arr[:,11].>0)), " lines are off and will be discarded (out of ", num_lines, ")")
    end

    lines = Array{Line}(undef, num_on)

    lit=0
    for i in lines_on
        @assert branch_arr[i,11] == 1  #should be on since we discarded all other
        lit += 1
        lines[lit] = Line(branch_arr[i, 1:13]...)
        #=
        if (lines[lit].angmin != 0 || lines[lit].angmax != 0) && (lines[lit].angmin>-360 || lines[lit].angmax<360)
        println("Voltage bounds on line ", i, " with angmin ", lines[lit].angmin, " and angmax ", lines[lit].angmax)
        error("Bounds of voltage angles are still to be implemented.")
        end
        =#

    end
    @assert lit == num_on

    #
    # load generators
    #
    # gen_arr = readdlm("data/" * case_name * ".gen")
    gen_arr = readdlm(case_name * ".gen")
    # costgen_arr = readdlm("data/" * case_name * ".gencost")
    costgen_arr = readdlm(case_name * ".gencost")
    num_gens = size(gen_arr,1)

    baseMVA=100

    @assert num_gens == size(costgen_arr,1)

    gens_on=findall(x->x!=0, gen_arr[:,8]); num_on=length(gens_on)
    if num_gens-num_on>0
        println("loaddata: ", num_gens-num_on, " generators are off and will be discarded (out of ", num_gens, ")")
    end

    generators = Array{Gener}(undef, num_on)
    i=0
    for git in gens_on
        i += 1
        generators[i] = Gener(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, Array{Int}(undef, 0)) #gen_arr[i,1:end]...)

        generators[i].bus      = gen_arr[git,1]
        generators[i].Pg       = gen_arr[git,2] / baseMVA
        generators[i].Qg       = gen_arr[git,3] / baseMVA
        generators[i].Qmax     = gen_arr[git,4] / baseMVA
        generators[i].Qmin     = gen_arr[git,5] / baseMVA
        generators[i].Vg       = gen_arr[git,6]
        generators[i].mBase    = gen_arr[git,7]
        generators[i].status   = gen_arr[git,8]
        @assert generators[i].status==1
        generators[i].Pmax     = gen_arr[git,9]  / baseMVA
        generators[i].Pmin     = gen_arr[git,10] / baseMVA
        generators[i].Pc1      = gen_arr[git,11]
        generators[i].Pc2      = gen_arr[git,12]
        generators[i].Qc1min   = gen_arr[git,13]
        generators[i].Qc1max   = gen_arr[git,14]
        generators[i].Qc2min   = gen_arr[git,15]
        generators[i].Qc2max   = gen_arr[git,16]
        generators[i].gentype  = costgen_arr[git,1]
        generators[i].startup  = costgen_arr[git,2]
        generators[i].shutdown = costgen_arr[git,3]
        generators[i].n        = costgen_arr[git,4]
        @assert(generators[i].n <= 3 && generators[i].n >= 2)
        if generators[i].gentype == 1
            generators[i].coeff = costgen_arr[git,5:end]
            error("Piecewise linear costs remains to be implemented.")
        else
            if generators[i].gentype == 2
                generators[i].coeff = costgen_arr[git,5:end]
                #println(generators[i].coeff, " ", length(generators[i].coeff), " ", generators[i].coeff[2])
            else
                error("Invalid generator cost model in the data.")
            end
        end
    end

    # build a dictionary between buses ids and their indexes
    busIdx = mapBusIdToIdx(buses)

    # set up the FromLines and ToLines for each bus
    FromLines,ToLines = mapLinesToBuses(buses, lines, busIdx)

    # generators at each bus
    BusGeners = mapGenersToBuses(buses, generators, busIdx)

    #println(generators)
    #println(bus_ref)
    return OPFData(buses, lines, generators, bus_ref, baseMVA, busIdx, FromLines, ToLines, BusGeners)
end

function opf_loaddata(case_name, lineOff=Line())
    if endswith(case_name, ".m") || endswith(case_name, ".raw")
        _opf_loaddata_mat(case_name, lineOff)
    else
        _opf_loaddata_legacy(case_name, lineOff)
    end
end

function computeAdmitances(lines, buses, baseMVA, device; VI=Array{Int}, VD=Array{Float64})
    nlines = length(lines)
    YffR=Array{Float64}(undef, nlines)
    YffI=Array{Float64}(undef, nlines)
    YttR=Array{Float64}(undef, nlines)
    YttI=Array{Float64}(undef, nlines)
    YftR=Array{Float64}(undef, nlines)
    YftI=Array{Float64}(undef, nlines)
    YtfR=Array{Float64}(undef, nlines)
    YtfI=Array{Float64}(undef, nlines)

    for i in 1:nlines
        @assert lines[i].status == 1
        Ys = 1/(lines[i].r + lines[i].x*im)
        #assign nonzero tap ratio
        tap = lines[i].ratio==0 ? 1.0 : lines[i].ratio

        #add phase shifters
        tap *= exp(lines[i].angle * pi/180 * im)

        Ytt = Ys + lines[i].b/2*im
        Yff = Ytt / (tap*conj(tap))
        Yft = -Ys / conj(tap)
        Ytf = -Ys / tap

        #split into real and imag parts
        YffR[i] = real(Yff); YffI[i] = imag(Yff)
        YttR[i] = real(Ytt); YttI[i] = imag(Ytt)
        YtfR[i] = real(Ytf); YtfI[i] = imag(Ytf)
        YftR[i] = real(Yft); YftI[i] = imag(Yft)
        #@printf("[%4d]  tap=%12.9f   %12.9f\n", i, real(tap), imag(tap));
    end

    nbuses = length(buses)
    YshR = Array{Float64}(undef, nbuses)
    YshI = Array{Float64}(undef, nbuses)
    for i in 1:nbuses
        YshR[i] = buses[i].Gs / baseMVA
        YshI[i] = buses[i].Bs / baseMVA
        #@printf("[%4d]   Ysh  %15.12f + %15.12f i \n", i, YshR[i], YshI[i])
    end

    @assert 0==length(findall(isnan.(YffR)))+length(findall(isinf.(YffR)))
    @assert 0==length(findall(isnan.(YffI)))+length(findall(isinf.(YffI)))
    @assert 0==length(findall(isnan.(YttR)))+length(findall(isinf.(YttR)))
    @assert 0==length(findall(isnan.(YttI)))+length(findall(isinf.(YttI)))
    @assert 0==length(findall(isnan.(YftR)))+length(findall(isinf.(YftR)))
    @assert 0==length(findall(isnan.(YftI)))+length(findall(isinf.(YftI)))
    @assert 0==length(findall(isnan.(YtfR)))+length(findall(isinf.(YtfR)))
    @assert 0==length(findall(isnan.(YtfI)))+length(findall(isinf.(YtfI)))
    @assert 0==length(findall(isnan.(YshR)))+length(findall(isinf.(YshR)))
    @assert 0==length(findall(isnan.(YshI)))+length(findall(isinf.(YshI)))

    if isa(device, KA.GPU)
        return copyto!(VD(undef, nlines), 1, YffR, 1, nlines),
        copyto!(VD(undef, nlines), 1, YffI, 1, nlines),
        copyto!(VD(undef, nlines), 1, YttR, 1, nlines),
        copyto!(VD(undef, nlines), 1, YttI, 1, nlines),
        copyto!(VD(undef, nlines), 1, YftR, 1, nlines),
        copyto!(VD(undef, nlines), 1, YftI, 1, nlines),
        copyto!(VD(undef, nlines), 1, YtfR, 1, nlines),
        copyto!(VD(undef, nlines), 1, YtfI, 1, nlines),
        YshR, YshI
    else
        return YffR, YffI, YttR, YttI, YftR, YftI, YtfR, YtfI, YshR, YshI
    end
end

# Builds a map from lines to buses.
# For each line we store an array with zero or one element containing
# the  'From' and 'To'  bus number.
function mapLinesToBuses(buses, lines, busDict)
    nbus = length(buses)
    FromLines = [Int[] for b in 1:nbus]
    ToLines   = [Int[] for b in 1:nbus]
    for i in 1:length(lines)
        busID = busDict[lines[i].from]
        @assert 1<= busID <= nbus
        push!(FromLines[busID], i)

        busID = busDict[lines[i].to]
        @assert 1<= busID  <= nbus
        push!(ToLines[busID], i)
    end

    return FromLines,ToLines
end

# Builds a mapping between bus ids and bus indexes
#
# Returns a dictionary with bus ids as keys and bus indexes as values
function mapBusIdToIdx(buses)
    dict = Dict{Int,Int}()
    for b in 1:length(buses)
        @assert !haskey(dict,buses[b].bus_i)
        dict[buses[b].bus_i] = b
    end
    return dict
end

# Builds a map between buses and generators.
# For each bus we keep an array of corresponding generators number (as array).
#
# (Can be more than one generator per bus)
function mapGenersToBuses(buses, generators,busDict)
    gen2bus = [Int[] for b in 1:length(buses)]
    for g in 1:length(generators)
        busID = busDict[ generators[g].bus ]
        #@assert(0==length(gen2bus[busID])) #at most one generator per bus
        push!(gen2bus[busID], g)
    end
    return gen2bus
end


text:
dotestmean   = true
doprintmeans = false
reltol       = .01          # required test accuracy
algs      = (Direct(),)

# run the given number of SSAs and return the mean
function runSSAs(jump_prob, Nsims, idx)
    Psamp = zeros(Int, Nsims)
    for i in 1:Nsims
        sol = solve(jump_prob, SSAStepper())
        Psamp[i] = sol[idx,end]
    end
    mean(Psamp)
end

function execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, idx, test_name)
    for method in algs
        jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
        avg_val = runSSAs(jump_prob, Nsims, idx)

        if dotestmean
            if doprintmeans
                println(test_name, ", method = ", typeof(method), ", mean = ", avg_val, ", act_mean = ", expected_avg)
            end
            @test abs(avg_val - expected_avg) < reltol * expected_avg
        end
    end

end

# full macro
Nsims = 32000
tf = .01
u0 = [200, 100, 150]
expected_avg = 84.876015624999994
rs = @reaction_network dtype begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
rates = [1., 2., .5, .75, .25]
prob = DiscreteProblem(u0, (0.0, tf), rates)
execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test (no syms)")

prob = DiscreteProblem(rs, u0, (0.0, tf), rates)
execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test (syms)")

# min_network macro
rs = @min_reaction_network dtype begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
addjumps!(rs)
prob = DiscreteProblem(u0, (0.0, tf), rates)
execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test, min network (no syms)")

prob = DiscreteProblem(rs, u0, (0.0, tf), rates)
execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test, min network (syms)")


text:
# Sinus test

function prepare_input(input_size)
	return linspace(0, 2pi, input_size)
end

function test_function(input_data)
	return sin.(input_data)
end


text:
module TestDocumentation

using Base.Test
using DataFrames
using TimeData

println("\n Running documentation tests\n")

using WorldBankDataTd

## single indicator, single country
gnp = WorldBankDataTd.wdi("NY.GNP.PCAP.CD", "BR")

## single indicator, multiple countries
gnp = WorldBankDataTd.wdi("NY.GNP.PCAP.CD", ["BR", "US", "DE"])

## multiple indicators, single country
data = WorldBankDataTd.wdi(["NY.GNP.PCAP.CD", "SP.DYN.LE00.IN"],
                         "BR")

## multiple indicators, multiple countries
data = WorldBankDataTd.wdi(["NY.GNP.PCAP.CD", "SP.DYN.LE00.IN"],
                         ["BR", "US", "DE"])

## same, but as DataFrame
data = WorldBankDataTd.wdi(["NY.GNP.PCAP.CD", "SP.DYN.LE00.IN"],
                         ["BR", "US", "DE"],
                         format = DataFrame)

## multiple indicators, multiple countries, additional information
data = WorldBankDataTd.wdi(["NY.GNP.PCAP.CD", "SP.DYN.LE00.IN"],
                         ["BR", "US", "DE"], extra = true)

data[1:5, :]

countryData = getWBMeta("countries")
countryData[1:5, :]

names(countryData)

size(countryData)

indicatorData = getWBMeta("indicators")
indicatorData[1:5, :]

names(indicatorData)

size(indicatorData)

WorldBankDataTd.country_cache[1:5, :]

WorldBankDataTd.indicator_cache[1:5, :]

res = search_wdi("countries", :name, r"united"i)
res

res = search_wdi("indicators", :description, r"gross national expenditure"i)
res[:name]

search_wdi("countries", :iso2c, r"TZ"i)
search_wdi("countries", :income, r"upper middle"i)
search_wdi("countries", :region, r"Latin America"i)
search_wdi("countries", :capital, r"^Ka"i)
search_wdi("countries", :lending, r"IBRD"i)
search_wdi("indicators", :name, r"gross national expenditure"i)
search_wdi("indicators", :description, r"gross national expenditure"i)
search_wdi("indicators", :source_database, r"Sustainable"i)
search_wdi("indicators", :source_organization,
           r"Global Partnership"i)[1:5, :]

data = wdi("NY.GNP.PCAP.CD", ["US","BR"], 1980, 2012, extra = true)
usData = chkDates(x-> x[:iso2c] .== "US", eachdate(data)) |>
         x -> asArr(x, Bool, false) |>
         x -> data[x[:], :]
usData

dfAS = wdi("EN.ATM.CO2E.KT", "AS")

end


text:
"""
Implementation of a Rauch-Tung-Striebel smoother based on
Ch.8 of Bayesian Filtering & Smoothing (Särkka, 2014).

Author: Wouter Kouw
Last update: 03-07-2021
"""

using Distributions
using Random

function fixedinterval_rts_smoother(observations,
                                    transition_matrix,
                                    emission_matrix,
                                    process_noise,
                                    measurement_noise,
                                    state0)
    """
    Fixed-interval Rauch-Tung-Striebel smoother (Th. 8.2)

    This filter is built for a linear Gaussian dynamical system with known
    transition coefficients, process and measurement noise.
    """

    # Dimensionality
    Dx = size(process_noise,1)
    Dy = size(measurement_noise,1)

    # Recast process noise to matrix
    if Dx == 1
        if typeof(process_noise) != Array{Float64,2}
            process_noise = reshape([process_noise], 1, 1)
        end
        if typeof(measurement_noise) != Array{Float64,2}
            measurement_noise = reshape([measurement_noise], 1, 1)
        end
    end

    # Time horizon
    time_horizon = length(observations)

    # Initialize estimate arrays
    mk = zeros(Dx, time_horizon)
    Pk = zeros(Dx, Dx, time_horizon)

    # Initial state prior
    m_0, P_0 = state0

    # Start previous state variable
    m_kmin = m_0
    P_kmin = P_0

    "Forward pass"
    for k = 1:time_horizon

        # Forward prediction step
        m_k_pred = transition_matrix*m_kmin
        P_k_pred = transition_matrix*P_kmin*transition_matrix' .+ process_noise

        # Forward update step
        v_k = observations[:,k] .- emission_matrix*m_k_pred
        S_k = emission_matrix*P_k_pred*emission_matrix' .+ measurement_noise
        K_k = P_k_pred*emission_matrix'*inv(S_k)
        m_k = m_k_pred .+ K_k*v_k
        P_k = P_k_pred .- K_k*S_k*K_k'

        # Store estimates
        mk[:,k] = m_k
        Pk[:,:,k] = P_k

        # Update previous state variable
        m_kmin = m_k
        P_kmin = P_k

    end

    # Initialize smoothing estimate arrays
    msk = zeros(Dx, time_horizon)
    Psk = zeros(Dx, Dx, time_horizon)

    # Smoothed estimates at time horizon
    msk[:,time_horizon] = mk[:,time_horizon]
    Psk[:,:,time_horizon] = Pk[:,:,time_horizon]

    "Backward pass"
    for k = time_horizon-1:-1:1

        # Backward prediction
        m_kplus = transition_matrix * mk[:,k]
        P_kplus = transition_matrix * Pk[:,:,k] * transition_matrix' .+ process_noise

        # Backward update step
        G_k = Pk[:,:,k]*transition_matrix' * inv(P_kplus)
        msk[:,k] = mk[:,k] + G_k*(msk[:,k+1] - m_kplus)
        Psk[:,:,k] = Pk[:,:,k] + G_k*(Psk[:,:,k+1] - P_kplus)*G_k'

    end

    return msk, Psk
end

max_stars_repo_name: wmkouw/Bayesian-smoothers

max_issues_repo_name: wmkouw/Bayesian-smoothers

max_forks_repo_name: wmkouw/Bayesian-smoothers


text:
# Solution to Q2a (implement a Linear module)
import Base

type Linear{T}
    # Parameters of the module
    w::Matrix{T}
    b::Vector{T}

    # (Accumulated) gradients of the parameters and input after backprop
    ∂w::Matrix{T}
    ∂b::Vector{T}
    ∂x::Matrix{T}

    # Buffer of the input to fprop, needed during bprop.
    x::Matrix{T}

    Linear(n_in, n_out ; σ=0) = n_in < 1 || n_out < 1 ? error("Invalid dimensions of Linear") : new(
        map(T, σ*randn(n_in, n_out)),  # w
        zeros(T, n_out),               # b

        zeros(T, n_in, n_out),  # ∂w
        zeros(T, n_out),        # ∂b

        zeros(T, 0, 0),  # ∂x
        zeros(T, 0, 0),  # x
    )
end

function fprop!{T}(l::Linear{T}, x::Matrix{T})
    # Remember the input for later use during bprop.
    l.x = x

    # Compute and return the output.
    return x * l.w .+ l.b'
end

function bprop!{T}(l::Linear{T}, ∂y::Matrix{T})
    # Accumulate the gradients wrt. W and b
    l.∂w[:,:] += l.x' * ∂y
    l.∂b[:] += sum(∂y, 1)'

    # Store and return the gradient wrt. the inputs
    return l.∂x = ∂y * l.w'
end

params(l::Linear) = (l.w, l.b)
grads(l::Linear) = (l.∂w, l.∂b)

Base.show(io::IO, l::Linear) = print(io, "Linear", size(l.w))