text:
module pAdicConj

#export conjugates_pAdic_log

#= XXX: qadics in flint:
   defining polynomial needs to be conway (currently) as this is the only constructor
   code assumes always the def. poly is sparse and has small (<=p) coefficients
   (at least for performance)
   I don't know if this implies its pointless to use them?
   At least the general use in Hecke is not quite clear.
=#

using Hecke

function (O::Generic.PolyRing{padic})(f::fmpz_poly)
  g = O()
  R = base_ring(O)
  for i=0:length(f)
    Hecke.setcoeff!(g, i, R(coeff(f, i)))
  end
  return g
end

mutable struct pAdic_root_ctx  #TODO: parametrize to allow non-simples as well
  H::Hecke.HenselCtx
  prec::Int
  Q#::Array{Generic.ResRing{fmpz_poly}, 1} #TODO: FlintQadicField (to be wrapped)
  pN::fmpz
  Qp::FlintPadicField
  Qpt
  cache::Dict{nf_elem, Array{Any, 1}}
  function pAdic_root_ctx(K::AnticNumberField, p::Int)
    d = Dict{Int, pAdic_root_ctx}()
    try
      d = _get_nf_conjugate_data_pAdic(K)
      @show "usig old ctx"
    catch e
      #check proper type of e
      @show e
      _set_nf_conjugate_data_pAdic(K, d)
    end

    if haskey(d, p)
      return d[p]
    end

    r = new()
    @show "new r"
    zx, x = PolynomialRing(FlintZZ, cached = false)
    f = zx(K.pol * denominator(K.pol))
    r.H = Hecke.HenselCtx(f, fmpz(p))
    Hecke.start_lift(r.H, 10)
    r.prec = 10
    r.pN = fmpz(p)^10
    r.Qp = Hecke.FlintPadicField(p, 10)
    r.Qpt = PolynomialRing(r.Qp, cached = false)[1]

    d[p] = r
    s = Hecke.factor_to_dict(r.H.LF)
    r.Q = [ResidueRing(r.Qpt, r.Qpt(f)) for f = keys(s)]
    return r
  end
end

t = Hecke.create_accessors(AnticNumberField, Dict{Int, pAdic_root_ctx}, get_handle())
global _get_nf_conjugate_data_pAdic = t[1]
global _set_nf_conjugate_data_pAdic = t[2]

function assure_precision(r::pAdic_root_ctx, n::Int)
  if r.prec < n
    Hecke.continue_lift(r.H, n)
    r.prec = n
    s = Hecke.factor_to_dict(r.H.LF)
    #think: order may change as this is a dictionary
    p = r.Qp.p #see warning above
    r.pN = fmpz(p)^n
    r.Qp = Hecke.FlintPadicField(p, n)
    r.Qpt = PolynomialRing(r.Qp, cached = false)[1]
    s = Hecke.factor_to_dict(r.H.LF)
    r.Q = [ResidueRing(r.Qpt, r.Qpt(f)) for f = keys(s)]
  end
end

function Hecke.prime(R::PadicField, i::Int)
  p = fmpz()
  ccall((:padic_ctx_pow_ui, :libflint), Cvoid, (Ref{fmpz}, Int, Ref{PadicField}), p, i, R)
  return p
end

function frobenius_gen(R, e::Int = 1)
  Qpt = base_ring(R)
  Qp = base_ring(Qpt)
  f = R.modulus
  #Frob is a lift of x -> x^(p^e) this is a root of f...
  fs = derivative(f)
  i = gen(R)^prime(Qp, e)
  #@show methods(inv, (typeof(i), ))
  o = inv(fs(i))
  pr = 1
  #TODO: use a lifting chain and adjust the precision
  while pr < Qp.prec_max
    i = i-f(i)*o
    o *= (2-o*fs(i))
    parent(o) #julia bug: without this, it dies
    pr *= 2
  end
  return i
end

function frobenius(a, g = frobenius_gen(parent(a)))
  return a.data(g) # should be compose_mod
end

function conjugates_pAdic(a::nf_elem, p::Int, n::Int = 10; orbit::Bool = true)
  r = pAdic_root_ctx(parent(a), p)

  c = []
  assure_precision(r, n)
  Zx = PolynomialRing(FlintZZ, cached = false)[1]
  for C = r.Q
    b = Zx(a*denominator(a))
    x = C(r.Qpt(b)*inv(r.Qp(denominator(a))))
#    di = invmod(denominator(a), r.pN)
#    x *= di
    push!(c, x)
    if orbit
      g = frobenius_gen(parent(x))
      for i=2:degree(parent(x).modulus)
        push!(c, frobenius(c[end], g))
      end
    end
  end
  return c
end

#= log(x)' = 1/x
   1/(1-x) = sum_0 x^i
   log(1-x) = sum x^i/i
   =#
function _valuation(a)
  return minimum(valuation(coeff(a.data, i)) for i=0:degree(a.data))
end

#TODO/ ToThink
#  - if v_p(1-a) = l, then v_p((1-a)^p) = l+1, so
#    by raising to a power of p, the series converges faster
#    aparently, optimal is p^sqrt(prec)
#    if p >> 0, I don't think this is fast
#  - pari: 1/2 log((1+x)/(1-x)) = atanh(x), thus
#    1/2 log(x) = atanh((x-1)/(x+1))
#    atanh(x) = sum x^(2l+a)/(2l+a)
#    converges faster.
#  - evaluation via Brent-Kung rather than naive
#  - the sum is not long enough for small primes as the denominator
#    will have valuations
#  - needs to be optimized...
#  - maybe use Teichmueller lift instead of raising to the ?-1 st power:
#    a = tb for t the Teichmueller and b=1(pi)
#    in particular, t is torsion, hence log(t) = 1
#    There are non-trivial, fast(?) methods to compute t
#    The naive is t = a^((p^n)^prec) which is even worse...
function _log(a)
  Qpt = base_ring(parent(a))
  Qp = base_ring(Qpt)
  #log a = log(1-(1-a))
  @assert _valuation(a) == 0
  ex = prime(Qp, degree(modulus(parent(a)))) -1
  a = a^ex
  # .. this guarantees a equiv 1 mod p, hence the power series converges
  a = 1-a
  @assert all(x->valuation(coeff(a.data, x)) > 0, 0:degree(a.data))
  n = Qp.prec_max
  pa = a
  l = a
  i = 2
  z = ceil(Int, log(Int(prime(base_ring(a.data))), n))
  while any(j-> coeff(pa.data, j).v <= n, 0:degree(pa.data))
    pa *= a
    for j=0:degree(pa.data)
      c = coeff(pa.data, j)
      c.N = min(n+z, c.N)
    end
    l += pa*inv(Qp(i))
    i += 1
  end
  return l*inv(Qp(ex))
end

function _log2(a)
  Qpt = base_ring(parent(a))
  Qp = base_ring(Qpt)
  #log a = log(1-(1-a))
  @assert _valuation(a) == 0
  ex = prime(Qp, degree(modulus(parent(a)))) -1
  a = a^ex
  # .. this guarantees a equiv 1 mod p, hence the power series converges
  a = 1-a
  @assert all(x->valuation(coeff(a.data, x)) > 0, 0:degree(a.data))
  n = Qp.prec_max
  pa = a
  l = a
  for i=2:n
    pa *= a
    l += pa*inv(Qp(i))
  end
  return l*inv(Qp(ex))
end

#= XXX: this takes only one value for each completion and returns the coeff vec
#     is this enough or do I need to work in a splitting field??
  Yep: we should use log(u^(1)) ... log(u^(n)) for all conjugates, but
    log commutes with frob and thus any relation in the (say) 1st conjugate automatically
    implies one (via frob) for the second
    We are looking for the Qp-rank, NOT the q-adic one!
    This is so neat!
=#
function _conjugates_pAdic_log(a::nf_elem, p::Int, n::Int = 10; normalise::Bool = false)
  r = 1
  try
    s = _get_nf_conjugate_data_pAdic(parent(a))
    if haskey(s, p)
      r = s[p]
    end
    if !isdefined(r, :cache)
      @show "create cache 1"
      r.cache = Dict{nf_elem, Array{Any, 1}}()
    end

    if haskey(r.cache, a)
      la = r.cache[a]
      if precision(coeff(la[1].data, 0)) >= n
        @show "retrieve"
        return la
      end
      @show "bad"
    else
      @show "new elt"
    end
  catch t
    @show "errror 1", t
  end
  c = conjugates_pAdic(a, p, n, orbit = false)
  if normalise && any(x -> _valuation(x) != 0, c)
    lp = prime_decomposition(maximal_order(parent(a)), p)
    b = prod([p[1].gen_two^valuation(a, p[1]) for p = lp])
    a *= inv(elem_in_nf(b))
    c = conjugates_pAdic(a, p, n, orbit = false)
  end
  l = []
  for x = c
    y = _log(x)
    push!(l, y)
  end
  r = _get_nf_conjugate_data_pAdic(parent(a))[p]
  if !isdefined(r, :cache)
    @show "create cache 2"
    r.cache = Dict{nf_elem, Array{Any, 1}}()
  end
  r.cache[a] = l
  return l
end

function conjugates_pAdic_log(a::nf_elem, p::Int, n::Int = 10; normalise::Bool = false)
  d = _conjugates_pAdic_log(a, p, n, normalise = normalise)
  l = []
  for y = d
    for i=0:degree(modulus(parent(y)))-1
      push!(l, coeff(y.data, i))
    end
  end
  return matrix(parent(l[1]), 1, length(l), l)
end
#TODO: caching padic logs in the number field?
#TODO: adjust the precision downwards

function conjugates_pAdic_log(a::FacElem{nf_elem, AnticNumberField}, p::Int, n::Int = 10)
  l = 1
  #if facelem is a (s) unit coprime to p, then one can/ could/should deal with
  #valuations differently: mult by suitable valuation elements - their
  #contribution must cancel
  #done
  for (k, v) = a.fac
    if l === 1
      l = v*conjugates_pAdic_log(k, p, n, normalise=true)
    else
      l += v*conjugates_pAdic_log(k, p, n, normalise=true)
    end
  end
  return l
end

#not working, no exp...
function conjugates_pAdic(a::FacElem{nf_elem, AnticNumberField}, p::Int, n::Int = 10)

  #if facelem is a (s) unit coprime to p, then one can/ could/should deal with
  #valuations differently: mult by suitable valuation elements - their
  #contribution must cancel
  #done
  l = _conjugates_pAdic_log(base_ring(a)(1), p, n, normalise = true)
  for (k, v) = a.fac
    iszero(v) && continue
    if l === 1
      l = v*_conjugates_pAdic_log(k, p, n, normalise=true)
    else
      l += v*_conjugates_pAdic_log(k, p, n, normalise=true)
    end
  end
  return [_exp(x) for x = l]
end

function conjugates_pAdic_log(a::NfOrdElem, p::Int, n::Int = 10)
  return conjugates_pAdic_log(a.elem_in_nf, p, 10)
end

#TODO: new name!
function getUnit(a::padic)
  u = fmpz()
  ccall((:fmpz_set, :libflint), Cvoid, (Ref{fmpz}, Ref{Int}), u, a.u)
  return u, a.v, a.N
end

function lift_reco(::FlintRationalField, a::padic; reco::Bool = false)
  if reco
    u, v, N = getUnit(a)
    R = parent(a)
    fl, c, d = rational_reconstruction(u, prime(R, N-v))
    !fl && return false
    if false && 2*max(nbits(c), nbits(d)) > nbits(prime(R, N-v)) -20 #arbitrary
      @show "bad"
      return false
    end
    @assert fl
    x = FlintQQ(c, d)
    if v < 0
      return x//prime(R, -v)
    else
      return x*prime(R, v)
    end
  else
    return lift(FlintQQ, a)
  end
end

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

import Base.*

function *(A::fmpz_mat, B::MatElem{padic})
  return matrix(base_ring(B), A) * B
end

#= kept for the comments

function mult_syzygies_units(a::Array{FacElem{nf_elem, AnticNumberField}, 1})
  p = next_prime(2^10) #experimentally, the runtime is dominated by
         # log which in case is dominated by the a^(p^n-1) in the 1st step
         # thus try p smaller..., ideally also try n smaller...
         # also, see comments at _log
  u = FacElem{nf_elem, AnticNumberField}[]
  la = [conjugates_pAdic_log(e, p, 300) for e = a] #can loose precision
    # needs to be traced
    # precision needs to be updated.
  n = ncols(la[1])
  Qp = base_ring(la[1])
  lu = matrix(Qp, 0, n, [])
  for i=1:length(la)
    if iszero(la[i])
      println("torsino at $i\n")
      continue
    end
    lv = vcat(lu, la[i])
    k = Hecke.left_kernel_basis(lv)
    @assert length(k) < 2
    if length(k) == 0
      println("new at $i")
      push!(u, a[i])
      lu = vcat(lu, la[i])
    else # length == 1 extend the module
      r = [lift_reco(FlintQQ, x, reco = true) for x = k[1]]
      #lift can fail if precision wrong.
      #or at least result is (can be) garbage
      #= bound on relation
        u has s independent units (independence is failproof p-adically)
        v is new and (probably, up to verification of the relation) dependent
        Then <u, v> = <e_1, .., e_s> where the e_'s are unknown.
        we have u = U(e_1, ..., e_n) for U in Z^(s x s)
                v = V(e_1, ..., e_n) vor V in Z^(1 x s)
        Over Q: v = x(u_1, ..., u_s) thus
            V(e_1, ..., e_s) = xU(e_1, ..., e_s)
        and, since the e_i are independent
            V = xU
        Cramer's rule
            x_i = det(U_1, ..., U_i-1, V, ... U_s) / det(U)
        and we need bounds for the determinants.
        As an example, det U
        det U = <e_1, ..., e_s> : <u_1, ..., u_s>
        and using universal lower bounds on the size of units (Dobrowski)
        and the successive minimal, we can get a lower bound on the
        regulator of <e_i>. Hadramat gives an upper bound on reg(<u_i>)
        (regulator in the sence of lattice disciminant)

        Think: can we use the p-adic regulator to help???
               possibly increase precision until we either have
               indepence or a relation
               ignore bounds?
      =#
      d = reduce(lcm, map(denominator, r))
      gamma = [FlintZZ(x*d) for x = r]
      #technically, this relations needs to be verified.
      #=
        we have a relation mod p^k, which means
          e = prod gamma[i]*u[i]
        has all logarithms == 0 mod p^k
        which should mean Norm(e) = product of the local norms
                                  = prod of exp(trace(components))
        For p large enough (larger than precision),
        val(exp(x)) = val(1-x) I think
        so
                                  >= p^(k* # primes) = p^l
        Now, either e is torsion (ie. logs identically 0) or
        not, but then arithmetic means:
        N(e) <= (T_2(e)^(1/2)/n)^n
        implies a non-trivial lower bound on the T_2:
        T_2(e) >= (np^(l/n))^2 =: B
        If T_2(e) < B, then e is torsion. B is MUCH larger than
        Dobrowski. Hence this can be evaluated with low precision.

        Not done.
      =#
      @assert reduce(gcd, gamma) == 1 # should be a primitive relation
      _, U = hnf_with_trafo(matrix(FlintZZ, length(r), 1, gamma))
      U = inv(U)
      U = sub(U, 1:nrows(U), 2:ncols(U))
      #new basis is the cols of U
      push!(u, a[i])
      u = Hecke._transform(u, U)
      lu = U'*lv
    end
  end
  return u
end

=#

function verify_gamma(a::Array{FacElem{nf_elem, AnticNumberField}, 1}, g::Array{fmpz, 1}, v::fmpz)
  #knowing that sum g[i] log(a[i]) == 0 mod v, prove that prod a[i]^g[i] is
  #torsion
  #= I claim N(1-a) > v^n for n the field degree:
   Let K be one of the p-adic fields involved, set b = a^g
   then log(K(b)) = 0 (v = p^l) by assumption
   so val(log(K(b))) >= l, but
   val(X) = 1/deg(K) val(norm(X)) for p-adics
   This is true for all completions involved, and sum degrees is n
 =#

  t = prod([a[i]^g[i] for i=1:length(a)])
  # t is either 1 or 1-t is large, norm(1-t) is div. by p^ln
  #in this case T2(1-t) is large, by the arguments above: T2>= (np^l)^2=:B
  # and, see the bottom, \|Log()\|_2^2 >= 1/4 arcosh((B-2)/2)^2
  B = ArbField(nbits(v)*2)(v)^2
  B = 1/2 *acosh((B-2)/2)^2
  p = Hecke.upper_bound(log(B)/log(parent(B)(2)), fmpz)
  b = conjugates_arb_log(t, -Int(p))
  global res = (B, b, t)
  @show (B > sum(x*x for x = b)) == istorsion_unit(t)[1]
  @assert (B > sum(x*x for x = b)) == istorsion_unit(t)[1]
  return B > sum(x*x for x = b)
end

function Base.setprecision(a::padic, N::Int)
  b = deepcopy(a)
  b.N = N
  @assert precision(b) == N
  return b
end

function mult_syzygies_units(A::Array{FacElem{nf_elem, AnticNumberField}, 1})
  p = next_prime(106) #experimentally, the runtime is dominated by
         # log which in case is dominated by the a^(p^n-1) in the 1st step
         # thus try p smaller..., ideally also try n smaller...
         # also, see comments at _log
  p = 10429
  u = FacElem{nf_elem, AnticNumberField}[]
  prec = 10  #needs to match the 10 used in the creation of the
             #ctx.
             #or prec needs to be honored and downwards corrected
  r1, r2 = signature(base_ring(parent(A[1])))
  r = r1+r2 -1
  n = degree(base_ring(parent(A[1])))
  la = conjugates_pAdic_log(A[1], p, prec)
  lu = zero_matrix(base_ring(la), 0, n)
  uu = []
  for a = A
    while true
      la = conjugates_pAdic_log(a, p, prec)
      if iszero(la)
        @assert verify_gamma([a], [fmpz(1)], fmpz(p)^prec)
        println("torsion found")
        break
      end
      lv = vcat(lu, la)
      #check_precision and change
      if false && any(x->precision(x) < prec, lv)
        println("loss of precision - not sure what to do")
        for i=1:rows(lv)
          for j = cols(lv) #seems to not do anything
            lv[i, j] = setprecision(lv[i, j], min_p)
            @assert precision(lv[i,j]) == min_p
          end
        end
      end
      k = Hecke.left_kernel_basis(lv)
      @assert length(k) < 2
      if length(k) == 0
        println("new ")
        push!(u, a)
        lu = vcat(lu, la)
        @assert length(u) <= r
      else # length == 1 extend the module
        @show s = [lift_reco(FlintQQ, x, reco = true) for x = k[1]]
        if any(x -> x === false, s)
          prec *= 2
          @show "increase prec to ", prec
          lu = vcat([conjugates_pAdic_log(x, p, prec) for x = u])
          continue
        end
        d = reduce(lcm, map(denominator, s))
        gamma = [FlintZZ(x*d) for x = s]
        @assert reduce(gcd, gamma) == 1 # should be a primitive relation
        if !verify_gamma(push!(copy(u), a), gamma, fmpz(p)^prec)
          prec *= 2
          @show "increase prec to ", prec
          lu = vcat([conjugates_pAdic_log(x, p, prec) for x = u])
          continue
        end
        @assert length(gamma) == length(u)+1
        gamma = vcat(gamma[1:length(u)], [0 for i=length(u)+1:r+length(uu)], [gamma[end]])
        push!(uu, (a, gamma))
      end
      break
    end
  end
  #=
    let u_1, .., u_n be units and
       <u_i | i> has rank s and
        r_i in Z^n be such that
          prod u_i^r_i = 1  (OK, sum of the logs is zero)
          rank <r_i | i> = s as well
    so the r_i form a Q-basis for the relations.
    Essentially, the gamma of above are the r_i
    Let [H|0] = [r_i|i]*T be the hnf with trafo, so T in Gl(n, Z)
    Then
      <u_i|i> = <[u_i|i] T>
      [r_i|i] * [u_i|i]^t = 0 (by construction)
      [r_i|i] T inv(T) [u[i] | i] = 0
      [H | 0]   [v_i | i] = 0
      so, since H is triangular(!!) v_1, ... v_n-s = 0
      and <u_i |i> = <v_n-s+1, ..., v_n>

    for the case of n=s+1 this is mostly the "normal" construction.
    Note: as a side, the relations do not have to be primitive.
      If they are, (and n=s+1), then H = 1
  =#

  for i=1:length(uu)-1
    append!(uu[i][2], zeros(FlintZZ, length(uu[end][2])-length(uu[i][2])))
  end
  U = matrix(FlintZZ, length(uu), length(uu[end][2]), reduce(vcat, [x[2] for x = uu]))
  _, U = hnf_with_trafo(U')
  if false
    U = inv(U)
    V = sub(U, 1:rows(U), 1:cols(U)-length(u))
    U = sub(U, 1:rows(U), cols(U)-length(u)+1:cols(U))
    #U can be reduced modulo V...
    Z = zero_matrix(FlintZZ, cols(V), cols(U))
    I = identity_matrix(FlintZZ, cols(U)) * p^(2*prec)
    k = base_ring(A[1])
    A = [ Z V'; I U']
    l = lll(A)
    U = sub(l, cols(V)+1:rows(l), cols(U)+1:cols(l))
    U = lll(U)
  else
    U = lll(U')
  end
  return Hecke._transform(vcat(u, FacElem{nf_elem,AnticNumberField}[FacElem(k(1)) for i=length(u)+1:r], [x[1] for x = uu]), U')
end

function non_torsion_lower_bound(R::NfOrd, B::Int = 2*degree(R))
  L = Hecke.enum_ctx_from_ideal(1*R, zero_matrix(FlintZZ, 0, 0))
  n = degree(R)
  i = B
  while true
    #lat enum is wrong. Or maybe rounding funnily:
    # try on wildanger_field(13, 17) and observe the vectors found...
    # in particular 1 is not (reliably) found
    @show s = Hecke.enum_ctx_short_elements(L, i*L.d)
    #remove torsion
    if nrows(s) > 5
      M = minimum([t2(R(collect(sub(s, i:i, 1:n)))) for i=1:nrows(s)])
      j = n-2
      return (n-j)/4*acosh((M-j)/(n-j))^2
    end
    @show i += 1
  end
end

function unit_lower_bound(R::NfOrd, B::Int = 2*degree(R))
  L = Hecke.enum_ctx_from_ideal(1*R, zero_matrix(FlintZZ, 0, 0))
  n = degree(R)
  i = B
  while true
    s = Hecke.enum_ctx_short_elements(L, i*L.d)
    #remove torsion!!! now that \pm 1 is actually found
    e = [R(collect(sub(s, i:i, 1:n))) for i=1:nrows(s)]
    u = [x for x in e if isunit(x)]
    if nrows(s) > 5
      if length(u) == 0
        R = parent(t2(e[1]))
        @show M = R(i)
      else
        @show M = minimum([t2(x) for x = u])
      end
      j = n-2
      return (n-j)/4*acosh((M-j)/(n-j))^2, u
    end
    @show i += 1
  end
end

function regulator_lower_bound(R::NfOrd, B::Int = 2*degree(R))
  Ms, _ = unit_lower_bound(R, B)
  r1, r2 = signature(R)
  r = r1+r2-1
  n = degree(R)
  return Ms^r * 2^r2 /n / Hecke.hermite_constant(r)
end

end


text:
"""
    get_buff_modifier(buff)

Compute the mulitplier associated with stat buffs (multiplied by 12 to return
an integer). As a result, the multiplier for no buff effect is 0. Inputs should
be between -4 and 4.

# Examples
```jldoctest
julia> get_buff_modifier(Int8(0))
12
"""
function get_buff_modifier(i::UInt32, agent::UInt8)
    a, d = agent == 0x01 ?
        (UInt8((i ÷ UInt32(13230))   % 9),
        UInt8((i ÷ UInt32(119070))   % 9)) :
        (UInt8((i ÷ UInt32(1071630)) % 9),
        UInt8((i ÷ UInt32(9644670))  % 9))
    a1, a2 = a > UInt8(3) ? (a, UInt8(4)) :
        (UInt8(4), UInt8(8) - a)
    d1, d2 = d > UInt8(3) ? (d, UInt8(4)) :
        (UInt8(4), UInt8(8) - d)
    return a1*d2, a2*d1
end

evaluate_fast_moves(state::DynamicIndividualState,
    static_state::StaticIndividualState, using_fm::Tuple{Bool, Bool}) =
    DynamicIndividualState(
        add_energy(damage(state[0x01],
            using_fm[0x02] ? calculate_damage(
                static_state[0x02].stats.attack,
                state.data,
                static_state[0x01],
                static_state[0x02].fast_move,
            ) : 0x0000), (using_fm[0x01] ?
            get_energy(static_state[0x01].fast_move) : Int8(0))),
        add_energy(damage(state[0x02],
            using_fm[0x01] ? calculate_damage(
                static_state[0x01].stats.attack,
                state.data,
                static_state[0x02],
                static_state[0x01].fast_move,
            ) : 0x0000), (using_fm[0x02] ?
            get_energy(static_state[0x02].fast_move) : Int8(0))),
        state.data)

"""
    evaluate_charged_move(state, static_state, cmp, move_id, charge,
        shielding, buffs_applied)

Takes in the dynamic state, the static state, the attacking agent, the move,
the charge, whether or not the opponent shields, and whether or not buffs are
applied (say in the case of a random buff move) and returns
the dynamic state after the charged move has occurred, with precisely one copy
"""
function evaluate_charged_move(state::DynamicIndividualState,
    static_state::StaticIndividualState, cmp::UInt32, move_id::UInt8,
    charge::UInt8, shielding::Bool)
    next_state = state
    agent = isodd(cmp) ? 0x01 : 0x02
    d_agent = get_other_agent(agent)
    data = next_state.data
    move = move_id == 0x01 ? static_state[agent].charged_move_1 :
        static_state[agent].charged_move_2

    buff_chance = move.buffChance
    if buff_chance == Int8(100)
        data = apply_buff(data, move, agent)
    elseif buff_chance != Int8(0)
        data += agent == 0x01 ? move_id == 0x01 ? UInt32(2205) : UInt32(4410) :
                             move_id == 0x01 ? UInt32(6615) : UInt32(8820)
    end

    attacking_team = subtract_energy(next_state[agent], get_energy(move))

    if shielding
        defending_team = damage(next_state[d_agent], 0x0001)
        data -= d_agent == 0x01 ? 1 : 3
    else
        defending_team = damage(next_state[d_agent], calculate_damage(
                static_state[agent].stats.attack,
                state.data,
                static_state[d_agent],
                move,
                charge = charge
            ))
    end
    return DynamicIndividualState(
        agent == 1 ? attacking_team : defending_team,
        agent == 2 ? attacking_team : defending_team,
        # go from cmp 4 (2 then 1) to cmp 1
        # or go from cmp 3 (1 then 2) to cmp 2
        # or go from cmp 2 to 0
        # or go from cmp 1 to 0
        data  - (cmp == 4 ? (get_hp(defending_team) != 0x0000 ?
                    1323 : 1764) :
                (cmp == 3 ? (get_hp(defending_team) != 0x0000 ?
                    441 : 1323) :
                (cmp == 2 ? 882 : 441)))
    )
end

function apply_buff(data::UInt32, move::ChargedMove, agent::UInt8)
    buffs1 = get_buffs(data, 1)
    buffs2 = get_buffs(data, 2)
     return agent == 1 ? data + Int32(13230) * clamp(
        get_atk(move.self_buffs), -Int8(buffs1[1]), 9 - Int8(buffs1[1])) +
        Int32(119070) * clamp(get_def(move.opp_buffs),
        -Int8(buffs1[2]), 9 - Int8(buffs1[2]))  +
        Int32(1071630) * clamp(get_atk(move.opp_buffs),
        -Int8(buffs2[1]), 9 - Int8(buffs2[1]))  +
        Int32(9644670) * clamp(get_def(move.self_buffs),
        -Int8(buffs2[2]), 9 - Int8(buffs2[2])) :
        data + Int32(13230) * clamp(get_atk(move.opp_buffs),
        -Int8(buffs1[1]), 9 - Int8(buffs1[1])) +
        Int32(119070) * clamp(get_def(move.self_buffs),
        -Int8(buffs1[2]), 9 - Int8(buffs1[2]))  +
        Int32(1071630) * clamp(get_atk(move.self_buffs),
        -Int8(buffs2[1]), 9 - Int8(buffs2[1]))  +
        Int32(9644670) * clamp(get_def(move.opp_buffs),
        -Int8(buffs2[2]), 9 - Int8(buffs2[2]))
end

function step_timers(state::DynamicIndividualState, fmCooldown1::Int8,
    fmCooldown2::Int8)
    fmPending = get_fast_moves_pending(state)
    data = state.data
    if fmCooldown1 != Int8(0)
        data += (UInt32(fmCooldown1) - fmPending[1]) * 9
    elseif fmPending[1] != UInt32(0)
        data -= 9
    end
    if fmCooldown2 != Int8(0)
        data += (UInt32(fmCooldown2) - fmPending[2]) * 63
    elseif fmPending[2] != UInt32(0)
        data -= 63
    end

    return DynamicIndividualState(state[0x01], state[0x02], data)
end

"""
    min_score(state, static_state)

Given the state and the static state (here just for starting hp values), compute
the PvPoke-like score that would occur if the first agent stopped attacking
altogether. This is currently only used in computing the final score, but it
could be used as strict bounds for α/β pruning, for example.
"""
function min_score(state::DynamicIndividualState,
    static_state::StaticIndividualState)
     return 0.5 * (static_state[0x02].stats.hitpoints - get_hp(state[0x02])) /
        static_state[0x02].stats.hitpoints
end

"""
    max_score(state, static_state)

Given the state and the static state (here just for starting hp values), compute
the PvPoke-like score that would occur if the second agent stopped attacking
altogether. This is currently only used in computing the final score, but it
could be used as strict bounds for α/β pruning, for example.
"""
function max_score(state::DynamicIndividualState,
    static_state::StaticIndividualState)
     return 0.5 + 0.5 * get_hp(state[0x01]) / static_state[0x01].stats.hitpoints
end

"""
    battle_score(state, static_state)

Given the state and the static state (here just for starting hp values), compute
the PvPoke-like score for the battle. Note that this can also be computed for
battles in progress, and thus differs from PvPoke's use cases
"""
function battle_score(state::DynamicIndividualState,
    static_state::StaticIndividualState)
    return min_score(state, static_state) + max_score(state, static_state) - 0.5
end

ext: jl


text:
# a simple LP example
#
#   maximize x + y
#
#   s.t. 50 x + 24 y <= 2400
#        30 x + 33 y <= 2100
#        x >= 45, y >= 5
#
#   solution: x = 45, y = 6.25, objv = 51.25

using CPLEX

env = CPLEX.Env()

# method = getparam(env, "Method")
# println("method = $method")

model = CPLEX.Model(env, "lp_01")
CPLEX.set_sense!(model, :Max)

# add variables
CPLEX.add_var!(model, 1.0, 45., Inf)  # x
CPLEX.add_var!(model, 1.0,  5., Inf)  # y

# add constraints
CPLEX.add_constr!(model, [50., 24.], '<', 2400.)
CPLEX.add_constr!(model, [30., 33.], '<', 2100.)

println(model)

# perform optimization
CPLEX.optimize!(model)

sol = CPLEX.get_solution(model)
println("soln = $(sol)")

objv = CPLEX.get_objval(model)
println("objv = $(objv)")

gc()  # test finalizers


text:
#Implements Future Value calculation as per Excel
#Type = 1 means payment at start of period
function FV(Rate=0,Nper=0,Pmt=0,Pv=0,Type=0)
    Fv = (-1*Pv*(1 + Rate)^Nper - Pmt*(1 + Rate*Type)*(((1 + Rate)^Nper - 1)/Rate))
	round(Fv, 2)
end

#######################################

#Implements Present Value calculation as per Excel
#Type = 1 means payment at start of period
function PV(Rate=0,Nper=0,Pmt=0,Fv=0,Type=0)
    Pv = (-1*Fv - Pmt*(1 + Rate*Type)*(((1 + Rate)^Nper - 1)/Rate))/(1 + Rate)^Nper
	round(Pv, 2)
end

#######################################

#Implements Payment calculation as per Excel
#Type = 1 means payment at start of period
function PMT(Rate=0,Nper=0,Pv=0,Fv=0,Type=0)
    Pmt = (-1*Fv - Pv*(1 + Rate)^Nper)/((1 + Rate*Type)*((1 + Rate)^Nper - 1)/Rate)
	round(Pmt, 2)
end

#######################################

#Implements Effective Annual Rate calculation as per Excel
function EFFECT(Nominal_rate,Npery)
	Effect = ((1 + Nominal_rate/Npery)^Npery - 1)
	round(Effect, 6)
end

#######################################

#Implements Net Present Value calculation as per Excel
#Add the option of adding inital payment
function NPV(Rate,Cash_Flow,InitialPayment=0)
	npv=0
	for i in 1:length(Cash_Flow)
		npv += Cash_Flow[i] / (1+Rate)^i
	end
	npv += InitialPayment
	round(npv, 2)
end

#######################################

#Implements Internal Rate of Return calculation as per Excel
function IRR(Cash_Flow,Guess=0.1)
	const limit=0.0000001
	initial=Cash_Flow[1]
	values=Cash_Flow[2:length(Cash_Flow)]

	rate=Guess*1.01
	oldRate=Guess

	while (abs((rate-oldRate)/rate) > limit)
		newRate=rate-NPV(rate,values,initial)*((rate-oldRate)/(NPV(rate,values,initial)-NPV(oldRate,values,initial)))
		oldRate=rate
		rate=newRate
	end

	round(rate,6)
end

max_stars_repo_name: arnaudamiel/Finance

max_issues_repo_name: arnaudamiel/Finance

max_forks_repo_name: arnaudamiel/Finance


text:
# Tests for the Arbiter package

using Arbiter
using FactCheck

include("test_arbiter.jl")
include("test_graph.jl")
include("test_scheduler.jl")
include("test_sync.jl")

FactCheck.exitstatus()


text:
@doc doc"""
    Sphere{N} <: Manifold

The unit sphere manifold $\mathbb S^n$ represented by $n+1$-Tuples, i.e. in by
vectors in $\mathbb R^{n+1}$ of unit length

# Constructor

    Sphere(n)

Generate the $\mathbb S^{n}\subset \mathbb R^{n+1}$
"""
struct Sphere{N} <: Manifold end

Sphere(n::Int) = Sphere{n}()

function get_basis(M::Sphere{N}, x, B::DiagonalizingOrthonormalBasis) where {N}
    A = zeros(N + 1, N + 1)
    A[1, :] = transpose(x)
    A[2, :] = transpose(B.v)
    V = nullspace(A)
    κ = ones(N)
    if !iszero(B.v)
        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
        V = cat(B.v / norm(M, x, B.v), V; dims = 2)
        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
    end
    vecs = [V[:, i] for i = 1:N]
    return PrecomputedDiagonalizingOrthonormalBasis(vecs, κ)
end

"""
    check_manifold_point(S, x; kwargs...)

Check whether `x` is a valid point on the [`Sphere`](@ref) `S`, i.e. is a vector
of length [`manifold_dimension`](@ref)`(S)+1` (approximately) of unit length.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(S::Sphere{N}, x; kwargs...) where {N}
    if size(x) != representation_size(S)
        return DomainError(
            size(x),
            "The point $(x) does not lie on $S, since its size is not $(N+1).",
        )
    end
    if !isapprox(norm(x), 1.0; kwargs...)
        return DomainError(
            norm(x),
            "The point $(x) does not lie on the sphere $(S) since its norm is not 1.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(S, x, v; kwargs... )

Check whether `v` is a tangent vector to `x` on the [`Sphere`](@ref) `S`, i.e.
after [`check_manifold_point`](@ref)`(S,x)`, `v` has to be of same dimension as `x`
and orthogonal to `x`.
The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(S::Sphere{N}, x, v; kwargs...) where {N}
    perr = check_manifold_point(S, x)
    perr === nothing || return perr
    if size(v) != representation_size(S)
        return DomainError(
            size(v),
            "The vector $(v) is not a tangent to a point on $S since its size does not match $(N+1).",
        )
    end
    if !isapprox(abs(dot(x, v)), 0.0; kwargs...)
        return DomainError(
            abs(dot(x, v)),
            "The vector $(v) is not a tangent vector to $(x) on $(S), since it is not orthogonal in the embedding.",
        )
    end
    return nothing
end

@doc doc"""
    distance(M::Sphere, x, y)

Compute the geodesic distance betweeen `x` and `y` on the [`Sphere`](@ref) `M`.
The formula is given by the (shorter) great arc length on the (or a) great circle
both `x` and `y` lie on.

````math
d_{\mathbb S^n}(x,y) = \operatorname{acos}(\langle x, y\rangle).
````
"""
distance(S::Sphere, x, y) = acos(clamp(dot(x, y), -1, 1))

@doc doc"""
    exp(M::Sphere, x, v)

Compute the exponential map from `x` into the tangent direction `v` on the [`Sphere`](@ref)
`M` by following the great arc eminating from `x` in direction `v` along with length of `v`.

````math
\exp_x v = \cos(\lVert v \rVert_x)x + \sin(\lVert v \rVert_x)\frac{v}{\lVert v \rVert_x},
````
where $\lVert v \rVert_x$ is the [`norm`](@ref norm(::Sphere,x,v)) on the
[`Sphere`](@ref) `M`.
"""
exp(::Sphere, ::Any...)

function exp!(M::Sphere, y, x, v)
    θ = norm(M, x, v)
    y .= cos(θ) .* x .+ usinc(θ) .* v
    return y
end

flat!(M::Sphere, v::CoTFVector, x, w::TFVector) = copyto!(v, w)

@doc doc"""
    injectivity_radius(M::Sphere[, x])

Return the injectivity radius for the [`Sphere`](@ref) `M`, which is globally $\pi$.

    injectivity_radius(M::Sphere, x, ::ProjectionRetraction)

Return the injectivity radius for the [`ProjectionRetraction`](@ref) on the
[`Sphere`](@ref), which is globally $\frac{\pi}{2}$.
"""
injectivity_radius(::Sphere, ::Any...) = π
injectivity_radius(::Sphere, ::Any, ::ProjectionRetraction) = π / 2

@doc doc"""
    inner(S::Sphere, x, w, v)

Compute the inner product of the two tangent vectors `w,v` from the tangent
plane at `x` on the sphere `S=`$\mathbb S^n$ using the restriction of the
metric from the embedding, i.e. $ (v,w)_x = v^\mathrm{T}w $.
"""
@inline inner(S::Sphere, x, w, v) = dot(w, v)

function get_vector(M::Sphere{N}, x, v, B::ArbitraryOrthonormalBasis) where {N}
    x[1] ≈ 1 && return vcat(0, v)
    xp1 = x .+ ntuple(i -> ifelse(i == 1, 1, 0), N + 1)
    v0 = vcat(0, v)
    return 2 * xp1 * dot(xp1, v0) / dot(xp1, xp1) - v0
end

@doc doc"""
    inverse_retract(M::Sphere, x, y, ::ProjectionInverseRetraction)

Compute the inverse of the projection based retraction on the [`Sphere`](@ref),
i.e. rearranging $x+v = y\lVert x+d \rVert_2$ yields
since $\langle x,v\rangle = 0$ and when $d_{\mathbb S^2}(x,y) \leq \frac{\pi}{2}$ that

````math
\operatorname{retr}_x^{-1}(y) = \frac{y}{\langle x, y \rangle} - x.
````
"""
inverse_retract(::Sphere, ::Any, ::Any, ::ProjectionInverseRetraction)

function inverse_retract!(::Sphere, v, x, y, ::ProjectionInverseRetraction)
    return (v .= y ./ dot(x, y) .- x)
end

@doc doc"""
    log(M::Sphere, x, y)

Compute the logarithmic map on the [`Sphere`](@ref) `M`, i.e. the tangent vector, whose
geodesic starting from `x` reaches `y` after time 1.
The formula reads for $x\neq -y$

````math
\log_x y = d_{\mathbb S^n}(x,y)
\frac{y-\langle x,y\rangle x}{\lVert y-\langle x,y\rangle x \rVert_2},
````
and a deterministic choice from the set of tangent vectors is returned if $x=-y$, i.e. for
opposite points.
"""
log(::Sphere, ::Any...)

function log!(S::Sphere, v, x, y)
    cosθ = dot(x, y)
    if cosθ ≈ -1 # appr. opposing points, return deterministic choice from set-valued log
        fill!(v, 0)
        if x[1] ≈ 1
            v[2] = 1
        else
            v[1] = 1
        end
        copyto!(v, v .- dot(x, v) .* x)
        v .*= π / norm(v)
    else
        cosθ = cosθ > 1 ? one(cosθ) : cosθ
        θ = acos(cosθ)
        v .= (y .- cosθ .* x) ./ usinc(θ)
    end
    return project_tangent!(S, v, x, v)
end

@doc doc"""
    manifold_dimension(S::Sphere)

Return the dimension of the manifold $\mathbb S^n$, i.e. $n$.
"""
manifold_dimension(S::Sphere{N}) where {N} = N

"""
    mean(
        S::Sphere,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolationWithinRadius(π/2);
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` using
[`GeodesicInterpolationWithinRadius`](@ref).
"""
mean(::Sphere, ::Any...)

function mean!(S::Sphere, y, x::AbstractVector, w::AbstractVector; kwargs...)
    return mean!(S, y, x, w, GeodesicInterpolationWithinRadius(π / 2); kwargs...)
end

@doc doc"""
    norm(M::Sphere, x, v)

Compute the length of the tangent vector `v` from the tangent space at `x` on the
[`Sphere`](@ref) `M`, which is the norm in the embedding, i.e.
````math
\lVert v \rVert_x = \lVert v \rVert_2.
````
"""
norm(M::Sphere, x, v) = norm(v)

"""
    normal_tvector_distribution(S::Sphere, x, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `x`.
"""
function normal_tvector_distribution(S::Sphere, x, σ)
    d = Distributions.MvNormal(zero(x), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(S), x, d, project_vector!, x)
end

@doc doc"""
    project_point(M::Sphere, x)

Project the point `x` from the embedding onto the [`Sphere`](@ref) `M`.

````math
\operatorname{proj}_{\mathbb S^n}(x) = \frac{x}{\lVert x \rVert_2}.
````
"""
project_point(::Sphere, ::Any...)

project_point!(S::Sphere, x) = (x ./= norm(x))

@doc doc"""
    project_tangent(M::Sphere, x, v)

Project the point `v` onto the tangent space at `x` on the [`Sphere`](@ref) `M`.

````math
\operatorname{proj}_{x}(v) = v - \langle x, v \rangle_x
````
"""
project_tangent(::Sphere, ::Any...)

project_tangent!(S::Sphere, w, x, v) = (w .= v .- dot(x, v) .* x)

@doc doc"""
    get_coordinates(M::Sphere, x, v, B::ArbitraryOrthonormalBasis)

Represent the tangent vector `v` at point `x` from a sphere `M` in
an orthonormal basis by rotating the vector `v` using rotation matrix
$2\frac{x_p x_p^\mathrm{T}}{x_p^\mathrm{T} x_p} - I$ where $x_p = x + (1, 0, \dots, 0)$.
"""
function get_coordinates(M::Sphere{N}, x, v, B::ArbitraryOrthonormalBasis) where {N}
    if isapprox(x[1], 1)
        return v[2:end]
    else
        xp1 = x .+ ntuple(i -> ifelse(i == 1, 1, 0), N + 1)
        return (2*xp1*dot(xp1, v)/dot(xp1, xp1)-v)[2:end]
    end
end

@doc doc"""
    representation_size(M::Sphere)

Return the size points on the [`Sphere`](@ref) `M` are represented as, i.e.
for the `n`-dimensional [`Sphere`](@ref) it is vectors of size `(n+1,)`.
"""
@generated representation_size(::Sphere{N}) where {N} = (N + 1,)

@doc doc"""
    retract(M::Sphere, x, y, ::ProjectionRetraction)

Compute the retraction that is based on projection, i.e.

````math
\operatorname{retr}_x(v) = \frac{x+v}{\lVert x+v \rVert_2}
````
"""
retract(::Sphere, ::Any, ::Any, ::ProjectionRetraction)

function retract!(M::Sphere, y, x, v, ::ProjectionRetraction)
    y .= x .+ v
    return project_point!(M, y)
end

sharp!(M::Sphere, v::TFVector, x, w::CoTFVector) = copyto!(v, w)

"""
    uniform_distribution(S::Sphere, x)

Uniform distribution on given sphere. Generated points will be of similar
type to `x`.
"""
function uniform_distribution(S::Sphere, x)
    d = Distributions.MvNormal(zero(x), 1.0)
    return ProjectedPointDistribution(S, d, project_point!, x)
end

@doc doc"""
    vector_transport_to(M::Sphere, x, v, y, ::ParallelTransport)

Compute the [`ParallelTransport`](@ref) on the [`Sphere`](@ref) `M`, which is given by

````math
P_{y\gets x}(v) = v - \frac{\langle \log_xy,v\rangle_x}{d^2_{\mathbb S^n}(x,y)}
\bigl(\log_xy + \log_yx \bigr).
````
"""
vector_transport_to(::Sphere, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(M::Sphere, vto, x, v, y, ::ParallelTransport)
    v_xy = log(M, x, y)
    vl = norm(M, x, v_xy)
    copyto!(vto, v)
    if vl > 0
        factor = 2 * dot(v, y) / (norm(x + y)^2)
        vto .-= factor .* (x .+ y)
    end
    return vto
end

@doc doc"""
    zero_tangent_vector(M::Sphere, x)

Return the zero tangent vector from the tangent space at `x` on the [`Sphere`](@ref) `M`,
which is the zero vector in the embedding.
"""
zero_tangent_vector(::Sphere, ::Any...)

zero_tangent_vector!(S::Sphere, v, x) = fill!(v, 0)


text:
## Add Linear constraints

# add single constraint

function add_constr!(model::Model, inds::IVec, coeffs::FVec, rel::Cchar, rhs::Float64)
    if rel == convert(Cchar,'>')
        rel = XPRS_GEQ
    elseif rel == convert(Cchar,'<')
        rel = XPRS_LEQ
    elseif rel == convert(Cchar,'=')
        rel = XPRS_EQ
    end

    length(inds) == length(coeffs) || error("Inconsistent argument dimensions.")

    ret = @xprs_ccall(addrows, Cint,(
        Ptr{Void}, #prob
        Cint,    #number new rowws
        Cint, #number nonzeros
        Ptr{Cchar}, #cstr type
        Ptr{Float64}, # rhs
        Ptr{Float64}, # range (size new row)
        Ptr{Cint}, # start (size newrow)
        Ptr{Cint}, # ind (size new nz)
        Ptr{Float64} # val (size newnz)
        ),
        model.ptr_model, 1, length(inds), Cchar[rel], Float64[rhs], C_NULL, Cint[0], inds-1, coeffs
    )

    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

function add_constr!(model::Model, inds::Vector, coeffs::Vector, rel::GChars, rhs::Real)
    add_constr!(model, ivec(inds), fvec(coeffs), cchar(rel), Float64(rhs))
end

function add_constr!(model::Model, coeffs::Vector, rel::GChars, rhs::Real)
    inds = find(coeffs)
    vals = coeffs[inds]
    add_constr!(model, inds, vals, rel, rhs)
end

# add multiple constraints

function add_constrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, rel::CVec, rhs::FVec)

    for i in 1:length(rel)
        if rel[i] == convert(Cchar,'>')
            rel[i] = XPRS_GEQ
        elseif rel[i] == convert(Cchar,'<')
            rel[i] = XPRS_LEQ
        elseif rel[i] == convert(Cchar,'=')
            rel[i] = XPRS_EQ
        end
    end

    m = length(cbegins)
    nnz = length(inds)
    (m == length(rel) == length(rhs) && nnz == length(coeffs)) || error("Inconsistent argument dimensions.")

    if m > 0
        ret = @xprs_ccall(addrows, Cint,(
            Ptr{Void}, #prob
            Cint,    #number new rowws
            Cint, #number nonzeros
            Ptr{Cchar}, #cstr type
            Ptr{Float64}, # rhs
            Ptr{Float64}, # range (size new row)
            Ptr{Cint}, # start (size newrow)
            Ptr{Cint}, # ind (size new nz)
            Ptr{Float64} # val (size newnz)
            ),
            model.ptr_model, m, nnz, rel, rhs, C_NULL, cbegins - 1, inds - 1, coeffs
        )

        if ret != 0
            throw(XpressError(model))
        end
    end
    nothing
end

function add_constrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, rel::GCharOrVec, rhs::Vector)
    add_constrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), cvecx(rel, length(cbeg)), fvec(rhs))
end

function add_constrs_t!(model::Model, At::SparseMatrixCSC{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs!(model, At.colptr[1:At.n], At.rowval, At.nzval, rel, b)
end

function add_constrs_t!(model::Model, At::Matrix{Float64}, rel::GCharOrVec, b::Vector)
    n, m = size(At)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, sparse(At), rel, b)
end

function add_constrs!(model::Model, A::CoeffMat, rel::GCharOrVec, b::Vector{Float64})
    m, n = size(A)
    (m == length(b) && n == num_vars(model)) || error("Incompatible argument dimensions.")
    add_constrs_t!(model, transpose(A), rel, b)
end

# add single range constraint

function add_rangeconstr!(model::Model, inds::IVec, coeffs::FVec, lb::Float64, ub::Float64)
    # b = ub
    r = ub - lb
    ret = @xprs_ccall(addrows, Cint,(
        Ptr{Void}, #prob
        Cint,    #number new rowws
        Cint, #number nonzeros
        Ptr{Cchar}, #cstr type
        Ptr{Float64}, # rhs
        Ptr{Float64}, # range (size new row)
        Ptr{Cint}, # start (size newrow)
        Ptr{Cint}, # ind (size new nz)
        Ptr{Float64} # val (size newnz)
        ),
        model.ptr_model, 1, length(inds), Cchar['R'], Float64[ub], Float64[r], Cint[0], inds-1, coeffs
    )
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

function add_rangeconstr!(model::Model, inds::Vector, coeffs::Vector, lb::Real, ub::Real)
    add_rangeconstr!(model, ivec(inds), fvec(coeffs), Float64(lb), Float64(ub))
end

function add_rangeconstr!(model::Model, coeffs::Vector, lb::Real, ub::Real)
    inds = find(coeffs)
    vals = coeffs[inds]
    add_rangeconstr!(model, inds, vals, lb, ub)
end

# add multiple range constraints

function add_rangeconstrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, lb::FVec, ub::FVec)
    m = length(cbegins)
    nnz = length(inds)
    (m == length(lb) == length(ub) && nnz == length(coeffs)) || error("Incompatible argument dimensions.")

    if m > 0
        r = ub - lb
        ret = @xprs_ccall(addrows, Cint,(
            Ptr{Void}, #prob
            Cint,    #number new rowws
            Cint, #number nonzeros
            Ptr{Cchar}, #cstr type
            Ptr{Float64}, # rhs
            Ptr{Float64}, # range (size new row)
            Ptr{Cint}, # start (size newrow)
            Ptr{Cint}, # ind (size new nz)
            Ptr{Float64} # val (size newnz)
            ),
            model.ptr_model, m, nnz, cvecx( convert(Cchar,'R') , m), ub, r, cbegins-1, inds-1, coeffs
        )

        if ret != 0
            throw(XpressError(model))
        end
    end
    nothing
end

function add_rangeconstrs!(model::Model, cbeg::Vector, inds::Vector, coeffs::Vector, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, ivec(cbeg), ivec(inds), fvec(coeffs), fvec(lb), fvec(ub))
end

function add_rangeconstrs_t!(model::Model, At::SparseMatrixCSC{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs!(model, At.colptr[1:At.n], At.rowval, At.nzval, lb, ub)
end

function add_rangeconstrs_t!(model::Model, At::Matrix{Float64}, lb::Vector, ub::Vector)
    add_rangeconstrs_t!(model, sparse(At), lb, ub)
end

function add_rangeconstrs!(model::Model, A::CoeffMat, lb::Vector, ub::Vector)
    m, n = size(A)
    (m == length(lb) == length(ub) && n == num_vars(model)) || error("Incompatible argument dimensions.")

    add_rangeconstrs_t!(model, transpose(A), lb, ub)
end

# prepare macros for retrieving numvars and num cols and num NNZ

function get_constrmatrix(model::Model)
    nnz = num_cnzs(model)
    m = num_constrs(model)
    n = num_vars(model)
    numnzP = Array(Cint, 1)
    cbeg = Array(Cint, m+1)
    cind = Array(Cint, nnz)
    cval = Array(Float64, nnz)
    ret = @xprs_ccall(getrows, Cint, (
                     Ptr{Void},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64},
                     Cint,
                     Ptr{Cint},
                     Cint,
                     Cint
                     ),
                     model.ptr_model,
                     cbeg,
                     cind,
                     cval,
                     nnz,
                     numnzP,
                     0,
                     m-1)
    if ret != 0
        throw(XpressError(model))
    end
    cbeg[end] = nnz
    I = Array(Int, nnz)
    J = Array(Int, nnz)
    V = Array(Float64, nnz)
    for i in 1:length(cbeg)-1
        for j in (cbeg[i]+1):cbeg[i+1]
            I[j] = i
            J[j] = cind[j]+1
            V[j] = cval[j]
        end
    end
    return sparse(I, J, V, m, n)
end

function add_sos!(model::Model, sostype::Symbol, idx::Vector{Int}, weight::Vector{Cdouble})
    ((nelem = length(idx)) == length(weight)) || error("Index and weight vectors of unequal length")
    (sostype == :SOS1) ? (typ = XPRS_SOS_TYPE1) : ( (sostype == :SOS2) ? (typ = XPRS_SOS_TYPE2) : error("Invalid SOS constraint type") )
    ret = @xprs_ccall(addsets, Cint, (
                     Ptr{Void},
                     Cint,
                     Cint,
                     Ptr{Cchar},
                     Ptr{Cint},
                     Ptr{Cint},
                     Ptr{Float64}
                     ),
                     model.ptr_model,
                     convert(Cint, 1),
                     convert(Cint, nelem),
                     Cchar[typ],
                     Cint[0,0],
                     convert(Vector{Cint}, idx-1),
                     weight)
    if ret != 0
        throw(XpressError(model))
    end
end

del_constrs!{T<:Real}(model::Model, idx::T) = del_constrs!(model, Cint[idx])
del_constrs!{T<:Real}(model::Model, idx::Vector{T}) = del_constrs!(model, convert(Vector{Cint},idx))
function del_constrs!(model::Model, idx::Vector{Cint})
    numdel = length(idx)
    ret = @xprs_ccall(delrows, Cint, (
                     Ptr{Void},
                     Cint,
                     Ptr{Cint}),
                     model.ptr_model, convert(Cint,numdel), idx.-1)
    if ret != 0
        throw(XpressError(model))
    end
end


text:
using Documenter
using DocumenterCitations
using NQCBase, NQCModels, NQCDistributions, NQCDynamics
using CubeLDFAModel, NNInterfaces

DocMeta.setdocmeta!(NQCDynamics, :DocTestSetup, :(using NQCDynamics); recursive=true)
DocMeta.setdocmeta!(NQCModels, :DocTestSetup, :(using NQCModels, Symbolics); recursive=true)
DocMeta.setdocmeta!(NQCBase, :DocTestSetup, :(using NQCBase); recursive=true)

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"), sorting=:nyt)

function find_all_files(directory)
    map(
        s -> joinpath(directory, s),
        sort(readdir(joinpath(@__DIR__, "src", directory)))
    )
end

@time makedocs(
    bib,
    sitename = "NQCDynamics.jl",
    modules = [NQCDynamics, NQCDistributions, NQCModels, NQCBase, CubeLDFAModel],
    strict = false,
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        canonical = "https://nqcd.github.io/NQCDynamics.jl/stable/",
        assets = ["assets/favicon.ico"],
        ansicolor = true,
        ),
    authors = "James Gardner and contributors.",
    pages = [
        "Introduction" => "index.md"
        "Getting started" => "getting_started.md"
        "Atoms" => "atoms.md"
        "Ensemble simulations" => "ensemble_simulations.md"
        "Saving and loading" => "saving_loading.md"
        "NQCModels.jl" => Any[
            "NQCModels/overview.md"
            "NQCModels/analyticmodels.md"
            "NQCModels/ase.md"
            "NQCModels/neuralnetworkmodels.md"
            "NQCModels/frictionmodels.md"
        ]
        "NQCDistributions.jl" => Any[
            "NQCDistributions/overview.md"
        ]
        "Initial conditions" => find_all_files("initialconditions")
        "Dynamics simulations" => Any[
            "dynamicssimulations/dynamicssimulations.md"
            find_all_files("dynamicssimulations/dynamicsmethods")
        ]
        "Examples" => find_all_files("examples")
        "Developer documentation" => find_all_files("devdocs")
        "API" => Any[
            "NQCBase" => find_all_files("api/NQCBase")
            "NQCModels" => find_all_files("api/NQCModels")
            "NQCDistributions" => find_all_files("api/NQCDistributions")
            "NQCDynamics" => find_all_files("api/NQCDynamics")
        ]
        "References" => "references.md"
    ])

if get(ENV, "CI", nothing) == "true"
    deploydocs(
        repo = "github.com/NQCD/NQCDynamics.jl",
        push_preview=true
    )
end


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using Surrogates
using DataFrames
using RCall
using Random
using Plots;

R"library(hmer)";

Random.seed!(123)
R"set.seed(123)";

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,γ) = p
    @inbounds begin
        du[1] = -β*S*I
        du[2] = β*S*I - γ*I
        du[3] = γ*I
    end
    nothing
end;

tmax = 10000.0
tspan = (0.0,tmax)
δt = 1.0;

n_train = 50 # Number of training samples
n_test = 1000; # Number of test samples

# Parameters are β, γ
lb = [0.00005, 0.1]
ub = [0.001, 1.0];

N = 1000.0
u0 = [990.0,10.0,0.0]
p = [0.0005,0.25]
prob_ode = ODEProblem(sir_ode!,u0,tspan,p);

sampler = LatinHypercubeSample();

θ = Surrogates.sample(n_train,lb,ub,sampler);

logit = (x) -> log(x/(1-x))
invlogit = (x) -> exp(x)/(exp(x)+1.0)
cb_ss = TerminateSteadyState()
logit_final_size = function(z)
  prob = remake(prob_ode;p=z)
  sol = solve(prob, ROS34PW3(),callback=cb_ss)
  fsp = sol[end][3]/N
  logit(fsp)
end;

lfs = logit_final_size.(θ);

training_df = DataFrame(θ)
rename!(training_df,["b","g"])
training_df[!,:lfs] = lfs;

input_ranges = Dict(:b => [lb[1], ub[1]], :g => [lb[2], ub[2]])
output_names = ["lfs"];

@rput training_df
@rput output_names
@rput input_ranges;

R"emulator <- emulator_from_data(training_df, output_names, input_ranges)"

θ_test = sample(n_test,lb,ub,sampler)
lfs_test = logit_final_size.(θ_test);

test_df = DataFrame(θ_test)
rename!(test_df,["b","g"])
test_df[!,:lfs] = lfs_test
@rput test_df; # copy to R

R"lfs_test_pred <- list()"
R"lfs_test_pred$mean <- emulator$lfs$get_exp(test_df)"
R"lfs_test_pred$unc <- emulator$lfs$get_cov(test_df)"
@rget lfs_test_pred;

scatter(invlogit.(lfs_test),
        invlogit.(lfs_test_pred[:mean]),
        xlabel = "Model final size",
        ylabel = "Surrogate final size",
        legend = false,
        title = "Test set")

β_grid = collect(lb[1]:0.00001:ub[1])
θ_eval = [(βᵢ,0.25) for βᵢ in β_grid]
eval_df = DataFrame(θ_eval)
rename!(eval_df,["b","g"])
@rput eval_df
R"lfs_eval <- list()"
R"lfs_eval$mean <- emulator$lfs$get_exp(eval_df)"
R"lfs_eval$unc <- emulator$lfs$get_cov(eval_df)"
@rget lfs_eval
fs_eval = invlogit.(lfs_eval[:mean])
fs_eval_uc = invlogit.(lfs_eval[:mean] .+ 1.96 .* sqrt.(lfs_eval[:unc]))
fs_eval_lc = invlogit.(lfs_eval[:mean] .- 1.96 .* sqrt.(lfs_eval[:unc]))
plot(β_grid,
     fs_eval,
     xlabel = "Infectivity parameter, β",
     ylabel = "Final size",
     label = "Model")
plot!(β_grid,
      invlogit.(logit_final_size.(θ_eval)),
      ribbon = (fs_eval .- fs_eval_lc, fs_eval_uc - fs_eval),
      label = "Surrogate",
      legend = :right)

γ_grid = collect(lb[2]:0.001:ub[2])
θ_eval = [(0.001,γᵢ) for γᵢ in γ_grid]
eval_df = DataFrame(θ_eval)
rename!(eval_df,["b","g"])
@rput eval_df
R"lfs_eval <- list()"
R"lfs_eval$mean <- emulator$lfs$get_exp(eval_df)"
R"lfs_eval$unc <- emulator$lfs$get_cov(eval_df)"
@rget lfs_eval
fs_eval = invlogit.(lfs_eval[:mean])
fs_eval_uc = invlogit.(lfs_eval[:mean] .+ 1.96 .* sqrt.(lfs_eval[:unc]))
fs_eval_lc = invlogit.(lfs_eval[:mean] .- 1.96 .* sqrt.(lfs_eval[:unc]))
plot(γ_grid,
     fs_eval,
     xlabel = "Recovery rate, γ",
     ylabel = "Final size",
     label = "Model")
plot!(γ_grid,
      invlogit.(logit_final_size.(θ_eval)),
      ribbon = (fs_eval .- fs_eval_lc, fs_eval_uc - fs_eval),
      label = "Surrogate")

obs = logit_final_size(p)
obs_lc = logit(0.99*invlogit(obs))
obs_uc = logit(1.01*invlogit(obs))
target = Dict(:lfs => [obs_lc, obs_uc])
@rput target

R"new_points <- generate_new_runs(emulator, $n_test, target, method = 'lhs', cutoff = 3)"
@rget new_points;

l = @layout [a b]
pl1 = histogram(new_points[!,:b],legend=false,xlim=(lb[1],ub[1]),bins=lb[1]:0.00005:ub[1],title="NROY values for β")
vline!(pl1,[p[1]])
pl2 = histogram(new_points[!,:g],legend=false,xlim=(lb[2],ub[2]),bins=lb[2]:0.05:ub[2],title="NROY values for γ")
vline!(pl2,[p[2]])
plot(pl1, pl2, layout = l)


text:
using LinearAlgebra
using SparseArrays
using Random

Random.seed!(0)

## Test IndGraph

m, n = (50, 100)
sm = n + div(n, 10)

function test_against_IndAffine(f, A, cd)
    m, n = size(A)
    if m > n
        return 0.0  # no variant for skinny case for now
    end

    T = typeof(A[1,1])

    ## TODO: remove when IndAffine will be revised
    if T <: Complex{Float64}
        return 0.0
    end

    B = ifelse(issparse(A), [A -SparseMatrixCSC(I, m, m)],  [A -Matrix(I, m, m)])
    # INIT IndAffine for the case
    faff = IndAffine(B, zeros(T, m))

    xy_aff, fx_aff = prox(faff, cd)
    @test faff(xy_aff) == 0.0

    # INIT testing function
    xy, fx = prox(f, cd)
    @test f(xy) == 0.0

    # test against IndAFfine
    if T <: Complex
        return 0.0# currently complex case has different mappings, though both valid
    else
        @test xy ≈ xy_aff
    end
    return 0.0
end

## First, do common tests
stuff = [
        Dict(
            "constr" => IndGraph,
            "params" => (
                    (sprand(m, n, 0.2),),
                    (sprand(Complex{Float64}, m, n, 0.2),),
                    (rand(m, n),),
                    (rand(Complex{Float64}, m, n),),
                    (rand(sm, n),),
                    (rand(Complex{Float64}, sm, n),),
            ),
            "args"   => (
                    randn(m + n),
                    randn(m + n)+im * randn(m + n),
                    randn(m + n),
                    randn(m + n)+im * randn(m + n),
                    randn(sm + n),
                    randn(sm + n)+im * randn(sm + n)
            )
      ),
]

for i = 1:length(stuff)
  constr = stuff[i]["constr"]

  if haskey(stuff[i], "wrong")
    for j = 1:length(stuff[i]["wrong"])
      wrong = stuff[i]["wrong"][j]
      @test_throws ErrorException constr(wrong...)
    end
  end

  for j = 1:length(stuff[i]["params"])
    params = stuff[i]["params"][j]
    x      = stuff[i]["args"][j]
    f = constr(params...)

    predicates_test(f)

##### argument split
    c = view(x, 1:f.n)
    d = view(x, f.n + 1:f.n + f.m)
    ax = zero(c)
    ay = zero(d)

##### just call f
    fx = call_test(f, x)

##### compute prox with default gamma
    y, fy = prox_test(f, x)

##### compute prox with random gamma
    gam = 5*rand()
    y, fy = prox_test(f, x, gam)

##### test calls to prox! with more signatures
    prox!(ax, ay, f, c, d)
    @test f(ax, ay) <= TOL_ASSERT
    ax_naive, ay_naive, fv_naive = ProximalOperators.prox_naive(f, c, d, 1.0)
    @test f(ax_naive, ay_naive) <= TOL_ASSERT

    prox!((ax, ay), f, (c, d))
    @test f((ax, ay)) <= TOL_ASSERT
    axy_naive, fv_naive = ProximalOperators.prox_naive(f, (c, d))
    @test f(axy_naive) <= TOL_ASSERT

##### test against IndAffine
    test_against_IndAffine(f, params[1], x)
  end
end


text:
import Base: isprime

## http://docs.sympy.org/dev/modules/ntheory.html?highlight=ntheory#sympy.ntheory.factor_
ntheory_sympy_methods = (
                         :sieve,
                         :prime, :primepi,
                         :isprime, # import
                         :nextprime, :prevprime, :primorial,
                         :trailing,
                         :multiplicity, :perfect_power,
                         :pollard_rho, :pollard_pm1,
                         :factorint, :primefactors, :divisors, :divisor_count,
                         :totient,
                         :binomial_coefficients, :binomial_coefficients_list,
                         :multinomial_coefficients,
                         :npartitions,
                         :n_order, :is_primitive_root, :is_quad_residue,
                         :legendre_symbol, :jacobi_symbol
                         )

"""

[primerange](http://docs.sympy.org/dev/modules/ntheory.html): Generate all prime numbers in the range [a, b).

This is a replication of `primes` to which we dispatch.

"""
function primerange(a::Sym,b)
    a = max(2,convert(Int, a))
    b = max(a+1, convert(Int, b))
    primes(a,b)
end
export primerange


text:
export build

"""
    parse_problem(terms::Tuple, solver::ForwardBackwardSolver)

Takes as input a tuple containing the terms defining the problem and the solver.

Returns a tuple containing the optimization variables and the problem terms
to be fed into the solver.

# Example

```julia
julia> x = Variable(4)
Variable(Float64, (4,))

julia> A, b = randn(10,4), randn(10);

julia> p = problem( ls(A*x - b ) , norm(x) <= 1 );

julia> StructuredOptimization.parse_problem(p, ForwardBackward());
```
"""
function parse_problem(terms::Tuple, solver::T) where T <: ForwardBackwardSolver
  x = extract_variables(terms)
  # Separate smooth and nonsmooth
  smooth, nonsmooth = split_smooth(terms)
  if is_proximable(nonsmooth)
    g = extract_proximable(x, nonsmooth)
    kwargs = Dict{Symbol, Any}(:g => g)
    if !isempty(smooth)
      if is_linear(smooth)
        f = extract_functions(smooth)
        A = extract_operators(x, smooth)
        kwargs[:A] = A
      else  # ??
        f = extract_functions_nodisp(smooth)
        A = extract_affines(x, smooth)
        f = PrecomposeNonlinear(f, A)
      end
      kwargs[:f] = f
    end
    return (x, kwargs)
  end
  error("Sorry, I cannot parse this problem for solver of type $(T)")
end

export solve

"""
    solve(terms::Tuple, solver::ForwardBackwardSolver)

Takes as input a tuple containing the terms defining the problem and the solver options.

Solves the problem returning a tuple containing the iterations taken and the build solver.

# Example

```julia
julia> x = Variable(4)
Variable(Float64, (4,))

julia> A, b = randn(10,4), randn(10);

julia> p = problem(ls(A*x - b ), norm(x) <= 1);

julia> solve(p, ForwardBackward());

julia> ~x
```
"""
function solve(terms::Tuple, solver::ForwardBackwardSolver)
    x, kwargs = parse_problem(terms, solver)
    x_star, it = solver(~x; kwargs...)
    ~x .= x_star
    return x, it
end


text:
let # empty
  matrix = CoverMatrix(Dict{Any, Set{Any}}())

  @test size(matrix) == (0, 0)
  @test exact_cover_exists(matrix)
  @test exact_cover_count(matrix) == 1
  @test isempty(find_exact_cover(matrix))
end

let # 1-element
  set = CoverMatrix(Dict(1 => Set(23)))

  @test size(set) == (1, 1)
  @test exact_cover_exists(set)
  @test exact_cover_count(set) == 1
  @test find_exact_cover(set) == Set([1])
end

let # twice the same 1-element subset
  set = CoverMatrix(Dict(
    1 => Set(23),
    2 => Set(23)
  ))

  @test size(set) == (2, 1)
  @test exact_cover_exists(set)
  @test exact_cover_count(set) == 2

  Set(exact_cover_producer(set)) == Set([Set(1), Set(2)])
end

let # single subset, 2 elements
  set = CoverMatrix(Dict(1 => Set([23, 42])))

  @test size(set) == (1, 2)
  @test exact_cover_exists(set)
  @test exact_cover_count(set) == 1
  @test find_exact_cover(set) == Set(1)
end

let # 2 subsets of 1 element each
  set = CoverMatrix(Dict(
    1 => Set(23),
    2 => Set(42)
  ))

  @test size(set) == (2, 2)
  @test exact_cover_exists(set)
  @test exact_cover_count(set) == 1
  @test find_exact_cover(set) == Set([1, 2])
end

let # trivial decision
  set = CoverMatrix(Dict(
    1 => Set(23),
    2 => Set([23, 42])
  ))

  @test size(set) == (2, 2)
  @test exact_cover_exists(set)
  @test exact_cover_count(set) == 1
  @test find_exact_cover(set) == Set(2)
end

let # impossible decision
  set = CoverMatrix(Dict(
    1 => Set([12, 23]),
    2 => Set([23, 42])
  ))

  @test size(set) == (2, 3)
  @test !exact_cover_exists(set)
  @test exact_cover_count(set) == 0
  @test find_exact_cover(set) == nothing
  @test set |> exact_cover_producer |> collect |> length == 0
end

let # more complex key types
  set = CoverMatrix(Dict(
    1 => Set(23),
    :symbol => Set(42),
    (:tagged, 1, 2) => Set(123),
    [1, 2, 3] => Set(321)
  ))

  @test size(set) == (4, 4)
  @test exact_cover_exists(set)
  @test exact_cover_count(set) == 1
  @test Set(find_exact_cover(set)) == Set(Any[1, :symbol, (:tagged, 1, 2), [1, 2, 3]])
end

let # more complex value types
  set = CoverMatrix(Dict(
    1 => Set(Any[(:tagged, 1)]),
    2 => Set(Any[(:tagged, 2), (:tagged, 2)]),
    3 => Set(Any[(:tagged, 3), (:tagged, 4)]),
    4 => Set(Any[5, :symbol, (:tagged, 6)]),
    5 => Set(Any[:symbol, (:tagged, 6)])
  ))

  @test size(set) == (5, 7)
  @test exact_cover_count(set) == 1
  @test Set(find_exact_cover(set)) == Set([1, 2, 3, 4])
end

let # user-defined value types
  immutable Point
    x::Integer
    y::Integer
  end

  set = CoverMatrix(Dict(
    1 => Set(Any[Point(1, 2), Point(2, 3)]),
    2 => Set(Any[Point(2, 3)]),
    3 => Set(Any[(1, 2)])
  ))

  @test size(set) == (3, 3)
  @test exact_cover_count(set) == 1
  @test Set(find_exact_cover(set)) == Set([1, 3])
end


text:
#set up a matrix that has the form we're looking for.
M = [1.0 2.0 3.0
     2.0 5.0 5.0
     2.0 2.0 6.0]
new_row = zeros(3)
SigmoidNumbers.get_unscaled_replacement_row!(new_row, M, 2, zeros(3), Quire(Float64))
@test new_row == [0, 1, -1]

#do it with a bigger matrix.
M = [1.0 2.0 3.0
     0.0 1.0 -1.0
     2.0 2.0 6.0]
new_row = zeros(3)
SigmoidNumbers.get_unscaled_replacement_row!(new_row, M, 3, zeros(3), Quire(Float64))
@test new_row == [0, 0, -2]

M = [1.0 2.0 3.0
     2.0 5.0 5.0
     2.0 2.0 6.0]
v = ones(3)

SigmoidNumbers.row_echelon!(M, v, Quire(Float64))
@test M == [1.0 2.0 3.0
            0.0 1.0 -1.0
            0.0 0.0 1.0]
@test v == [1.0, -1.0, 1.5]

r = ones(3)
SigmoidNumbers.solve_row_echelon!(r, M, v, Quire(Float64))
@test r == [-4.5, 0.5, 1.5]

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
# # Hydraulic functions

# This tutorial shows how to specify the hydraulic functions
# used in Richard's equation. In particular,
# we show how to choose the formalism for matric potential and hydraulic
# conductivity, and how to make the hydraulic conductivity account for
# the presence of ice as well as the temperature dependence of the
# viscosity of liquid water.

# # Preliminary setup

# External modules
using Plots

# ClimateMachine modules
using ClimateMachine
using ClimateMachine.Land
using ClimateMachine.Land.SoilWaterParameterizations

FT = Float32;
# # Specifying a hydraulics model
# ClimateMachine's Land model allows the user to pick between two hydraulics models,
# that of van Genuchten [vanGenuchten1980](@cite) or that of Brooks and Corey, see [BrooksCorey1964](@cite) or [Corey1977](@cite). The
# same model is consistently used for the matric potential
# and hydraulic conductivity.

# The van Genuchten model requires two free parameters, `α` and `n`.
# A third parameter, `m`, is computed from `n`. Of these, only `α` carries
# units, of inverse meters. The Brooks and Corey model also uses
# two free parameters, `ψ_b`, the magnitude of the matric potential at saturation,
#  and a constant `M`. `ψ_b` carries units of meters. These parameter sets are stored in
# either the [`vanGenuchten`](@ref ClimateMachine.Land.SoilWaterParameterizations.vanGenuchten) or the
# [`BrooksCorey`](@ref ClimateMachine.Land.SoilWaterParameterizations.BrooksCorey)
# hydraulics model structures (more details below). These parameters are enough to compute the matric potential.

# The hydraulic conductivity
# requires an additional parameter, `Ksat` (m/s), which is the hydraulic conductivity
# in saturated soil. This parameter is
# not stored in the hydraulics model, but rather as part of the
# [`WaterParamFunctions`](@ref ClimateMachine.Land.WaterParamFunctions), which stores
# other parameters needed for the soil water modeling.

# Below we show how to create two concrete examples of these hydraulics models,
# for sandy loam ([Bonan19a](@cite)). Note that the parameters chosen are a function of soil type,
# and that the parameters are converted to type `FT` internally.
vg_α = 7.5 # m^-1
vg_n = 1.89
hydraulics = vanGenuchten(FT; α = vg_α, n = vg_n);

ψ_sat = 0.218 # m
Mval = 0.2041
hydraulics_bc = BrooksCorey(FT; ψb = ψ_sat, m = Mval);
# # Matric Potential
# The matric potential `ψ` reflects the negative pressure of water
# in unsaturated soil. The negative pressure (suction) of water arises
# because of adhesive forces between water and soil.

# The van Genuchten expression for matric potential is
# ``
# ψ = -\frac{1}{α} S_l^{-1/(nm)}\times (1-S_l^{1/m})^{1/n},
# ``

# and the Brooks and Corey expression is
# ``
# ψ = -ψ_b S_l^{-M}.
# ``

# Here `S_l` is the effective saturation of liquid water, `θ_l/ν`, where `ν` is
# porosity of the soil. We generally neglect the residual pore space in the CliMA model,
# but the user can set the parameter in the
# [`WaterParamFunctions`](@ref ClimateMachine.Land.WaterParamFunctions) structure if it is
# desired.

# In the CliMA code, we use [multiple dispatch](https://en.wikipedia.org/wiki/Multiple_dispatch).
# With multiple dispatch, a function can have many
# ways of executing (called methods), depending on the *type* of the
# variables passed in. A simple example of multiple dispatch is the division operation.
# Integer division takes two numbers as input, and returns an integer - ignoring the decimal.
# Float division takes two numbers as input, and returns a floating point number, including the decimal.
# In Julia, we might write these as:

# ```julia
# function division(a::Int, b::Int)
#      return floor(Int, a/b)
# end
# ```
# ```julia
# function division(a::Float64, b::Float64)
#      return a/b
# end
# ```

# We can see that `division` is now a function with two methods.

# ```julia
# julia> division
# division (generic function with 2 methods)
# ```

# Now, using the same function signature, we can carry out integer
# division or floating point division, depending on the types of the
# arguments:

# ```julia
# julia> division(1,2)
# 0
#
# julia> division(1.0,2.0)
# 0.5
# ```

# Here are more pertinent examples:
# Based on our choice of `FT = Float32`,

# ```julia
# julia> typeof(hydraulics)
# vanGenuchten{Float32,Float32,Float32,Float32}
# ```

# but meanwhile,

# ```julia
# julia> typeof(hydraulics_bc)
# BrooksCorey{Float32,Float32,Float32}
# ```

# The function `matric_potential` will execute different methods
# depending on if we pass a hydraulics model of type `vanGenuchten` or
# `BrooksCorey`. In both cases, it will return the correct value
# for `ψ`.

# Let's plot the matric potential as a function of the effective saturation `S_l = θ_l/ν`,
# which can range from zero to one.
S_l = FT.(0.01:0.01:0.99)
ψ = matric_potential.(Ref(hydraulics), S_l)
ψ_bc = matric_potential.(Ref(hydraulics_bc), S_l)
plot(
    S_l,
    log10.(-ψ),
    xlabel = "effective saturation",
    ylabel = "Log10(|ψ|)",
    label = "van Genuchten",
)
plot!(S_l, log10.(-ψ_bc), label = "Brooks and Corey")
savefig("bc_vg_matric_potential.png")
# ![](bc_vg_matric_potential.png)

# The steep slope in
# `ψ` near saturated and completely dry soil are part of the reason
# why Richard's equation is such a challenging numerical problem.

# # Hydraulic conductivity
# The hydraulic conductivity is a more complex function than the matric potential,
# as it depends on the temperature of the water, the volumetric ice fraction, and
# the volumetric liquid water fraction. It also depends on the hydraulics model
# chosen.

# We represent the hydraulic conductivity `K` as the product of four factors:
# `Ksat`, an impedance factor (which accounts for the effect of ice on conductivity)
# a viscosity factor (which accounts for the effect of temperature on the
# viscosity of liquid water, and how that in turn affects conductivity)
# and a moisture factor (which accounts for the effect of liquid water, and is determined by the hydraulics model).
# We are going to calculate `K = Ksat × viscosity factor × impedance factor × moisture factor`.
# In the code, each of these factors is
# computed by a function with multiple methods, except for `Ksat`.
# Like we defined new type
# classes for `vanGenuchten` and `BrooksCorey`, we also created new type classes
# for the impedance choice, the viscosity choice, and the moisture choice.

# The function [`viscosity_factor`](@ref ClimateMachine.Land.SoilWaterParameterizations.viscosity_factor)
# takes as arguments the temperature of the soil and the
# viscosity model desired, and returns the factor `k_v` by which the hydraulic conductivity is scaled.
# One option is to account for this effect:

# ``
# k_v = e^{γ (T-T_{\rm ref})}
# ``

# where γ = 0.0264/K and ``T_{\rm ref}`` = 288K.

# For example, at the freezing point of water, using the default values
# for γ and T_ref, viscosity reduces the conductivity by a third:
viscous_effect_model = TemperatureDependentViscosity{FT}();
viscosity_factor(viscous_effect_model, FT(273.15))

# The other option is to ignore this effect:

# ``
# k_v = 1
# ``

# This is the default approach.
no_viscous_effect_model = ConstantViscosity{FT}();
viscosity_factor(no_viscous_effect_model, FT(273.15))

# Very similarly, the function
# [`impedance_factor`](@ref ClimateMachine.Land.SoilWaterParameterizations.impedance_factor)
# takes as arguments the liquid water and ice
# volumetric fractions in the soil, as well as the impedance model being used, and returns
# the factor `k_i` by which the hydraulic conductivity is scaled.
# One option is to account for this effect:

# ``
# k_i = 10^{-Ω f_i},
# ``

# where `Ω = 7` is an empirical factor and
# `f_i` is the ratio of the volumetric
# ice fraction to total volumetric water fraction  ([Lundin1990](@cite)).

# For example, with ``\theta_i = \theta_l``, or f_i = 0.5, ice reduces the conductivity by over 1000x.
impedance_effect_model = IceImpedance{FT}();
impedance_factor(impedance_effect_model, FT(0.5))

# The other option is to ignore this effect:

# ``
# k_i = 1
# ``

# This is the default approach.
no_impedance_effect_model = NoImpedance{FT}();
impedance_factor(no_impedance_effect_model, FT(0.5))

# As for the moisture dependence of hydraulic conductivity, it can also be either
# independent of moisture, or dependent on moisture. If it is dependent on moisture,
# the specific function evaluated is dictated by the hydraulics model.
# The [`moisture_factor`](@ref ClimateMachine.Land.SoilWaterParameterizations.moisture_factor)
# for the van Genuchten model is (denoting it as ``k_m``)

# ``
#  k_m = \sqrt{S_l}[1-(1-S_l^{1/m})^m]^2,
# ``

# for ``S_l < 1``,

# and for the Brooks and Corey model it is

# ``
# k_m = S_l^{2M+3},
# ``

# also for ``S_l<1``. When ``S_l\geq 1``, ``k_m = 1`` for each model.

# Let's put all these factors together now. Below
# we choose additional parameters, consistent with the hydraulics parameters
# for sandy loam ([Bonan19a](@cite)), and show how hydraulic conductivity varies with
# liquid water content, in the case without ice impedance or temperature effects.
Ksat = FT(4.42 / (3600 * 100))
T = FT(0.0)
f_i = FT(0.0)
K =
    hydraulic_conductivity.(
        Ref(Ksat),
        Ref(impedance_factor(NoImpedance{FT}(), f_i)),
        Ref(viscosity_factor(ConstantViscosity{FT}(), T)),
        moisture_factor.(Ref(MoistureDependent{FT}()), Ref(hydraulics), S_l),
    );

# Let's also compute `K` when we include the effects of temperature
# and ice on the hydraulic conductivity.
# In the cases where a
# [`TemperatureDependentViscosity`](@ref ClimateMachine.Land.SoilWaterParameterizations.TemperatureDependentViscosity)
# or
# [`IceImpedance`](@ref ClimateMachine.Land.SoilWaterParameterizations.IceImpedance)
# type is passed, the correct factors are calculated,
# based on the temperature `T` and volumetric ice fraction `θ_i`.

T = FT(273.15)
S_i = FT(0.1); # = θ_i/ν
# The total volumetric water fraction cannot
# exceed unity, so the effective liquid water saturation
# should have a max of 1-S_i.
S_l_accounting_for_ice = FT.(0.01:0.01:(0.99 - S_i))
f_i = S_i ./ (S_l_accounting_for_ice .+ S_i)
K_w_factors =
    hydraulic_conductivity.(
        Ref(Ksat),
        impedance_factor.(Ref(NoImpedance{FT}()), f_i),
        Ref(viscosity_factor(ConstantViscosity{FT}(), T)),
        moisture_factor.(
            Ref(MoistureDependent{FT}()),
            Ref(hydraulics),
            S_l_accounting_for_ice,
        ),
    );
plot(
    S_l,
    log10.(K),
    xlabel = "total effective saturation, (θ_i+θ_l)/ν",
    ylabel = "Log10(K)",
    label = "Base case",
    legend = :bottomright,
)
plot!(
    S_l_accounting_for_ice .+ S_i,
    log10.(K_w_factors),
    label = "θ_i = 0.1, T = 273.15",
)
savefig("T_ice_K.png")
# ![](T_ice_K.png)

# If the user is not considering phase transitions
# and does not add in Freeze/Thaw source terms, the default is for zero
# ice in the model, for all time and space. In this case the ice impedance
# factor evaluates to 1 regardless of which type is passed.

# # Other features
# The user also has the choice of making the conductivity constant by choosing
# [`MoistureIndependent`](@ref ClimateMachine.Land.SoilWaterParameterizations.MoistureIndependent)
# along with
# [`ConstantViscosity`](@ref ClimateMachine.Land.SoilWaterParameterizations.ConstantViscosity)
# and
# [`NoImpedance`](@ref ClimateMachine.Land.SoilWaterParameterizations.NoImpedance).
# This is useful for debugging!
no_moisture_dependence = MoistureIndependent{FT}()
K_constant =
    hydraulic_conductivity.(
        Ref(Ksat),
        Ref(FT(1.0)),
        Ref(FT(1.0)),
        moisture_factor.(Ref(no_moisture_dependence), Ref(hydraulics), S_l),
    );
# ```julia
# julia> unique(K_constant)
# 1-element Array{Float32,1}:
#  1.2277777f-5
# ```

# Note that choosing this option does not mean the matric potential
# is constant, as a hydraulics model is still required and employed.

# And, lastly, you might also find it helpful in debugging
# to be able to turn off the flow of water by setting `Ksat = 0`.

# # References
# - [vanGenuchten1980](@cite)
# - [BrooksCorey1964](@cite)
# - [Corey1977](@cite)
# - [Lundin1990](@cite)
# - [Bonan19a](@cite)


text:
module TestInterface

using Functors
using Random
using Test

function test_interface(rng::AbstractRNG, lik, out_dist, D_in=1; functor_args=())
    N = 10
    T = Float64 # TODO test Float32 as well
    f, fs = if D_in == 1
        f = randn(rng, T)
        fs = randn(rng, T, 10)
        f, fs
    else
        f = randn(rng, T, D_in)
        fs = [randn(rng, T, D_in) for _ in 1:N]
        f, fs
    end
    # Check if likelihood produces the correct distribution
    # and is sampleable
    @test lik(f) isa out_dist
    @test_nowarn rand(rng, lik(f))

    # Check if the likelihood samples are of correct length
    @test length(rand(rng, lik(fs))) == N

    # Check if functor works properly
    xs, re = Functors.functor(lik)
    @test lik == re(xs)
    if isempty(functor_args)
        @test xs === ()
    else
        @test keys(xs) == functor_args
    end

    return nothing
end

@doc raw"""
    test_interface([rng::AbstractRNG,] lik, out_dist; functor_args=())

This function provides unified method to check the interface of the various likelihoods
defined. It checks if the likelihood produces the correct distribution, length of likelihood
samples is correct and if the functor works as intended.

...
## Arguments
- `lik`: the likelihood to test the interface of
- `out_dist::Type{<:Distribution}`: the type of distribution the likelihood should return
- `D_in::Int=1` : The input dimension of the likelihood

## Keyword arguments
- `functor_args=()`: a collection of symbols of arguments to match functor parameters with.
...
"""
function test_interface(lik, out_dist, D_in=1; kwargs...)
    return test_interface(Random.GLOBAL_RNG, lik, out_dist, D_in; kwargs...)
end

end