text:
const BigFloats = Union{BigFloat,Complex{BigFloat}}

if VERSION < v"0.7-"
    import Base.FFTW: fft, fft!, rfft, irfft, ifft, conv, dct, idct, dct!, idct!,
                        plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
                        plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct
else
    import FFTW: fft, fft!, rfft, irfft, ifft, dct, idct, dct!, idct!,
                        plan_fft!, plan_ifft!, plan_dct!, plan_idct!,
                        plan_fft, plan_ifft, plan_rfft, plan_irfft, plan_dct, plan_idct
    import DSP: conv
end

# The following implements Bluestein's algorithm, following http://www.dsprelated.com/dspbooks/mdft/Bluestein_s_FFT_Algorithm.html
# To add more types, add them in the union of the function's signature.
function fft(x::Vector{T}) where T<:BigFloats
    n = length(x)
    ispow2(n) && return fft_pow2(x)
    ks = range(zero(real(T)),stop=n-one(real(T)),length=n)
    Wks = exp.((-im).*convert(T,π).*ks.^2 ./ n)
    xq, wq = x.*Wks, conj([exp(-im*convert(T,π)*n);reverse(Wks);Wks[2:end]])
    return Wks.*conv(xq,wq)[n+1:2n]
end

function fft!(x::Vector{T}) where T<:BigFloats
    x[:] = fft(x)
    return x
end

# add rfft for BigFloat, by calling fft
#  this creates ToeplitzMatrices.rfft, so avoids changing rfft

rfft(v::Vector{T}) where T<:BigFloats = fft(v)[1:div(length(v),2)+1]
function irfft(v::Vector{T},n::Integer) where T<:BigFloats
    @assert n==2length(v)-1
    r = Vector{Complex{BigFloat}}(undef, n)
    r[1:length(v)]=v
    r[length(v)+1:end]=reverse(conj(v[2:end]))
    real(ifft(r))
end

ifft(x::Vector{T}) where {T<:BigFloats} = conj!(fft(conj(x)))/length(x)
function ifft!(x::Vector{T}) where T<:BigFloats
    x[:] = ifft(x)
    return x
end

function conv(u::StridedVector{T}, v::StridedVector{T}) where T<:BigFloats
    nu,nv = length(u),length(v)
    n = nu + nv - 1
    np2 = nextpow(2,n)
    append!(u,zeros(T,np2-nu)),append!(v,zeros(T,np2-nv))
    y = ifft_pow2(fft_pow2(u).*fft_pow2(v))
    #TODO This would not handle Dual/ComplexDual numbers correctly
    y = T<:Real ? real(y[1:n]) : y[1:n]
end

# This is a Cooley-Tukey FFT algorithm inspired by many widely available algorithms including:
# c_radix2.c in the GNU Scientific Library and four1 in the Numerical Recipes in C.
# However, the trigonometric recurrence is improved for greater efficiency.
# The algorithm starts with bit-reversal, then divides and conquers in-place.
function fft_pow2!(x::Vector{T}) where T<:BigFloat
    n,big2=length(x),2one(T)
    nn,j=n÷2,1
    for i=1:2:n-1
        if j>i
            x[j], x[i] = x[i], x[j]
            x[j+1], x[i+1] = x[i+1], x[j+1]
        end
        m = nn
        while m ≥ 2 && j > m
            j -= m
            m = m÷2
        end
        j += m
    end
    logn = 2
    while logn < n
        θ=-big2/logn
        wtemp = sinpi(θ/2)
        wpr, wpi = -2wtemp^2, sinpi(θ)
        wr, wi = one(T), zero(T)
        for m=1:2:logn-1
            for i=m:2logn:n
                j=i+logn
                mixr, mixi = wr*x[j]-wi*x[j+1], wr*x[j+1]+wi*x[j]
                x[j], x[j+1] = x[i]-mixr, x[i+1]-mixi
                x[i], x[i+1] = x[i]+mixr, x[i+1]+mixi
            end
            wr = (wtemp=wr)*wpr-wi*wpi+wr
            wi = wi*wpr+wtemp*wpi+wi
        end
        logn = logn << 1
    end
    return x
end

function fft_pow2(x::Vector{Complex{T}}) where T<:BigFloat
    y = interlace(real(x),imag(x))
    fft_pow2!(y)
    return complex.(y[1:2:end],y[2:2:end])
end
fft_pow2(x::Vector{T}) where {T<:BigFloat} = fft_pow2(complex(x))

function ifft_pow2(x::Vector{Complex{T}}) where T<:BigFloat
    y = interlace(real(x),-imag(x))
    fft_pow2!(y)
    return complex.(y[1:2:end],-y[2:2:end])/length(x)
end

function dct(a::AbstractVector{Complex{BigFloat}})
	N = big(length(a))
    c = fft([a; flipdim(a,1)])
    d = c[1:N]
    d .*= exp.((-im*big(pi)).*(0:N-1)./(2*N))
    d[1] = d[1] / sqrt(big(2))
    lmul!(inv(sqrt(2N)), d)
end

dct(a::AbstractArray{BigFloat}) = real(dct(complex(a)))

function idct(a::AbstractVector{Complex{BigFloat}})
	N = big(length(a))
    b = a * sqrt(2*N)
    b[1] = b[1] * sqrt(big(2))
    shift = exp.(-im * 2 * big(pi) * (N - big(1)/2) * (0:(2N-1)) / (2 * N))
    b = [b; 0; -flipdim(b[2:end],1)] .* shift
    c = ifft(b)
    flipdim(c[1:N],1)
end

idct(a::AbstractArray{BigFloat}) = real(idct(complex(a)))

dct!(a::AbstractArray{T}) where {T<:BigFloats} = (b = dct(a); a[:] = b)
idct!(a::AbstractArray{T}) where {T<:BigFloats} = (b = idct(a); a[:] = b)

# dummy plans
struct DummyFFTPlan{T,inplace} <: Plan{T} end
struct DummyiFFTPlan{T,inplace} <: Plan{T} end
struct DummyrFFTPlan{T,inplace} <: Plan{T} end
struct DummyirFFTPlan{T,inplace} <: Plan{T} end
struct DummyDCTPlan{T,inplace} <: Plan{T} end
struct DummyiDCTPlan{T,inplace} <: Plan{T} end

for (Plan,iPlan) in ((:DummyFFTPlan,:DummyiFFTPlan),
                     (:DummyrFFTPlan,:DummyirFFTPlan),
                     (:DummyDCTPlan,:DummyiDCTPlan))
   @eval begin
       Base.inv(::$Plan{T,inplace}) where {T,inplace} = $iPlan{T,inplace}()
       Base.inv(::$iPlan{T,inplace}) where {T,inplace} = $Plan{T,inplace}()
    end
end

for (Plan,ff,ff!) in ((:DummyFFTPlan,:fft,:fft!),
                      (:DummyiFFTPlan,:ifft,:ifft!),
                      (:DummyrFFTPlan,:rfft,:rfft!),
                      (:DummyirFFTPlan,:irfft,:irfft!),
                      (:DummyDCTPlan,:dct,:dct!),
                      (:DummyiDCTPlan,:idct,:idct!))
    @eval begin
        *(p::$Plan{T,true}, x::StridedArray{T,N}) where {T,N} = $ff!(x)
        *(p::$Plan{T,false}, x::StridedArray{T,N}) where {T,N} = $ff(x)
        function LAmul!(C::StridedVector, p::$Plan, x::StridedVector)
            C[:] = $ff(x)
            C
        end
    end
end

plan_fft!(x::Vector{T}) where {T<:BigFloats} = DummyFFTPlan{Complex{BigFloat},true}()
plan_ifft!(x::Vector{T}) where {T<:BigFloats} = DummyiFFTPlan{Complex{BigFloat},true}()

#plan_rfft!{T<:BigFloats}(x::Vector{T}) = DummyrFFTPlan{Complex{BigFloat},true}()
#plan_irfft!{T<:BigFloats}(x::Vector{T},n::Integer) = DummyirFFTPlan{Complex{BigFloat},true}()
plan_dct!(x::Vector{T}) where {T<:BigFloats} = DummyDCTPlan{T,true}()
plan_idct!(x::Vector{T}) where {T<:BigFloats} = DummyiDCTPlan{T,true}()

plan_fft(x::Vector{T}) where {T<:BigFloats} = DummyFFTPlan{Complex{BigFloat},false}()
plan_ifft(x::Vector{T}) where {T<:BigFloats} = DummyiFFTPlan{Complex{BigFloat},false}()
plan_rfft(x::Vector{T}) where {T<:BigFloats} = DummyrFFTPlan{Complex{BigFloat},false}()
plan_irfft(x::Vector{T},n::Integer) where {T<:BigFloats} = DummyirFFTPlan{Complex{BigFloat},false}()
plan_dct(x::Vector{T}) where {T<:BigFloats} = DummyDCTPlan{T,false}()
plan_idct(x::Vector{T}) where {T<:BigFloats} = DummyiDCTPlan{T,false}()

function interlace(a::Vector{S},b::Vector{V}) where {S<:Number,V<:Number}
    na=length(a);nb=length(b)
    T=promote_type(S,V)
    if nb≥na
        ret=zeros(T,2nb)
        ret[1:2:1+2*(na-1)]=a
        ret[2:2:end]=b
        ret
    else
        ret=zeros(T,2na-1)
        ret[1:2:end]=a
        if !isempty(b)
            ret[2:2:2+2*(nb-1)]=b
        end
        ret
    end
end


text:
################################################################################
#
#  Some types
#
################################################################################

# Given maximal orders O_1, ..., O_k in the algebra A over a field K, this
# helps to construct field extensions of K, whose maximal order is a submodule
# of one of the O_i.
mutable struct FieldOracle{S, T, U, M}
  algebra::S
  maximal_orders::Vector{T}

  hnf_basis_pmats::Vector{<: PMat}
  small_elements::Vector{U}
  fields::Vector{Vector{M}}
  phase::Int
  rounds::Int
  last_generated_field::Vector{Int}
  last_returned_field::Vector{Int}
  rand_coeff_bound::Int

  function FieldOracle{S, T, U, M}(A::S, orders::Vector{T}) where { S, T, U, M }
    z = new{S, T, U, M}()
    z.algebra = A
    z.maximal_orders = orders

    if A isa AbsAlgAss{nf_elem}
      # basis_pmatrices of orders are not required to be in HNF, so we compute them here
      z.hnf_basis_pmats = Vector{typeof(basis_pmatrix(orders[1], copy = false))}(undef, length(orders))
      for i = 1:length(orders)
        z.hnf_basis_pmats[i] = pseudo_hnf(basis_pmatrix(orders[i], copy = false), :lowerleft, true)
      end
    end

    z.small_elements = Vector{U}()
    for O in orders
      append!(z.small_elements, small_elements(O))
    end

    z.fields = Vector{Vector{M}}(undef, length(orders))

    for i = 1:length(orders)
      z.fields[i] = Vector{M}()
    end

    z.phase = 1
    # Phase 1: Construct fields using all possible sums of elements in
    #          small_elements with <= z.rounds summands.
    # Phase 2: Use random linear combinations of elements in small_elements.

    z.rounds = 10 # Magic number
    z.last_generated_field = zeros(Int, z.rounds)
    z.last_returned_field = zeros(Int, length(orders))
    z.rand_coeff_bound = Int(ceil(log(length(z.small_elements), 10^6))) # So we get around 10^6 random fields. Hopefully that's enough.

    return z
  end
end

function FieldOracle(A::AbsAlgAss{nf_elem}, orders::Vector{<: AlgAssRelOrd})
  return FieldOracle{typeof(A), typeof(orders[1]), elem_type(A), NfRelToAbsAlgAssMor}(A, orders)
end

function FieldOracle(A::AbsAlgAss{fmpq}, orders::Vector{<: AlgAssAbsOrd})
  return FieldOracle{typeof(A), typeof(orders[1]), elem_type(A), NfAbsToAbsAlgAssMor}(A, orders)
end

mutable struct NormCache{S, T, U, M, T2, T3}
  algebra::S
  maximal_orders::Vector{T}
  n::Int # n^2 == dim(algebra)

  field_oracle::FieldOracle{S, T, U, M}

  a::T2
  primes::Vector{T3}
  valuations::Vector{Int}
  fac_elem_mon::FacElemMon{S}
  partial_solutions::Vector{Dict{Set{Int}, Vector{FacElem{U, S}}}}
  solutions_mod_units::Vector{Dict{FacElem{U, S},  GrpAbFinGenElem}}

  # Only used if S <: AbsAlgAss{fmpq}
  norm_minus_one::Vector{U}

  # These fields are only set if S <: AbsAlgAss{nf_elem}
  UktoOk::MapUnitGrp
  GtoUk::Vector{GrpAbFinGenMap}
  GtoUk_surjective::BitVector
  fields_in_product::Vector{Vector{Tuple{NfRelToAbsAlgAssMor, NfToNfRel}}}

  function NormCache{S, T, U, M, T2, T3}(A::S, orders::Vector{T}, a::T2) where { S <: AbsAlgAss{nf_elem}, T, U, M, T2 <: NfAbsOrdElem, T3}
    primes = collect(keys(factor(a*parent(a))))
    vals = [ valuation(a, p) for p in primes ]
    z = NormCache{S, T, U, M, T2, T3}(A, orders, parent(a), primes, vals)
    z.a = a
    return z
  end

  function NormCache{S, T, U, M, T2, T3}(A::S, orders::Vector{T}, a::T2) where { S <: AbsAlgAss{fmpq}, T, U, M, T2 <: fmpz, T3}
    primes = collect(keys(factor(a).fac))
    vals = [ valuation(a, p) for p in primes ]
    z = NormCache{S, T, U, M, T2, T3}(A, orders, primes, vals)
    z.a = a
    return z
  end

  function NormCache{S, T, U, M, T2, T3}(A::S, orders::Vector{T}, Ok::NfAbsOrd, primes::Vector{T3}, valuations::Vector{Int}) where { S <: AbsAlgAss{nf_elem}, T, U, M, T2, T3 }
    z = new{S, T, U, M, T2, T3}()
    z.algebra = A
    z.maximal_orders = orders
    z.n = isqrt(dim(A))
    @assert z.n^2 == dim(A)

    z.field_oracle = FieldOracle(A, orders)

    Uk, UktoOk = unit_group(Ok)
    z.UktoOk = UktoOk
    z.primes = primes
    z.valuations = valuations
    z.fac_elem_mon = FacElemMon(A)
    z.partial_solutions = Vector{Dict{Set{Int}, Vector{FacElem{U, S}}}}(undef, length(orders))
    z.solutions_mod_units = Vector{Dict{FacElem{U, S}, GrpAbFinGenElem}}(undef, length(orders))
    z.GtoUk = Vector{GrpAbFinGenMap}(undef, length(orders))
    z.GtoUk_surjective = falses(length(orders))
    z.fields_in_product = Vector{Vector{Tuple{NfRelToAbsAlgAssMor, NfToNfRel}}}(undef, length(orders))
    for i = 1:length(orders)
      z.partial_solutions[i] = Dict{Set{Int}, Vector{FacElem{U, S}}}()
      z.partial_solutions[i][Set{Int}()] = [ z.fac_elem_mon() ]
      z.solutions_mod_units[i] = Dict{FacElem{U, S}, GrpAbFinGenElem}()
      z.GtoUk[i] = hom([ GrpAbFinGen(fmpz[])() ], [ Uk() ])
      z.fields_in_product[i] = Vector{Tuple{NfRelToAbsAlgAssMor, NfToNfRel}}()
    end

    return z
  end

  function NormCache{S, T, U, M, T2, T3}(A::S, orders::Vector{T}, primes::Vector{T3}, valuations::Vector{Int}) where { S <: AbsAlgAss{fmpq}, T, U, M, T2, T3 }
    z = new{S, T, U, M, T2, T3}()
    z.algebra = A
    z.maximal_orders = orders
    z.n = isqrt(dim(A))
    @assert z.n^2 == dim(A)

    z.field_oracle = FieldOracle(A, orders)
    z.norm_minus_one = Vector{U}(undef, length(orders))

    z.primes = primes
    z.valuations = valuations
    z.fac_elem_mon = FacElemMon(A)
    z.partial_solutions = Vector{Dict{Set{Int}, Vector{FacElem{U, S}}}}(undef, length(orders))
    z.solutions_mod_units = Vector{Dict{FacElem{U, S}, GrpAbFinGenElem}}(undef, length(orders))
    for i = 1:length(orders)
      z.partial_solutions[i] = Dict{Set{Int}, Vector{FacElem{U, S}}}()
      z.partial_solutions[i][Set{Int}()] = [ z.fac_elem_mon() ]
      z.solutions_mod_units[i] = Dict{FacElem{U, S}, GrpAbFinGenElem}()
    end

    return z
  end
end

function NormCache(A::AbsAlgAss{nf_elem}, orders::Vector{<: AlgAssRelOrd}, a::NfAbsOrdElem)
  return NormCache{typeof(A), typeof(orders[1]), elem_type(A), NfRelToAbsAlgAssMor, typeof(a), ideal_type(parent(a))}(A, orders, a)
end

function NormCache(A::AbsAlgAss{nf_elem}, orders::Vector{<: AlgAssRelOrd}, Ok::NfAbsOrd, primes::Vector{<: NfAbsOrdIdl}, valuations::Vector{Int})
  return NormCache{typeof(A), typeof(orders[1]), elem_type(A), NfRelToAbsAlgAssMor, elem_type(Ok), ideal_type(Ok)}(A, orders, Ok, primes, valuations)
end

function NormCache(A::AbsAlgAss{fmpq}, orders::Vector{<: AlgAssAbsOrd}, a::fmpz)
  return NormCache{typeof(A), typeof(orders[1]), elem_type(A), NfAbsToAbsAlgAssMor, fmpz, fmpz}(A, orders, a)
end

function NormCache(A::AbsAlgAss{fmpq}, orders::Vector{<: AlgAssAbsOrd}, primes::Vector{fmpz}, valuations::Vector{Int})
  return NormCache{typeof(A), typeof(orders[1]), elem_type(A), NfAbsToAbsAlgAssMor, fmpz, fmpz}(A, orders, primes, valuations)
end

################################################################################
#
#  Integral norm equations
#
################################################################################

# Returns the solution of the norm equation and the number of the order, in which
# the solution was found.
function norm_equation(orders::Vector{<: AlgAssAbsOrd}, a::fmpz)
  A = algebra(orders[1])
  if iszero(a)
    return A(), 1
  end
  d, i = norm_equation_fac_elem(orders, a)
  return evaluate(d), i
end

function norm_equation_fac_elem(orders::Vector{<: AlgAssAbsOrd}, a::fmpz)
  A = algebra(orders[1])
  @assert !iszero(a) # We cannot represent 0 as a FacElem
  if isone(a)
    return one(A), 1
  end
  NC = NormCache(A, orders, a)
  while true
    # We test the orders alternately until we find a solution in one of them
    for i = 1:length(orders)
      b, d = _norm_equation_absolute(NC, i)
      if b
        return d, i
      end
    end
  end
end

# This tests max_num_fields fields. If max_num_fields == 0, we try until we
# find a solution (although at some point NC.field_oracle does not return any
# new fields).
function _norm_equation_absolute(NC::NormCache, order_num::Int; max_num_fields::Int = 10)
  A = NC.algebra
  n = NC.n
  primes = NC.primes
  vals = NC.valuations
  full_set = Set(1:length(primes))
  partial_solutions = NC.partial_solutions[order_num]
  solutions_mod_units = NC.solutions_mod_units[order_num]

  m = 0
  while max_num_fields == 0 || m < max_num_fields
    m += 1
    KtoA = get_new_field(NC.field_oracle, order_num)
    K = domain(KtoA)
    OK = maximal_order(K)

    good_primes = __neq_find_good_primes(NC, OK)

    if !isempty(good_primes)
      if degree(K) != n
        n2 = degree(K)//n
        vals2 = [ Int(n2*vals[i]) for i = 1:length(vals) ]
      else
        vals2 = vals
      end

      remaining_primes = Set{Int}()
      # First search for solutions for single primes
      for p in good_primes
        sols = __neq_sunit(K, [ primes[p] ], [ vals2[p] ])
        if !isempty(sols)
          # We better evaluate s here: Some factors may not be integral which can
          # become a problem since the multiplication in A is non-commutative.
          # Also, multiplications and computing inverses in K are cheap.
          add_partial_solutions!(NC, order_num, Set(p), [ NC.fac_elem_mon(KtoA(evaluate(s))) for s in sols ])
        else
          push!(remaining_primes, p)
        end
      end

      if !isempty(remaining_primes)
        # Now the primes together, for which have not found a solution yet
        sols = __neq_sunit(K, [ primes[i] for i in remaining_primes ], [ vals2[i] for i in remaining_primes ])
        if !isempty(sols)
          add_partial_solutions!(NC, order_num, remaining_primes, [ NC.fac_elem_mon(KtoA(evaluate(s))) for s in sols ])
        elseif length(remaining_primes) < length(good_primes)
          # If this also failed, we test all primes together
          sols = __neq_sunit(K, [ primes[i] for i in good_primes ], [ vals2[i] for i in good_primes ])
          if !isempty(sols)
            add_partial_solutions!(NC, order_num, good_primes, [ NC.fac_elem_mon(KtoA(evaluate(s))) for s in sols ])
          end
        end
      end
    end

    if degree(K) == n
      UK, mUK = unit_group_fac_elem(OK)
      i = findfirst(x -> norm(mUK(x)) == -1, gens(UK))
      if i != nothing
        NC.norm_minus_one[order_num] = KtoA(evaluate(mUK(UK[i])))
      end
    end

    if !haskey(partial_solutions, full_set)
      continue
    end

    sols = partial_solutions[full_set]
    for s in sols
      if normred(s) != NC.a
        if !isassigned(NC.norm_minus_one, order_num)
          continue
        end
        s *= NC.norm_minus_one[order_num]
      end
      return true, s
    end
  end
  return false, NC.fac_elem_mon()
end

################################################################################
#
#  Integral norm equations - valuations only
#
################################################################################

# A special case needed for principal_generator_eichler.
# Finds a in O such that v_{p_i}(normred(a)) == v_i where p_i = primes[i] and
# v_i = valuations[i] and such that v_q(normred(a)) == 0 for all other primes q,
# assuming that such an element exists.
function _norm_equation_valuations_only(O::AlgAssAbsOrd, primes::Vector{fmpz}, valuations::Vector{Int})
  @assert !isempty(primes)
  @assert length(primes) == length(valuations)
  A = algebra(O)
  NC = NormCache(A, [ O ], primes, valuations)
  n = NC.n
  full_set = Set(1:length(primes))
  partial_solutions = NC.partial_solutions[1]

  while true
    KtoA = get_new_field(NC.field_oracle, 1, no_restriction = true)
    K = domain(KtoA)
    OK = maximal_order(K)

    good_primes = __neq_find_good_primes(NC, OK)

    if !isempty(good_primes)
      cache = Vector{Any}(undef, 2) # Used in __neq_find_sol_in_order
      if degree(K) != n
        n2 = degree(K)//n
        vals2 = [ Int(n2*valuations[i]) for i = 1:length(valuations) ]
      else
        vals2 = valuations
      end

      remaining_primes = Set{Int}()
      # First search for solutions for single primes
      for p in good_primes
        if haskey(partial_solutions, Set(p))
          continue
        end
        b, s = __neq_find_sol_in_order(O, KtoA, [ primes[p] ], [ vals2[p] ], cache)
        if b
          add_partial_solutions!(NC, 1, Set(p), [ NC.fac_elem_mon(s) ])
        else
          push!(remaining_primes, p)
        end
      end

      if !isempty(remaining_primes) && !haskey(partial_solutions, remaining_primes)
        # Now the primes together, for which have not found a solution yet
        b, s = __neq_find_sol_in_order(O, KtoA, [ primes[i] for i in remaining_primes ], [ vals2[i] for i in remaining_primes ], cache)
        if b
          add_partial_solutions!(NC, 1, remaining_primes, [ NC.fac_elem_mon(s) ])
        elseif length(remaining_primes) < length(good_primes)
          # If this also failed, we test all primes together
          b, s = __neq_find_sol_in_order(O, KtoA, [ primes[i] for i in good_primes ], [ vals2[i] for i in good_primes ], cache)
          if b
            add_partial_solutions!(NC, 1, good_primes, [ NC.fac_elem_mon(s) ])
          end
        end
      end
    end

    if haskey(partial_solutions, full_set)
      s = first(partial_solutions[full_set])
      return s
    end
  end
end

# Finds any a \in O \cap K such that v_p(nr(a)) = vals[i] for p = primes[i].
function __neq_find_sol_in_order(O::AlgAssAbsOrd, KtoA::NfAbsToAbsAlgAssMor, primes::Vector{fmpz}, vals::Vector{Int}, cache::Vector{Any})
  A = algebra(O)
  K = domain(KtoA)
  sols = __neq_sunit(K, primes, vals)
  if isempty(sols)
    return false, A()
  end

  K = domain(KtoA)
  sols_eval = Vector{elem_type(K)}()
  for s in sols
    s_eval = evaluate(s)
    t = KtoA(s_eval)
    if t in O
      return true, t
    end
    push!(sols_eval, s_eval)
  end

  # None of the found solutions lies in O.
  # We assume that equation_order(K) \subseteq O and try to adjust one of the
  # elements of sols by a unit.
  OK = maximal_order(K)
  UK, mUK = unit_group(OK)
  OE = equation_order(K)
  if !isassigned(cache, 1)
    cache[1] = OO_mod_F_mod_O_mod_F(OE)
  end
  G, GtoQ, OKtoQ = cache[1]
  sols2 = Vector{elem_type(K)}()
  for s in s_eval
    sinQ = OKtoQ(OK(s))
    if !isinvertible(sinQ)[1]
      push!(s, sols2)
    end
    # s is coprime to the conductor

    g = GtoQ\sinQ
    h = hom(UK, G, [ GtoQ\(OKtoQ(mUK(UK[i]))) for i = 1:ngens(UK) ])
    b, u = haspreimage(h, g)
    if b
      Q, toQ = quo(UK, kernel(h)[1])
      u = toQ\(toQ(u)) # Reduce the coefficient size (hopefully)
      s = s*inv(elem_in_nf(mUK(u), copy = false))
      return true, KtoA(s)
    end
  end

  if !isassigned(cache, 2)
    UE, mUE = unit_group(OE)
    UEinUK = [ mUK\(OK(elem_in_nf(mUE(UE[i]), copy = false))) for i = 1:ngens(UE) ]
    cache[2] = quo(UK, UEinUK)
  end
  Q, toQ = cache[2]
  for (i, g) in enumerate(Q)
    u = mUK(toQ\g)
    for s in sols2
      su = KtoA(s*u)
      if su in O
        return true, su
      end
    end
  end
  return false, A()
end

################################################################################
#
#  Some helpers
#
################################################################################

# Let S be the set of primes lying over primes.
# This finds a set of representatives of the S-units of K modulo units of O_K,
# whose norm has the valuations vals at the primes primes.
function __neq_sunit(K::AnticNumberField, primes::Vector{fmpz}, vals::Vector{Int})
  OK = maximal_order(K)
  primes_in_K = prime_ideals_over(OK, primes)
  class_group(K)
  SK, mSK = sunit_mod_units_group_fac_elem(primes_in_K)
  M = zero_matrix(FlintZZ, length(primes), ngens(SK))
  N = zero_matrix(FlintZZ, length(primes_in_K), ngens(SK))
  b = zero_matrix(FlintZZ, length(primes), 1)
  for c = 1:ngens(SK)
    g = mSK(SK[c])
    ng = norm(g)
    for r = 1:length(primes)
      M[r, c] = valuation(ng, primes[r])
      b[r, 1] = vals[r]
    end
    for r = 1:length(primes_in_K)
      N[r, c] = valuation(g, primes_in_K[r])
    end
  end

  sols = try
    solve_mixed(M, b, N)
  catch e
    if e isa UndefVarError
      error("This function requires Polymake.jl (original error message follows)")
    else
      rethrow(e)
    end
  end
  return elem_type(codomain(mSK))[ mSK(SK(s)) for s in sols ]
end

function __neq_find_good_primes(NC::NormCache, OK::NfAbsOrd)
  n = NC.n
  m = degree(nf(OK))
  mn = m//n
  good_primes = Set{Int}()
  for j = 1:length(NC.primes)
    v = NC.valuations[j]*mn
    if !isone(denominator(v))
      continue
    end
    p = NC.primes[j]
    pdec = prime_decomposition(OK, p)
    for (q, e) in pdec
      if q.splitting_type[2] <= v
        push!(good_primes, j)
        break
      end
    end
  end
  return good_primes
end

ext: jl


text:
# Orlando.jl -
module Orlando

using Flux
#using TightBinding # currently a bit too heavy + doesn't play well
using Nabla
using LinearAlgebra

println("    \"Life, it sings, or croons rather, like a kettle on a hob, Life, life what art thou?\" - Virginia Woolf, Orlando: a Biography")

# package code goes here

end # module


text:
"""
    NOCB <: Imputor

Fills in missing data using the Next Observation Carried Backward (NOCB) approach.
"""
struct NOCB <: Imputor end

"""
    impute!(imp::NOCB, ctx::Context, data::AbstractVector)

Iterates backwards through the `data` and fills missing data with the next
existing observation.

WARNING: missing elements at the tail of the array may not be imputed if there is no
existing observation to carry backward. As a result, this method does not guarantee
that all missing values will be imputed.

# Usage
```

```
"""
function impute!(imp::NOCB, ctx::Context, data::AbstractVector)
    end_idx = findlast(ctx, data) - 1
    for i in end_idx:-1:1
        if ismissing(ctx, data[i])
            data[i] = data[i+1]
        end
    end

    return data
end


text:
"""
Measure the memory bandwidth using the following streaming kernels with corresponding data access pattern (Notation: S - store, L - load, WA - write allocate). All variables are vectors, `s` is a scalar:

    init (S1, WA): Initilize an array: `a = s`. Store only.
    sum (L1): Vector reduction: `s += a`. Load only.
    copy (L1, S1, WA): Classic memcopy: `a = b`.
    update (L1, S1): Update vector: `a = a * scalar`. Also load + store but without write allocate.
    triad (L2, S1, WA): Stream triad: `a = b + c * scalar`.
    daxpy (L2, S1): Daxpy: `a = a + b * scalar`.
    striad (L3, S1, WA): Schoenauer triad: `a = b + c * d`.
    sdaxpy (L3, S1): Schoenauer triad without write allocate: `a = a + b * c`.

Keyword arguments:
* `N` (default: `120_000_000`): length of vectors
* `nthreads` (default: `Threads.nthreads()`): number of Julia threads to use
* `niter` (default: `10`): # of times we repeat the measurement
* `alignment` (default: `64`): array alignment
* `verbose` (default: `false`): print result table + thread information etc.
* `write_allocate` (default: `false`): include write allocate compensation factors
"""
function bwbench(;
    N::Integer=120_000_000,
    niter::Integer=10,
    verbose::Bool=false,
    nthreads::Integer=Threads.nthreads(),
    alignment::Integer=64,
    write_allocate::Bool=false,
)
    # check arguments
    1 ≤ N || throw(ArgumentError("N must be ≥ 1."))
    1 ≤ niter || throw(ArgumentError("niter must be ≥ 1."))
    ispow2(alignment) || throw(ArgumentError("alignment $alignment is not a power of 2"))
    alignment ≥ sizeof(Float64) ||
        throw(ArgumentError("alignment $alignment is not a multiple of $(sizeof(Float64))"))
    1 ≤ nthreads ≤ Threads.nthreads() || throw(
        ArgumentError(
            "nthreads $nthreads must ≥ 1 and ≤ $(Threads.nthreads()). If you want more threads, start Julia with a higher number of threads.",
        ),
    )

    # compute thread_indices
    # use only first `nthreads` threads, i.e. @threads for tid in 1:nthreads + manual splitting
    Nperthread = floor(Int, N / nthreads)
    rest = rem(N, nthreads)
    thread_indices = collect(Iterators.partition(1:N, Nperthread))
    if rest != 0
        # last thread compensates for the nonzero remainder
        thread_indices[end-1] = thread_indices[end-1].start:thread_indices[end].stop
    end

    # (maybe) initialize LIKWID Marker API
    LIKWID_init()

    # allocate data
    a = allocate(Float64, N, alignment)
    b = allocate(Float64, N, alignment)
    c = allocate(Float64, N, alignment)
    d = allocate(Float64, N, alignment)
    scalar = 3.0

    # initialize data in parallel (important for NUMA / first-touch policy)
    @threads :static for tid in 1:nthreads
        @inbounds for i in thread_indices[tid]
            a[i] = 2.0
            b[i] = 2.0
            c[i] = 0.5
            d[i] = 1.0
        end
    end

    # print information
    if verbose
        nthreads > 1 && println("Threading enabled, using $nthreads (of $(Threads.nthreads())) Julia threads")
        alloc = 4.0 * sizeof(Float64) * N * 1.0e-06
        println("Total allocated datasize: $(alloc) MB")

        # if Sys.islinux()
        #     @threads :static for i in 1:nthreads
        #         println("\tThread $i running on core $(get_core_id()).")
        #     end
        # end
    end

    # (maybe) register LIKWID Marker regions
    @threads :static for i in 1:nthreads
        LIKWID_register("INIT")
        LIKWID_register("COPY")
        LIKWID_register("UPDATE")
        LIKWID_register("TRIAD")
        LIKWID_register("DAXPY")
        LIKWID_register("STRIAD")
        LIKWID_register("SDAXPY")
    end

    # perform measurement
    times = _run_benchmark(init_kernel, copy_kernel, update_kernel, triad_kernel, daxpy_kernel, striad_kernel, sdaxpy_kernel,
                            a, b, c, d, scalar; niter=niter, thread_indices=thread_indices, nthreads=nthreads)

    # analysis / table output of results
    results = DataFrame(;
        Function=String[],
        var"Rate (MB/s)"=Float64[],
        var"Rate (MFlop/s)"=Float64[],
        var"Avg time"=Float64[],
        var"Min time"=Float64[],
        var"Max time"=Float64[],
    )
    for j in 1:NBENCH
        # ignore the first run because of compilation
        mintime = @views minimum(times[j, 2:end])
        maxtime = @views maximum(times[j, 2:end])
        avgtime = @views mean(times[j, 2:end])
        if write_allocate
            bytes = BENCHMARKS[j].words * BENCHMARKS[j].write_alloc_factor * sizeof(Float64) * N
        else
            bytes = BENCHMARKS[j].words * sizeof(Float64) * N
        end
        flops = BENCHMARKS[j].flops * N
        data_rate = 1.0e-06 * bytes / mintime
        flop_rate = 1.0e-06 * flops / mintime
        push!(
            results, [BENCHMARKS[j].label, data_rate, flop_rate, avgtime, mintime, maxtime]
        )
    end
    verbose && pretty_table(results)

    # validation
    validate(a, b, c, d, N, niter)

    # (maybe) close LIKWID Marker API
    LIKWID_close()
    return results
end

function _run_benchmark(init_kernel, copy_kernel, update_kernel, triad_kernel, daxpy_kernel, striad_kernel, sdaxpy_kernel,
                        a, b, c, d, scalar; niter, nthreads, kwargs...)
    times = zeros(NBENCH, niter)
    for k in 1:niter
        LIKWID_start("INIT", nthreads)
        times[1, k] = @elapsed init_kernel(b, scalar; nthreads=nthreads, kwargs...)
        LIKWID_stop("INIT", nthreads)

        LIKWID_start("COPY", nthreads)
        times[2, k] = @elapsed copy_kernel(c, a; nthreads=nthreads, kwargs...)
        LIKWID_stop("COPY", nthreads)

        LIKWID_start("UPDATE", nthreads)
        times[3, k] = @elapsed update_kernel(a, scalar; nthreads=nthreads, kwargs...)
        LIKWID_stop("UPDATE", nthreads)

        LIKWID_start("TRIAD", nthreads)
        times[4, k] = @elapsed triad_kernel(a, b, c, scalar; nthreads=nthreads, kwargs...)
        LIKWID_stop("TRIAD", nthreads)

        LIKWID_start("DAXPY", nthreads)
        times[5, k] = @elapsed daxpy_kernel(a, b, scalar; nthreads=nthreads, kwargs...)
        LIKWID_stop("DAXPY", nthreads)

        LIKWID_start("STRIAD", nthreads)
        times[6, k] = @elapsed striad_kernel(a, b, c, d; nthreads=nthreads, kwargs...)
        LIKWID_stop("STRIAD", nthreads)

        LIKWID_start("SDAXPY", nthreads)
        times[7, k] = @elapsed sdaxpy_kernel(a, b, c; nthreads=nthreads, kwargs...)
        LIKWID_stop("SDAXPY", nthreads)
    end
    return times
end

function validate(a, b, c, d, N, niter)
    # reproduce initialization
    aj = 2.0
    bj = 2.0
    cj = 0.5
    dj = 1.0

    # now execute timing loop
    scalar = 3.0

    for k in 1:niter
        bj = scalar
        cj = aj
        aj = aj * scalar
        aj = bj + scalar * cj
        aj = aj + scalar * bj
        aj = bj + cj * dj
        aj = aj + bj * cj
    end

    aj = aj * N
    bj = bj * N
    cj = cj * N
    dj = dj * N

    asum = 0.0
    bsum = 0.0
    csum = 0.0
    dsum = 0.0

    for i in 1:N
        asum += a[i]
        bsum += b[i]
        csum += c[i]
        dsum += d[i]
    end

    epsilon = 1.e-8

    if abs(aj - asum) / asum > epsilon
        @warn("Failed Validation on array a[]: Expected $aj, Observed $asum")
        return false
    elseif abs(bj - bsum) / bsum > epsilon
        @warn("Failed Validation on array b[]: Expected $bj, Observed $bsum")
        return false
    elseif abs(cj - csum) / csum > epsilon
        @warn("Failed Validation on array c[]: Expected $cj, Observed $csum")
        return false
    elseif abs(dj - dsum) / dsum > epsilon
        @warn("Failed Validation on array d[]: Expected $dj, Observed $dsum")
        return false
    else
        return true
    end
end

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
import StatsFuns.logsumexp
import Statistics.var

logmeanexp(a) = logsumexp(a) - log(length(a))

function testEstimator(estimator, trials)
  vals::Vector{Float64} = Vector{Float64}(undef, trials)
  totalTime = @elapsed for trial in 1:trials
    vals[trial] = estimator()
  end
  m::Float64 = logmeanexp(vals)
  v::Float64 = var(exp.(vals .- m))
  println("time per trial = ", totalTime / trials )
  println("log(mean) ≈ ", m, ", relative variance ≈ ", v)
end


text:
##############################################################################
##
## Allocations for AllocatedDogleg
##
##############################################################################

type AllocatedDogleg{Tx1, Tx2, Tx3, Tx4, Ty1, Ty2} <: AbstractAllocatedOptimizer
    δgn::Tx1
    δgr::Tx2
    δx::Tx3
    dtd::Tx4
    ftrial::Ty1
    fpredict::Ty2

    function AllocatedDogleg(δgn, δgr, δx, dtd, ftrial, fpredict)
        length(δgn) == length(δgr) || throw(DimensionMismatch("The lengths of δgn and δgr must match."))
        length(δgn) == length(δx) || throw(DimensionMismatch("The lengths of δgn and δx must match."))
        length(δgn) == length(dtd) || throw(DimensionMismatch("The lengths of δgn and dtd must match."))
        length(ftrial) == length(fpredict) || throw(DimensionMismatch("The lengths of ftrial and fpredict must match."))
        new(δgn, δgr, δx, dtd, ftrial, fpredict)
    end
end

function AllocatedDogleg{Tx1, Tx2, Tx3, Tx4, Ty1, Ty2}(δgn::Tx1, δgr::Tx2, δx::Tx3, dtd::Tx4, ftrial::Ty1, fpredict::Ty2)
    AllocatedDogleg{Tx1, Tx2, Tx3, Tx4, Ty1, Ty2}(δgn, δgr, δx, dtd, ftrial, fpredict)
end

function AbstractAllocatedOptimizer{Tx, Ty, Tf, TJ, Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::Dogleg)
    AllocatedDogleg(_zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x),
    _zeros(nls.y), _zeros(nls.y))
end

##############################################################################
##
## Method for AllocatedDogleg
##
##############################################################################
macro Allocateddoglegtrace()
    quote
        if tracing
            update!(tr,
                    iter,
                    ssr,
                    maxabs_gr,
                    store_trace,
                    show_trace,
                    show_every)
        end
    end
end

const MIN_Δ = 1e-16 # maximum trust region radius
const MAX_Δ = 1e16 # minimum trust region radius
const MIN_STEP_QUALITY = 1e-3
const MIN_DIAGONAL = 1e-6
const MAX_DIAGONAL = 1e32
const DECREASE_THRESHOLD = 0.25
const INCREASE_THRESHOLD = 0.75

function optimize!{Tx, Ty, Tf, TJ, Tg, Toptimizer <: AllocatedDogleg, Tsolver}(
    anls::LeastSquaresProblemAllocated{Tx, Ty, Tf, TJ, Tg, Toptimizer, Tsolver};
    xtol::Number = 1e-8, ftol::Number = 1e-8, grtol::Number = 1e-8,
    iterations::Integer = 1_000, Δ::Number = 1.0, store_trace = false, show_trace = false, show_every = 1)

     δgn, δgr, δx, dtd = anls.optimizer.δgn, anls.optimizer.δgr, anls.optimizer.δx, anls.optimizer.dtd
     ftrial, fpredict = anls.optimizer.ftrial, anls.optimizer.fpredict
     x, fcur, f!, J, g! = anls.x, anls.y, anls.f!, anls.J, anls.g!

    # initialize
    reuse = false
    f_calls,  g_calls, mul_calls = 0, 0, 0
    converged, x_converged, f_converged, gr_converged, converged =
        false, false, false, false, false
    f!(x, fcur)
    f_calls += 1
    ssr = sumabs2(fcur)
    maxabs_gr = Inf

    iter = 0
    eTx, eTy = eltype(x), eltype(fcur)

    tr = OptimizationTrace()
    tracing = store_trace || show_trace
    @Allocateddoglegtrace
    while !converged && iter < iterations
        iter += 1
        check_isfinite(x)
        # compute step
        if !reuse
            #update gradient
            g!(x, J)
            g_calls += 1
            colsumabs2!(dtd, J)
            clamp!(dtd, MIN_DIAGONAL, MAX_DIAGONAL)

            if iter == 1
                wnorm_x = wnorm(x, dtd)
                if wnorm_x > 0
                    Δ *= wnorm_x
                end
            end
            # compute (opposite) gradient
            Ac_mul_B!(one(eTx), J, fcur, zero(eTx), δgr)
            mul_calls += 1
            maxabs_gr = maxabs(δgr)
            wnorm_δgr = wnorm(δgr, dtd)

            # compute Cauchy point
            map!((x, y) -> x * sqrt(y), δgn, δgr, dtd)
            A_mul_B!(one(eTy), J, δgn, zero(eTy), fpredict)
            mul_calls += 1
            α = wnorm_δgr^2 / sumabs2(fpredict)

            # compute Gauss Newton step δgn
            fill!(δgn, zero(eTx))
            δgn, ls_iter = A_ldiv_B!(δgn, J, fcur, anls.solver)
            mul_calls += ls_iter
            wnorm_δgn = wnorm(δgn, dtd)
        end
        # compute δx
        if wnorm_δgn <= Δ
            #  Case 1. The Gauss-Newton step lies inside the trust region
            copy!(δx, δgn)
            wnorm_δx = wnorm_δgn
        elseif wnorm_δgr * α >= Δ
            # Case 2. The Cauchy point and the Gauss-Newton steps lie outside
            # the trust region.
            # rescale Cauchy step within the trust region and return
            copy!(δx, δgr)
            scale!(δx, Δ / wnorm_δgr)
            wnorm_δx = Δ
        else
            # Case 3. The Cauchy point is inside the trust region nd the
            # Gauss-Newton step is outside
            b_dot_a = α * wdot(δgr, δgn, dtd)
            a_squared_norm = (α * wnorm_δgr)^2
            b_minus_a_squared_norm =
                  a_squared_norm - 2 * b_dot_a + wnorm_δgn^2
            c =  b_dot_a - a_squared_norm
            d = sqrt(c^2 + b_minus_a_squared_norm * (Δ^2 - a_squared_norm))
            β = (c <= 0) ? (d - c)/ b_minus_a_squared_norm : (Δ^2 - a_squared_norm) / (d + c)
            copy!(δx, δgn)
            scale!(δx, β)
            axpy!(α * (1 - β), δgr, δx)
            wnorm_δx = wnorm(δx, dtd)
        end

        #update x
        axpy!(-one(eTx), δx, x)
        f!(x, ftrial)
        f_calls += 1

        # trial ssr
        trial_ssr = sumabs2(ftrial)

        # predicted ssr
        A_mul_B!(one(eTx), J, δx, zero(eTx), fpredict)
        mul_calls += 1
        axpy!(-one(eTy), fcur, fpredict)
        predicted_ssr = sumabs2(fpredict)

        ρ = (ssr - trial_ssr) / (ssr - predicted_ssr)
        x_converged, f_converged, gr_converged, converged =
            assess_convergence(δx, x, maxabs_gr, ssr, trial_ssr, xtol, ftol, grtol)

        if ρ >= MIN_STEP_QUALITY
            # Successful iteration
            reuse = false
            copy!(fcur, ftrial)
            ssr = trial_ssr
        else
            # unsucessful iteration
            reuse = true
            axpy!(one(eTx), δx, x)
        end

        if ρ < DECREASE_THRESHOLD
           Δ = max(MIN_Δ, Δ * 0.5)
        elseif ρ > INCREASE_THRESHOLD
           Δ = max(Δ, 3.0 * wnorm_δx)
       end
       @Allocateddoglegtrace
    end
    LeastSquaresResult("Dogleg", x, ssr, iter, converged,
                        x_converged, xtol, f_converged, ftol, gr_converged, grtol, tr,
                        f_calls, g_calls, mul_calls)
end


text:
#=
---------------------------------------------------------------------
This function computes the integer arrays FACE, mapL, and mapR, which
allows us to loop through the faces of the grid in order to perform flux integrals.

Written by F.X. Giraldo on July 14, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

function create_face(iside,intma,Nface,Np)

    #global arrays
    face=zeros(Int64,4,Nface)
    mapL=zeros(Int64,2,Np,4)
    mapR=zeros(Int64,2,Np,4)

    #local arrays
    inode=zeros(Int64,4,1)
    jnode=zeros(Int64,4,1)

    #Construct Boundary Pointer
    inode[1]=1
    inode[2]=Np
    inode[3]=Np
    inode[4]=1
    jnode[1]=1
    jnode[2]=1
    jnode[3]=Np
    jnode[4]=Np

    #Construct IMAP arrays
    for l=1:Np

       #eta=-1
       mapL[1,l,1]=l
       mapL[2,l,1]=1
       mapR[1,l,1]=Np+1-l
       mapR[2,l,1]=1

       #ksi=+1
       mapL[1,l,2]=Np
       mapL[2,l,2]=l
       mapR[1,l,2]=Np
       mapR[2,l,2]=Np+1-l

       #eta=+1
       mapL[1,l,3]=Np+1-l
       mapL[2,l,3]=Np
       mapR[1,l,3]=l
       mapR[2,l,3]=Np

       #ksi=-1
       mapL[1,l,4]=1
       mapL[2,l,4]=Np+1-l
       mapR[1,l,4]=1
       mapR[2,l,4]=l
    end #l

    #loop thru the sides
    for i=1:Nface

       ip1=iside[1,i]
       ip2=iside[2,i]
       iel=iside[3,i]
       ier=iside[4,i]

       #check for position on Left Element
       for j=1:4
          j1=j
          j2=j+1
          if (j2 > 4)
             j2=1
          end #j2
          jp1=intma[inode[j1],jnode[j1],iel]
          jp2=intma[inode[j2],jnode[j2],iel]

          if (ip1 == jp1 && ip2 == jp2)
             face[1,i]=j
             break #leave J loop
          end #ip1
       end #j

       #check for position on Right Element
       if (ier > 0)
          for j=1:4
             j1=j
             j2=j+1
             if (j2 > 4)
                j2=1
             end #j2
             jp1=intma[inode[j1],jnode[j1],ier]
             jp2=intma[inode[j2],jnode[j2],ier]

             if (ip1 == jp2 && ip2 == jp1)
                face[2,i]=j
                break #leave J loop
             end #ip1
          end #j
       end #ier

       #Store Elements into FACE
       face[3,i]=iel
       face[4,i]=ier
    end #i

    return(face,mapL,mapR)

end


text:
export WeightedSoftmaxExplorer

using Random
using StatsBase: sample, Weights
using Flux: softmax

"""
    WeightedSoftmaxExplorer(;rng=Random.GLOBAL_RNG)

See also: [`WeightedExplorer`](@ref)
"""
struct WeightedSoftmaxExplorer{R<:AbstractRNG} <: AbstractExplorer
    rng::R
end

function WeightedSoftmaxExplorer(; rng = Random.GLOBAL_RNG)
    WeightedSoftmaxExplorer(rng)
end

(s::WeightedSoftmaxExplorer)(values::AbstractVector{T}) where {T} =
    sample(s.rng, Weights(softmax(values), one(T)))

function (s::WeightedSoftmaxExplorer)(values::AbstractVector{T}, mask) where {T}
    values[.!mask] .= typemin(T)
    s(values)
end


text:
module Equate
# =======

using DataFrames: DataFrame
using Statistics: mean, std, cov, cor, var
using StatsBase: kurtosis, skewness
using Distributions: Poisson, Normal, pdf, cdf
using Optim: optimize, BFGS
using GLM: glm, predict, @formula, LogLink, coef
using StatsModels: TableRegressionModel
using StatsBase
import StatsBase: coef
using Printf
import Base: copy, show
import Bootstrap: draw!

export
    # Basic function
    freqtab,

    # Equate functions
    Equipercentile,
    Mean,
    Linear,
    Tucker,
    LevineCongeneric,
    LevineObservedScore,
    LevineTrueScore,
    ChainedLinear,
    FrequencyEstimation,
    BraunHolland,
    ChainedEquipercentile,

    # Smoothing functions
    presmoothing,
    KernelSmoothing,

    # Utility functions
    coef,
    ExpandTable,
    expandtable,
    SummaryStats,

    # Support functions
    round2,
    LogLinearFormula,
    ObservableStats,
    BandwidthPenalty,
    EstBandwidth,

    # SEE
    recalculate!,
    initialize!,
    copy,
    show,
    draw!,

    # struct
    NEATEquateResult,
    SGEquateResult

abstract type EquateDesign end
"""
    EG
Equivalent (single) group design.
If the population of 2 tests is identical, the test design will be called the single group design.
"""
abstract type EG <: EquateDesign end
"""
    NEAT
Non-Equivalent group with Anchor Test design.
"""
abstract type NEAT <: EquateDesign end
abstract type EquateMethod end
abstract type SGEquateMethod <: EquateMethod end
abstract type NEATEquateMethod <: EquateMethod end

include("freqtab.jl")
include("SGdesign.jl")
include("NEATdesign.jl")
include("Smoothing.jl")
include("coef.jl")
include("ExpandTable.jl")
include("show.jl")
include("SEE.jl")

"""
    A Julia package for test equating.

Major features:

*Single Group (SG)* design

- `Linear` provides the linear equating.
- `Equipercentile` provides the equipercentile equating.

*Non Equivalent group with Anchor Test (NEAT) design*

- `Tucker`
- `ChainedLinear`
- `ChainedEquipercentile`
- `BraunHolland`
- `FrequencyEstimation`

and more...

According to von Davier, Holland, and Thayer (2004), the test equating has separate 5 steps.

- First step, *Pre-smoothing*.
- Second step is *Estimating the score probabilities*.
- Third step is *Continuization*.
- Fourth step is *Equating*.
- Fifth, last step is *Calculating the Standard Error of Equating(SEE)*.

"""
Equate
#-----------------
end # module


text:
import DandelionSlack.post
import Base.==

facts("WebApi query") do
    @fact toquery(RtmStart(Token(utf8("foo")), utf8("1"), utf8("2"), utf8("3"))) --> Dict(
        "token" => "foo",
        "simple_latest" => "1",
        "no_unreads" => "2",
        "mpim_aware" => "3")

    @fact toquery(RtmStart(Token(utf8("foo")), utf8("1"), utf8("2"), Nullable{UTF8String}())) --> Dict(
        "token" => "foo",
        "simple_latest" => "1",
        "no_unreads" => "2")

    @fact toquery(RtmStart(Token(utf8("foo")), utf8("1"), Nullable{UTF8String}(), Nullable{UTF8String}())) -->
        Dict("token" => "foo", "simple_latest" => "1")
end

# Define a fake method, so we don't have to use an actual Slack method,
# because they're complex.
immutable FakeMethod
    token::UTF8String
end

immutable FakeMethodResponse
    url::UTF8String
end

DandelionSlack.getresponsetype(::Type{FakeMethod}) = FakeMethodResponse
DandelionSlack.method_name(::Type{FakeMethod}) = "fake.method"

==(a::FakeMethodResponse, b::FakeMethodResponse) = a.url == b.url

# A mock response from Slack.
type MockHttpResponse <: AbstractHttpResponse
    status::Int
    text::AbstractString
end

# A test case defines a request, a mock http response, and an expected status
# and deserialized response. The test method executes the request and checks
# actual response against the expected.
type MethodTestCase
    test::AbstractString
    request::FakeMethod
    http_response::MockHttpResponse
    expected_status::Status
    expected_response::Nullable{FakeMethodResponse}
end

# A mock HTTP object. The "post" method defined on it mocks the Request post method.
type MockHttp <: DandelionSlack.AbstractHttp
    rv::Any
end

function post(m::MockHttp, url::AbstractString; query=Dict{AbstractString, Any}())
    m.rv
end

function test_successful_method(tc::MethodTestCase)
    mock_http = MockHttp(tc.http_response)
    @fact DandelionSlack.makerequest(tc.request, mock_http) -->
        (tc.expected_status, get(tc.expected_response))
end

function test_method_error(tc::MethodTestCase)
    mock_http = MockHttp(tc.http_response)
    @fact_throws DandelionSlack.RequestException DandelionSlack.makerequest(tc.request, mock_http)
end

function test_http_call_error(tc::MethodTestCase)
    mock_http = MockHttp(tc.http_response)
    @fact_throws DandelionSlack.HttpException DandelionSlack.makerequest(tc.request, mock_http)
end

facts("Slack method tests") do
    testcase_dir = "test/testcases"
    testcase_dir_entries = readdir(testcase_dir)
    istestcase = x -> isfile(joinpath(testcase_dir, string(x))) && endswith(x, ".json")
    testcases = filter(istestcase, testcase_dir_entries)

    for testcase_file in testcases
        context("Test case $(testcase_file)") do
            tc_json = JSON.parsefile(joinpath(testcase_dir, string(testcase_file)))
            tc = DandelionSlack.deserialize(MethodTestCase, tc_json)

            if tc.test == "MethodCall"
                test_successful_method(tc)
            elseif tc.test == "MethodCallError"
                test_method_error(tc)
            elseif tc.test == "HttpCallError"
                test_http_call_error(tc)
            else
                throws(Exception("Unexpected test type $(tc.test)"))
            end
        end
    end
end


text:
# Computes the SVD of a 2-by-2 matrix M, returning U, S, and V such that A == U*S*V'.
# Reference:
# Blinn, J. (1996). Consider the lowly 2 x 2 matrix. IEEE Computer Graphics and Applications, 16(2), 82-88.
# https://scicomp.stackexchange.com/questions/8899/robust-algorithm-for-2-times-2-svd
function svd2x2(M::AbstractArray)
    E = (M[1,1] + M[2,2]) / 2
    F = (M[1,1] - M[2,2]) / 2
    G = (M[2,1] + M[1,2]) / 2
    H = (M[2,1] - M[1,2]) / 2
    Q = sqrt(E^2 +H^2)
    R = sqrt(F^2 + G^2)
    sx = Q + R
    sy = Q - R
    a₁ = atan2(G,F)
    a₂ = atan2(H,E)
    θ = (a₂ - a₁) / 2
    ϕ = (a₂ + a₁) / 2
    s = sign(sy)
    U = SMatrix{2,2,Float64}(cos(ϕ), sin(ϕ), -s*sin(ϕ), s*cos(ϕ))
    S = SMatrix{2,2,Float64}(sx, 0.0, 0.0, abs(sy))
    V = SMatrix{2,2,Float64}(cos(θ), -sin(θ), sin(θ), cos(θ))
    U,S,V
end

#  Computes the Moore-Penrose pseudoinverse for a 2-by-2 matrix M by only inverting
#  singular values above the threshold `tol = sqrt(eps(real(float(one(eltype(M))))))`
function pinv2x2(M::AbstractArray)
     tol = sqrt(eps(real(float(one(eltype(M))))))
     U,S,V = svd2x2(M)
     D = SMatrix{2,2,Float64}( S[1,1] > tol ?  1/S[1,1] : 0.0 , 0.0, 0.0, S[2,2] > tol ?  1/S[2,2] : 0.0 )
     U*D*V'
end


text:
start_init = time()

include("distributed.jl")
#add_procs() #can be restored after package registration

start_pkg = time()

#@everywhere using Pkg
#@everywhere Pkg.activate(".")

@everywhere using PowerModelsSecurityConstrained

include("goc_challenge1_huristic.jl")

println("package load time: $(time() - start_pkg)")

println("script startup time: $(time() - start_init)")

function MyJulia1(InFile1::String, InFile2::String, InFile3::String, InFile4::String, TimeLimitInSeconds::Int64, ScoringMethod::Int64, NetworkModel::String)
    println("running MyJulia1")
    println("  $(InFile1)")
    println("  $(InFile2)")
    println("  $(InFile3)")
    println("  $(InFile4)")
    println("  $(TimeLimitInSeconds)")
    println("  $(ScoringMethod)")
    println("  $(NetworkModel)")

    startup_time = 60

    compute_solution1(InFile1, InFile2, InFile3, InFile4, TimeLimitInSeconds-startup_time, ScoringMethod, NetworkModel)
end


text:
@inline loc_lower(xi, Δx, n, x) = @inbounds max(min(unsafe_trunc(Int, (xi - x[1]) / Δx) + 1, n - 1), 1)

"""
    interp1d(xi::FT, x, y, col) where {FT<:AbstractFloat}

perform 1D linear interpolation.
"""
@inline function interp1d(xi::FT, x, y, col) where {FT <: AbstractFloat}
    @inbounds Δx = x[2] - x[1]
    n = length(x)
    loc = loc_lower(xi, Δx, n, x)
    @inbounds factor = (xi - x[loc]) / Δx
    return @inbounds (y[loc, col] * (FT(1) - factor) + y[loc + 1, col] * factor)
end

"""
    interp2d(
        fη1::FT,
        fη2::FT,
        ftemp::FT,
        coeff::FTA3D,
        igpt::I,
        jη1::I,
        jη2::I,
        jtemp::I,
    ) where {FT<:AbstractFloat,FTA3D<:AbstractArray{FT,3},I<:Int}

Perform 2D linear interpolation.

`fminor[1, 1] = (FT(1) - fη1) * (1-ftemp)`

`fminor[2, 1] = fη1 * (1-ftemp)`

`fminor[1, 2] = (FT(1) - fη2) * ftemp`

`fminor[2, 2] = fη2 * ftemp`
"""
@inline function interp2d(
    fη1::FT,
    fη2::FT,
    ftemp::FT,
    coeff::FTA3D,
    igpt::I,
    jη1::I,
    jη2::I,
    jtemp::I,
) where {FT <: AbstractFloat, FTA3D <: AbstractArray{FT, 3}, I <: Int}
    return @inbounds (FT(1) - fη1) * (1 - ftemp) * coeff[igpt, jη1, jtemp] +
                     fη1 * (1 - ftemp) * coeff[igpt, jη1 + 1, jtemp] +
                     (FT(1) - fη2) * ftemp * coeff[igpt, jη2, jtemp + 1] +
                     fη2 * ftemp * coeff[igpt, jη2 + 1, jtemp + 1]
end

"""
    interp3d(
        jη1::I,
        jη2::I,
        fη1::FT,
        fη2::FT,
        jtemp::I,
        ftemp::FT,
        jpresst::I,
        fpress::FT,
        coeff::FTA4D,
        igpt::I,
        s1::FT = FT(1),
        s2::FT = FT(1),
    ) where {FT<:AbstractFloat,FTA4D<:AbstractArray{FT,4},I<:Int}

Perform 3D linear interpolation.

`fmajor[1, 1, itemp] = (FT(1) - fpress) * fminor[1, itemp]`

`fmajor[2, 1, itemp] = (FT(1) - fpress) * fminor[2, itemp]`

`fmajor[1, 2, itemp] = fpress * fminor[1, itemp]`

`fmajor[2, 2, itemp] = fpress * fminor[2, itemp]`

where,

`fminor[1, itemp=1] = (FT(1) - fη1) * (1-ftemp)`

`fminor[2, itemp=1] = fη1 * (1-ftemp)`

`fminor[1, itemp=2] = (FT(1) - fη2) * ftemp`

`fminor[2, itemp=2] = fη2 * ftemp`

"""
@inline function interp3d(
    jη1::I,
    jη2::I,
    fη1::FT,
    fη2::FT,
    jtemp::I,
    ftemp::FT,
    jpresst::I,
    fpress::FT,
    coeff::FTA4D,
    igpt::I,
    s1::FT = FT(1),
    s2::FT = FT(1),
) where {FT <: AbstractFloat, FTA4D <: AbstractArray{FT, 4}, I <: Int}
    omftemp = FT(1) - ftemp
    omfpress = FT(1) - fpress
    omfη1 = FT(1) - fη1
    omfη2 = FT(1) - fη2
    return @inbounds s1 * (
        omfpress *
        (omftemp * (omfη1 * coeff[igpt, jη1, jpresst - 1, jtemp] + fη1 * coeff[igpt, jη1 + 1, jpresst - 1, jtemp])) +
        fpress * (omftemp * (omfη1 * coeff[igpt, jη1, jpresst, jtemp] + fη1 * coeff[igpt, jη1 + 1, jpresst, jtemp]))
    ) +
                     s2 * (
        omfpress * (
            ftemp *
            (omfη2 * coeff[igpt, jη2, jpresst - 1, jtemp + 1] + fη2 * coeff[igpt, jη2 + 1, jpresst - 1, jtemp + 1])
        ) +
        fpress *
        (ftemp * (omfη2 * coeff[igpt, jη2, jpresst, jtemp + 1] + fη2 * coeff[igpt, jη2 + 1, jpresst, jtemp + 1]))
    )
end
"""
    increment!(
        op::TwoStream{FT},
        τ2,
        ssa2,
        g2,
        glay,
        gcol,
        igpt,
    ) where {FT<:AbstractFloat}

Increment TwoStream optical properties `op` for layer `glay`
and column `gcol` with optical thickness `τ2`, single scattering
albedo `ssa2` and symmetry parameter `g2`.
"""
function increment!(op::TwoStream{FT}, τ2, ssa2, g2, glay, gcol, igpt) where {FT <: AbstractFloat}
    τ1, ssa1, g1 = op.τ[glay, gcol], op.ssa[glay, gcol], op.g[glay, gcol]
    τ = τ1 + τ2
    ssa = τ1 * ssa1 + τ2 * ssa2
    ssag = (τ1 * ssa1 * g1 + τ2 * ssa2 * g2) / max(eps(FT), ssa)
    ssa /= max(eps(FT), τ)

    op.τ[glay, gcol] = τ
    op.ssa[glay, gcol] = ssa
    op.g[glay, gcol] = ssag
    return nothing
end
"""
    delta_scale(τ, ssa, g)

delta-scale two stream optical properties.
"""
function delta_scale(τ, ssa, g)
    FT = typeof(τ)
    f = g * g
    wf = ssa * f
    τ_s = (FT(1) - wf) * τ
    ssa_s = (ssa - wf) / max(eps(FT), FT(1) - wf)
    g_s = (g - f) / max(eps(FT), FT(1) - f)
    return (τ_s, ssa_s, g_s)
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
export get_belief,get_qval!,Policy,read_policy,evaluate

type Policy
    alpha       :: Matrix{Float64}
    actions     :: Matrix{Float64}
    nactions    :: Int64
    qvals       :: Vector{Float64}

    function Policy(alpha::Matrix{Float64}, actions::Matrix{Float64})
        return new(alpha, actions, size(actions, 2), zeros(size(actions, 2)))
    end # function Policy
end

function read_policy(actions::Matrix{Float64}, alpha::Matrix{Float64})
    return Policy(alpha, actions)
end # function read_policy

function evaluate(policy::Policy, belief::SparseMatrixCSC{Float64,Int64})
    fill!(policy.qvals, 0.0)
    get_qval!(policy, belief)
    return copy(policy.qvals)
end # function evaluate

function get_qval!(policy::Policy, belief::SparseMatrixCSC{Float64, Int64})
    fill!(policy.qvals, 0.0)
    for iaction in 1:policy.nactions
        for ib in 1:length(belief.rowval)
            policy.qvals[iaction] += belief.nzval[ib] * policy.alpha[belief.rowval[ib], iaction]
        end # for b
    end # for iaction
    #println(policy.qvals)
end # function get_qval!

function get_belief(pstate::Vector{Float64}, grid::RectangleGrid,interp::Bool=false,drl::Bool=false,XandY::Bool=false)
    belief = spzeros(NSTATES, 1)
    if drl
        if XandY
            belief = spzeros(NSTATES_drl_xandy,1)
        else
            belief = spzeros(NSTATES_drl,1)
        end
    end
    indices, weights = interpolants(grid, pstate)
    if !interp
        largestWeight = 0;
        largestIndex = 0;
        for i = 1:length(weights)
            if weights[i]>largestWeight
                largestWeight = weights[i]
                largestIndex = indices[i]
            end
        end
        indices = largestIndex
        weights = 1.0
    end
    for i = 1:length(indices)
        belief[indices[i]] = weights[i]
    end # for i
    return belief
end # function get_belief

function belief_states(r,th,psi_int,v_own,v_int,tau,pa,deltaR, deltaPsi,deltaV, deltaTheta,nnet)
    belief = zeros(9,7)
    pastActions = pas
    if nnet
        pastActions = pasTrue
    end
    belief[1,:] = [r,th,deg2rad(psi_int),v_own,v_int,tau,pastActions[pa]];
    belief[2,:] = [r+deltaR,th,deg2rad(psi_int),v_own,v_int,tau,pastActions[pa]];

    rTemp = r-deltaR
    if rTemp < 0
        rTemp = 0
    end
    belief[3,:] = [rTemp,th,deg2rad(psi_int),v_own,v_int,tau,pastActions[pa]];

    psiTemp = psi_int+deltaPsi
    if psiTemp>180.0
        psiTemp-=360.0
    end
    belief[4,:] = [r,th,deg2rad(psiTemp),v_own,v_int,tau,pastActions[pa]];

    psiTemp = psi_int-deltaPsi
    if psiTemp<-180.0
        psiTemp+=360.0
    end
    belief[5,:] = [r,th,deg2rad(psiTemp),v_own,v_int,tau,pastActions[pa]];
    belief[6,:] = [r,th,deg2rad(psi_int),v_own,v_int+deltaV,tau,pastActions[pa]];

    vTemp = v_int-deltaV
    if vTemp < 0
        vTemp = 0
    end
    belief[7,:] = [r,th,deg2rad(psi_int),v_own,vTemp,tau,pastActions[pa]];

    thTemp = th+deltaTheta
    if thTemp>pi
        thTemp-=2*pi
    end
    belief[8,:] = [r,thTemp,deg2rad(psi_int),v_own,v_int,tau,pastActions[pa]];

    thTemp = th-deltaTheta
    if thTemp<-pi
        thTemp+=2*pi
    end
    belief[9,:] = [r,thTemp,deg2rad(psi_int),v_own,v_int,tau,pastActions[pa]];

    return belief

end #function belief_states


text:
"""
    GreenNode(head, span)
    GreenNode(head, children...)

A "green tree" in Roslyn (C# compiler) terminology is a lossless syntax tree
which overlays all the source text. The most basic properties of a green tree
are that:

* Nodes cover a contiguous span of bytes in the text
* Sibling nodes are ordered in the same order as the text

As implementation choices, we choose that:

* Nodes are immutable and don't know their parents or absolute position, so can
  be cached and reused
* Nodes are homogenously typed at the language level so they can be stored
  concretely, with the `head` defining the node type. Normally this would
  include a "syntax kind" enumeration, but it can also include flags and record
  information the parser knew about the layout of the child nodes.
* For simplicity and uniformity, leaf nodes cover a single token in the source.
  This is like rust-analyzer, but different from Roslyn where leaves can
  include syntax trivia.

Design principles:
* Tree should remember what the lexer and parser knew about the source code
* Be position-independent so nodes can be interned and reused
* Be a low level textural overlay which is language independent.

Design alternatives to explore:
* Maybe allow some loss of local parser state if it can be derived again
  quickly? Particularly in the ordering of children.
* Store strings for tokens? (Surprisingly, rust-analyzer does this. It could be
  efficient if the strings or nodes are interned for the parsing session?)
* Never construct this tree? Instead serialize it to Vector{UInt8} in an
  efficient but compact format? Could this be more flexible with storing parser
  state and beat the interning approach? We could also store the source tokens
  in the serialization and discard the source text. (Caveat - unclear that this
  could deal with incremental parsing...)
"""
struct GreenNode{Head}
    head::Head
    span::UInt32
    args::Union{Tuple{},Vector{GreenNode{Head}}}
end

function GreenNode(head::Head, span::Integer) where {Head}
    GreenNode{Head}(head, span, ())
end

function GreenNode(head::Head, span::Integer, args::Vector{GreenNode{Head}}) where {Head}
    GreenNode{Head}(head, span, args)
end

function GreenNode(head::Head, args::GreenNode{Head}...) where {Head}
    span = sum(x.span for x in args)
    GreenNode{Head}(head, span, GreenNode{Head}[args...])
end

# Accessors / predicates
haschildren(node::GreenNode) = !(node.args isa Tuple{})
children(node::GreenNode)    = node.args
span(node::GreenNode)        = node.span
head(node::GreenNode)        = node.head

# Predicates
is_trivia(node::GreenNode) = is_trivia(node.head)
is_error(node::GreenNode)  = is_error(node.head)

Base.summary(node::GreenNode) = summary(node.head)

# Pretty printing
function _show_green_node(io, node, indent, pos, str, show_trivia)
    if !show_trivia && is_trivia(node)
        return
    end
    posstr = "$(lpad(pos, 6)):$(rpad(pos+span(node)-1, 6)) │"
    is_leaf = !haschildren(node)
    if is_leaf
        line = string(posstr, indent, summary(node))
    else
        line = string(posstr, indent, '[', summary(node), "]")
    end
    if !is_trivia(node) && is_leaf
        line = rpad(line, 40) * "✔"
    end
    if is_error(node)
        line = rpad(line, 41) * "✘"
    end
    if is_leaf && !isnothing(str)
        line = string(rpad(line, 43), ' ', repr(str[pos:prevind(str, pos + span(node))]))
    end
    line = line*"\n"
    if is_error(node)
        printstyled(io, line, color=:light_red)
    else
        print(io, line)
    end
    if !is_leaf
        new_indent = indent*"  "
        p = pos
        for x in children(node)
            _show_green_node(io, x, new_indent, p, str, show_trivia)
            p += x.span
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", node::GreenNode)
    _show_green_node(io, node, "", 1, nothing, true)
end

function Base.show(io::IO, ::MIME"text/plain", node::GreenNode, str::AbstractString; show_trivia=true)
    _show_green_node(io, node, "", 1, str, show_trivia)
end


text:
struct LinearCombinationPotential <: MixedPotential
    potentials::Vector{ArbitraryPotential}
    coefficients::Vector{AbstractFloat}
end
LinearCombinationPotential(p::ArbitraryPotential) = LinearCombinationPotential([p], [1])
LinearCombinationPotential(p::LinearCombinationPotential) = p

Base.:+(p::ArbitraryPotential) = +1 * p
Base.:-(p::ArbitraryPotential) = -1 * p
Base.:+(p1::ArbitraryPotential, p2::ArbitraryPotential) = LinearCombinationPotential(p1) + LinearCombinationPotential(p2)
Base.:-(p1::ArbitraryPotential, p2::ArbitraryPotential) = p1 + -p2
Base.:*(a::Real, p::ArbitraryPotential) = a * LinearCombinationPotential(p)
Base.:/(p::ArbitraryPotential, a::Real) = (1.0 / a) * p

Base.:+(mp1::LinearCombinationPotential, mp2::LinearCombinationPotential) = LinearCombinationPotential(append!(mp1.potentials, mp2.potentials), append!(mp1.coefficients, mp2.coefficients))
Base.:*(a::Real, mp::LinearCombinationPotential) = LinearCombinationPotential(mp.potentials, a * mp.coefficients)

apply_linear_combination(f::Function, sys::AbstractSystem, mp::LinearCombinationPotential) = sum(c * f(sys, p) for (p, c) in zip(mp.potentials, mp.coefficients))

function energy_and_force(sys::AbstractSystem, mp::LinearCombinationPotential)
    nts = [energy_and_force(sys, p) for p in mp.potentials]
    e = sum( c * nt.e for (nt, c) in zip(nts, mp.coefficients) )
    f = sum( c * nt.f for (nt, c) in zip(nts, mp.coefficients) )
    (e = e, f = f)
end

potential_energy(sys::AbstractSystem, mp::LinearCombinationPotential) = apply_linear_combination(potential_energy, sys, mp)
force(sys::AbstractSystem, mp::LinearCombinationPotential) = apply_linear_combination(force, sys, mp)
virial(sys::AbstractSystem, mp::LinearCombinationPotential) = apply_linear_combination(virial, sys, mp)
virial_stress(sys::AbstractSystem, mp::LinearCombinationPotential) = apply_linear_combination(virial_stress, sys, mp)


text:
# check speed of new and old sparse filter

using Revise
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pandas
using Seaborn
using Dates
using Logging

# using Profile
# using ProfileView
using TimerOutputs

const to = TimerOutput()

block_size=[4,4,4]
#
segmentation = LarSurf.generate_segmentation567(20)
LarSurf.grid_get_surf_Fvec_larmodel(segmentation)
segmentation = LarSurf.generate_segmentation567(1)
Flin1, larmodel1 = LarSurf.grid_get_surf_Fvec_larmodel(segmentation)
Flin2, larmodel2 = LarSurf.grid_get_surf_Fvec_larmodel_old(segmentation)
V1, top1 = larmodel1
V2, top2 = larmodel2
(VV1, EV1, FV1, CV1) = top1
(VV2, EV2, FV2, CV2) = top2

Flin_block1 = LarSurf.__grid_get_surface_Fchar_per_block(segmentation, block_size)
Flin_block2 = LarSurf.__grid_get_surface_Fchar_per_block_old_implementation(segmentation, block_size)

for i=1:size(Flin_block1,1)
    if Flin_block1[i] == Flin_block2[i]
    else
        if (Flin_block2[i] == 2) & (Flin_block1[i]==0)
        else
            error("Flin implementation does not fit, $i")
        end
    end
end

@timeit to "size10" begin
    for i=1:5
        @timeit to "full" begin
            @timeit to "new" LarSurf.grid_get_surf_Fvec_larmodel(segmentation)
            @timeit to "old" LarSurf.grid_get_surf_Fvec_larmodel_old(segmentation)
        end
        @timeit to "per block" begin
            @timeit to "new" LarSurf.__grid_get_surface_Fchar_per_block(segmentation, block_size)
            @timeit to "old" LarSurf.__grid_get_surface_Fchar_per_block_old_implementation(segmentation, block_size)
        end
    end
end
segmentation = LarSurf.generate_segmentation567(20)
@timeit to "size20" begin
    for i=1:5
        @timeit to "new" LarSurf.grid_get_surf_Fvec_larmodel(segmentation)
        @timeit to "old" LarSurf.grid_get_surf_Fvec_larmodel_old(segmentation)
    end
end

to
#


text:
# In the card game poker, a hand consists of five cards and are ranked, from
# lowest to highest, in the following way:
#
# High Card: Highest value card.
# One Pair: Two cards of the same value.
# Two Pairs: Two different pairs.
# Three of a Kind: Three cards of the same value.
# Straight: All cards are consecutive values.
# Flush: All cards of the same suit.
# Full House: Three of a kind and a pair.
# Four of a Kind: Four cards of the same value.
# Straight Flush: All cards are consecutive values of same suit.
# Royal Flush: Ten, Jack, Queen, King, Ace, in same suit.
# The cards are valued in the order:
# 2, 3, 4, 5, 6, 7, 8, 9, 10, Jack, Queen, King, Ace.
#
# If two players have the same ranked hands then the rank made up of the highest
# value wins; for example, a pair of eights beats a pair of fives (see example 1
# below). But if two ranks tie, for example, both players have a pair of queens,
# then highest cards in each hand are compared (see example 4 below); if the
# highest cards tie then the next highest cards are compared, and so on.
#
# Consider the following five hands dealt to two players:
#
# Hand  Player 1          Player 2           Winner
# 1     5H 5C 6S 7S KD    2C 3S 8S 8D TD     Player 2
#       Pair of Fives     Pair of Eights
#
# 2     5D 8C 9S JS AC    2C 5C 7D 8S QH     Player 1
#       Highest card Ace  Highest card Queen
#
# 3     2D 9C AS AH AC    3D 6D 7D TD QD     Player 2
#       Three Aces        Flush with Diamonds
#
# 4     4D 6S 9H QH QC    3D 6D 7H QD QS     Player 1
#       Pair of Queens    Pair of Queens
#       Highest card Nine Highest card Seven
#
# 5     2H 2D 4C 4D 4S     3C 3D 3S 9S 9D    Player 1
#       Full House         Full House
#       With Three Fours   with Three Threes
#
# The file, poker.txt, contains one-thousand random hands dealt to two players.
# Each line of the file contains ten cards (separated by a single space): the
# first five are Player 1's cards and the last five are Player 2's cards. You
# can assume that all hands are valid (no invalid characters or repeated cards),
# each player's hand is in no specific order, and in each hand there is a clear
# winner.
#
# How many hands does Player 1 win?

using ProjectEulerSolutions

# Define a unique score for every hand such that a winning hand will score
# higher than a losing hand.
function score_hand(hand::BitArray{2})::Float32
    sum_across_suits = sum(hand, dims=2)[:]
    is_flush = maximum(sum(hand, dims=1)) == 5

    single_index = findall(sum_across_suits .== 1)
    pair_index = findall(sum_across_suits .== 2)

    no_multiple_numbers = maximum(sum_across_suits) == 1
    nonzero_numbers = findall(sum_across_suits .> 0)

    max_value = nonzero_numbers[end]

    if no_multiple_numbers && nonzero_numbers[end] - nonzero_numbers[1] == 4
        if is_flush # straight or royal flush
            return 100000 + max_value
        else # straight
            return 50000 + max_value
        end
    elseif is_flush # flush
        # Check score by largest value, then next largest, etc
        return 60000 + sum([0.0001, 0.01, 1, 14, 14 * 14] .* nonzero_numbers)
    end
    if maximum(sum_across_suits) == 4 # 4 of a kind
        return 80000 + 14 * findall(sum_across_suits .== 4) + single_index[1]
    end
    if maximum(sum_across_suits) == 3
        three_index = findall(sum_across_suits .== 3)
        if length(pair_index) == 1 # Full house
            return 70000 + 14 * three_index[1] + pair_index[1]
        else # 3 of a kind
            return 40000 + 14 * single_index[2] + single_index[1]
        end
    end
    if length(pair_index) == 2 # 2 pair
        return 10000 + sum([14, 14*14] .* pair_index) + single_index[1]
    elseif length(pair_index) == 1 # 1 pair
        return 4000 + 14 * 14 * pair_index[1] + sum([0.01, 1, 14] .* single_index)
    else  # high card
        return sum([0.0001, 0.01, 1, 14, 14 * 14] .* single_index)
    end
end

function p054solution(path::String="")::Integer
    if !isfile(path)
        if path != ""
            println("Cannot find path: '", path, "'")
        end
        return -1
    end

    # Fill in card indexes into hand
    card_inds = Dict{String, Tuple{Integer, Integer}}()
    for (i, vi) in enumerate(['2', '3', '4', '5', '6', '7', '8', '9', 'T', 'J', 'Q', 'K', 'A'])
        for (j, sj) in enumerate(['C', 'S', 'H', 'D'])
            card_inds[vi * sj] = (i, j)
        end
    end

    count = 0
    open(path) do file
        for row in eachline(file)
            inds = map(x -> card_inds[x], split(row, ' '))
            hand1 = falses(13,4)
            hand2 = falses(13,4)
            for i in 1:5
                hand1[inds[i]...] = true
                hand2[inds[i+5]...] = true
            end
            if score_hand(hand1) > score_hand(hand2)
                count += 1
            end
        end
    end

    return count
end

p054 = Problems.Problem(p054solution)

input_path = "data/p054_poker.txt"

Problems.benchmark(p054, input_path)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
module BenchRandomIncrements

using BenchmarkTools
using RecordArrays

const CACHE = Ref{Any}()

@inline function unsafe_inc!(counters::RecordArray, i)
    @inbounds counters.a[i] += 1
end

@inline function unsafe_inc!(counters, i)
    x = @inbounds counters[end]
    x = merge(x, (; a = x.a + 1))
    @inbounds counters[i] = x
end

@noinline function inc!(counters, indices)
    for i in indices
        unsafe_inc!(counters, i)
    end
end

function generate(nindices = 2^10, ncounters = 10)
    return (
        counters_base = fill((a = 0, b = 0, c = 0, d = 0), ncounters),
        counters_rec = RecordArrays.fill((a = 0, b = 0, c = 0, d = 0), ncounters),
        indices = rand(1:ncounters, nindices),
    )
end

function setup()
    CACHE[] = generate()

    suite = BenchmarkGroup()

    suite["base"] = @benchmarkable(
        inc!(counters, indices),
        setup = begin
            counters = CACHE[].counters_base::$(typeof(CACHE[].counters_base))
            indices = CACHE[].indices::$(typeof(CACHE[].indices))
        end,
    )

    suite["rec"] = @benchmarkable(
        inc!(counters, indices),
        setup = begin
            counters = CACHE[].counters_rec::$(typeof(CACHE[].counters_rec))
            indices = CACHE[].indices::$(typeof(CACHE[].indices))
        end,
    )

    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
using Random
using Test

function cubic_hermite_test()
    # Need at least two samples:
    xs = ones(1)
    ys = ones(1)
    dydxs = ones(1)
    @test_throws DomainError CubicHermite(xs, ys, dydxs)

    # All vectors must be the same size:
    xs = ones(2)
    xs[2] = 2
    ys = ones(2)
    dydxs = ones(3)
    @test_throws DomainError CubicHermite(xs, ys, dydxs)

    xs = ones(5)
    for i = 2:length(xs)
        xs[i] = i
    end
    ys = ones(5)
    dydxs = zeros(5)
    ch = CubicHermite(xs, ys, dydxs)
    for i = 1:length(xs)-1
        x = Float64(i)
        @test ch(x) == 1.0
        @test gradient(ch, x) == 0.0
        @test hessian(ch, x) == 0.0
        x = x + 0.25
        @test ch(x) == 1.0
        @test gradient(ch, x) == 0.0
        @test hessian(ch, x) == 0.0
        x = x + 0.25
        @test ch(x) == 1.0
        @test gradient(ch, x) == 0.0
        @test hessian(ch, x) == 0.0
        x = x + 0.25
        @test ch(x) == 1.0
        @test gradient(ch, x) == 0.0
        @test hessian(ch, x) == 0.0
    end
    # Ensure that the ight endpoint query doesn't read past the end of the array:
    @test ch(5.0) == 1.0
    @test gradient(ch, 5.0) == 0.0
    @test hessian(ch, 5.0) == 0.0

    # Now linear functions:
    a = 7.2
    b = 9.6
    for i = 1:length(xs)
        ys[i] = a * xs[i] + b
        dydxs[i] = a
    end
    ch = CubicHermite(xs, ys, dydxs)
    for i = 1:length(xs)-1
        x = Float64(i)
        @test ch(x) ≈ a * x + b
        @test gradient(ch, x) ≈ a
        @test abs(hessian(ch, x)) < 3e-14
        x = x + 0.25
        @test ch(x) ≈ a * x + b
        @test gradient(ch, x) ≈ a
        @test abs(hessian(ch, x)) < 3e-14
        x = x + 0.25
        @test ch(x) ≈ a * x + b
        @test gradient(ch, x) ≈ a
        @test abs(hessian(ch, x)) < 3e-14
        x = x + 0.25
        @test ch(x) ≈ a * x + b
        @test gradient(ch, x) ≈ a
        @test abs(hessian(ch, x)) < 3e-14
    end
    @test ch(last(xs)) ≈ a * last(xs) + b
    @test gradient(ch, last(xs)) ≈ a
    @test abs(hessian(ch, last(xs))) < 3e-14

    # Now the interpolation condition:
    xs = zeros(50)
    ys = zeros(50)
    dydxs = zeros(50)
    xs[1] = 0.0
    ys[1] = randn()
    dydxs[1] = randn()
    for i = 2:50
        xs[i] = xs[i-1] + rand() + 0.1
        ys[i] = randn()
        dydxs[i] = randn()
    end

    ch = CubicHermite(xs, ys, dydxs)

    for i = 1:50
        @test ch(xs[i]) ≈ ys[i]
        @test gradient(ch, xs[i]) ≈ dydxs[i]
    end

    # Now quadratics:
    a = 1.0 / 8
    b = -3
    c = -2
    for i = 1:50
        ys[i] = a * xs[i] * xs[i] + b * xs[i] + c
        dydxs[i] = 2 * a * xs[i] + b
    end
    ch = CubicHermite(xs, ys, dydxs)
    for i = 1:200
        x = rand() * last(xs)
        @test ch(x) ≈ a * x * x + b * x + c
        @test gradient(ch, x) ≈ 2 * a * x + b
        @test hessian(ch, x) ≈ 2 * a
    end
    x = last(xs)
    @test ch(x) ≈ a * x * x + b * x + c
    @test gradient(ch, x) ≈ 2 * x * a + b
    @test hessian(ch, x) ≈ 2 * a
    # Cannot extrapolate:
    @test_throws DomainError ch(x + 0.1)
    @test_throws DomainError ch(xs[1] - 0.1)
end

@testset "Hermite" begin
    cubic_hermite_test()
end


text:
function get_a(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_co2ab(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getCo2AB", RealFieldElement, ())
end

function get_field_auxiliary_elements(obj::FieldForceModelContext)
    return jcall(obj, "getFieldAuxiliaryElements", FieldAuxiliaryElements, ())
end

function get_mean_motion(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getMeanMotion", RealFieldElement, ())
end

function get_mu(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getMu", RealFieldElement, ())
end

function get_oo_bpo(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getOoBpo", RealFieldElement, ())
end

function get_oo_mu(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getOoMU", RealFieldElement, ())
end

function get_ooa(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getOOA", RealFieldElement, ())
end

function get_ooab(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getOOAB", RealFieldElement, ())
end

function get_ton2a(obj::FieldAbstractGaussianContributionContext)
    return jcall(obj, "getTon2a", RealFieldElement, ())
end


text:
using 2022_01_comext
using Test

@testset "2022_01_comext.jl" begin
    # Write your tests here.
end


text:
# # LOCAL DEFINITION OF SCITYPE

# This is to avoid making ScientificTypes a dependency.

function scitype(c::CategoricalValue)
    nc = length(levels(c.pool))
    return ifelse(c.pool.ordered, OrderedFactor{nc}, Multiclass{nc})
end

# # CLASSES

"""
    classes(x)

Return, as a `CategoricalVector`, all the categorical elements with
the same pool as `CategoricalValue` `x` (including `x`), with an
ordering consistent with the pool. Note that `x in classes(x)` is
always true.

Not to be confused with `levels(x.pool)`. See the example below.

Also, overloaded for `x` a `CategoricalArray`, `CategoricalPool`, and
for views of `CategoricalArray`.

**Private method.*

    julia>  v = categorical([:c, :b, :c, :a])
    4-element CategoricalArrays.CategoricalArray{Symbol,1,UInt32}:
     :c
     :b
     :c
     :a

    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c

    julia> x = v[4]
    CategoricalArrays.CategoricalValue{Symbol,UInt32} :a

    julia> classes(x)
    3-element CategoricalArrays.CategoricalArray{Symbol,1,UInt32}:
     :a
     :b
     :c

    julia> levels(x.pool)
    3-element Array{Symbol,1}:
     :a
     :b
     :c

"""
classes(p::CategoricalPool) = [p[i] for i in 1:length(p)]
classes(x::CategoricalValue) = classes(CategoricalArrays.pool(x))
classes(v::CategoricalArray) = classes(CategoricalArrays.pool(v))
classes(v::SubArray{<:Any, <:Any, <:CategoricalArray}) = classes(parent(v))

# # CATEGORICAL VALUES TO INTEGERS

"""
   int(x)

The positional integer of the `CategoricalValue` `x`, in the ordering
defined by the pool of `x`. The type of `int(x)` is the reference type
of `x` (which differentiates this method from
`CategoricalArrays.levelcode`).

    int(X::CategoricalArray)
    int(W::AbstractArray{<:CategoricalValue})

Broadcasted versions of `int`.

    julia> v = categorical([:c, :b, :c, :a])
    julia> levels(v)
    3-element Array{Symbol,1}:
     :a
     :b
     :c
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001

See  [`decoder`](@ref) on how to invert the `int` operation.
"""
int(x) = throw(
    DomainError(x, "Can only convert categorical elements to integers. "))

int(x::Missing)       = missing
int(x::AbstractArray) = int.(x)

# first line is no good because it promotes type to larger integer type:
# int(x::CategoricalValue) = CategoricalArrays.levelcode(x)
int(x::CategoricalValue) = CategoricalArrays.refcode(x)

# # INTEGERS BACK TO CATEGORICAL VALUES

struct CategoricalDecoder{V,R}
    classes::CategoricalVector{V, R, V, CategoricalValue{V,R}, Union{}}
end

"""
    d = decoder(x)

A callable object for decoding the integer representation of a
`CategoricalValue` sharing the same pool as the `CategoricalValue`
`x`. Specifically, one has `d(int(y)) == y` for all `y` in the same
pool as `x`. One can also call `d` on integer arrays, in which case
`d` is broadcast over all elements.

    julia> v = categorical([:c, :b, :c, :a])
    julia> int(v)
    4-element Array{UInt32,1}:
     0x00000003
     0x00000002
     0x00000003
     0x00000001
    julia> d = decoder(v[3])
    julia> d(int(v)) == v
    true

*Warning:* There is no guarantee that `int(d(u)) == u` will always holds.

See also: [`int`](@ref).

"""
decoder(x) = CategoricalDecoder(classes(x))

(d::CategoricalDecoder{V,R})(i::Integer) where {V,R} =
    CategoricalValue{V,R}(d.classes[i])
(d::CategoricalDecoder)(a::AbstractArray{<:Integer}) = d.(a)

## TRANSFORMING BETWEEN CATEGORICAL ELEMENTS AND RAW VALUES

err_missing_class(c) =  DomainError("Value `$c` not in pool")

function _transform(pool, x)
    ismissing(x) && return missing
    x in levels(pool) || throw(err_missing_class(x))
    return pool[get(pool, x)]
end

_transform(pool, X::AbstractArray) = broadcast(x -> _transform(pool, x), X)

"""
    transform(e::Union{CategoricalElement,CategoricalArray,CategoricalPool},  X)

Transform the specified object `X` into a categorical version, using
the pool contained in `e`. Here `X` is a raw value (an element of
`levels(e)`) or an `AbstractArray` of such values.

```julia
v = categorical(["x", "y", "y", "x", "x"])
julia> transform(v, "x")
CategoricalValue{String,UInt32} "x"

julia> transform(v[1], ["x" "x"; missing "y"])
2×2 CategoricalArray{Union{Missing, Symbol},2,UInt32}:
 "x"       "x"
 missing   "y"

**Private method.**

"""
transform(e::Union{CategoricalArray, CategoricalValue},
                            arg) = _transform(CategoricalArrays.pool(e), arg)
transform(e::CategoricalPool, arg) =
    _transform(e, arg)