text:
###############################################################################
#
#  Ray Class Group: Ctx
#
###############################################################################

mutable struct ctx_rayclassgrp
  order::NfOrd
  class_group_map::MapClassGrp #The class group mod n map
  n::Int #the n for n_quo
  diffC::fmpz #exponent of the full class group, divided by n
  units::Vector{FacElem{nf_elem, AnticNumberField}}
  princ_gens::Vector{FacElem{nf_elem, AnticNumberField}}

  computed::Vector{Tuple{Dict{NfOrdIdl, Int}, Bool, MapRayClassGrp}}

  function ctx_rayclassgrp()
    z = new()
    return z
  end
end

function rayclassgrp_ctx(O::NfOrd, expo::Int)

  C1, mC1 = class_group(O)
  valclass = ppio(exponent(C1), fmpz(expo))[1]
  C, mC = Hecke.n_part_class_group(mC1, expo)
  ctx = ctx_rayclassgrp()
  ctx.order = O
  ctx.class_group_map = mC
  ctx.n = expo
  ctx.diffC = Int(divexact(exponent(C1), exponent(C)))
  return ctx

end

function assure_elements_to_be_eval(ctx::ctx_rayclassgrp)
  if isdefined(ctx, :units)
    return nothing
  end
  OK = ctx.order
  U, mU = unit_group_fac_elem(OK)
  mC = ctx.class_group_map
  _assure_princ_gen(mC)
  n = ctx.n
  units = FacElem{nf_elem, AnticNumberField}[_preproc(mU(U[i]), fmpz(n)) for i = 1:ngens(U)]
  princ_gens = FacElem{nf_elem, AnticNumberField}[_preproc(mC.princ_gens[i][2], fmpz(n)) for i = 1:length(mC.princ_gens)]
  ctx.units = units
  ctx.princ_gens = princ_gens
  return nothing
end

###############################################################################
#
#  Ray Class Group: Fields function
#
###############################################################################

function empty_ray_class(m::NfOrdIdl)
  O = order(parent(m))
  X = DiagonalGroup(Int[])

  local exp
  let O = O
    function exp(a::GrpAbFinGenElem)
      return FacElem(Dict(ideal(O,1) => fmpz(1)))
    end
  end

  local disclog
  let X = X
    function disclog(J::Union{NfOrdIdl, FacElem{NfOrdIdl}})
      return X(Int[])
    end
  end

  mp = Hecke.MapRayClassGrp()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , exp, disclog)
  mp.defining_modulus = (m, InfPlc[])
  return X,mp

end

function ray_class_group_quo(m::NfOrdIdl, y1::Dict{NfOrdIdl,Int}, y2::Dict{NfOrdIdl,Int}, inf_plc::Vector{InfPlc}, ctx::ctx_rayclassgrp; check::Bool = true, GRH::Bool = true)

  mC = ctx.class_group_map
  C = domain(mC)
  diffC = ctx.diffC

  if isempty(y1) && isempty(y2) && isempty(inf_plc)
    local exp_c
    let mC = mC
      function exp_c(a::GrpAbFinGenElem)
        return FacElem(Dict(mC(a) => 1))
      end
    end
    return class_as_ray_class(quo(C, ctx.n, false)[1], mC, exp_c, m)
  end

  O = ctx.order
  K = nf(O)
  n_quo = ctx.n

  lp = merge(max, y1, y2)

  powers = Vector{Tuple{NfOrdIdl, NfOrdIdl}}()
  quo_rings = Tuple{NfOrdQuoRing, Hecke.AbsOrdQuoMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[]
  groups_and_maps = Tuple{GrpAbFinGen, Hecke.GrpAbFinGenToAbsOrdQuoRingMultMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[]
  for (PP, ee) in lp
    if isone(ee)
      dtame = Dict{NfOrdIdl, Int}(PP => 1)
      dwild = Dict{NfOrdIdl, Int}()
    else
      dwild = Dict{NfOrdIdl, Int}(PP => ee)
      if haskey(y1, PP)
        dtame = Dict{NfOrdIdl, Int}(PP => 1)
      else
        dtame = Dict{NfOrdIdl, Int}()
      end
    end
    QQ = PP^ee

    push!(powers, (PP, QQ))
    QQQ, mQQQ = quo(O, QQ)
    if ee == 1
      QQQ.factor = dtame
    else
      QQQ.factor = dwild
    end
    push!(quo_rings, (QQQ, mQQQ))
    push!(groups_and_maps, _mult_grp_mod_n(quo_rings[end][1], dtame, dwild, n_quo))
  end

  if isempty(groups_and_maps)
    nG = 0
    expon = fmpz(1)
  else
    nG = sum(ngens(x[1]) for x in groups_and_maps)
    expon = lcm([exponent(x[1]) for x in groups_and_maps])
  end

  if iseven(n_quo)
    p = InfPlc[ x for x in inf_plc if isreal(x) ]
  else
    p = InfPlc[]
  end
  H, lH = log_infinite_primes(O, p)
  expon = lcm(expon, exponent(H))

  if exponent(C)*expon < n_quo && check
    return empty_ray_class(m)
  end

  assure_elements_to_be_eval(ctx)
  exp_class, Kel = find_coprime_representatives(mC, m, lp)

  nU = length(ctx.units)
  # We construct the relation matrix and evaluate units and relations with the class group in the quotient by m
  # Then we compute the discrete logarithms

  R = zero_matrix(FlintZZ, 2*(ngens(C)+nG+ngens(H))+nU, ngens(C)+ngens(H)+nG)
  for i = 1:ncols(R)
    R[i+ngens(C)+nG+ngens(H)+nU, i] = n_quo
  end
  for i=1:ngens(C)
    R[i, i] = C.snf[i]
  end
  ind = 1
  for s = 1:length(quo_rings)
    G = groups_and_maps[s][1]
    @assert issnf(G)
    for i = 1:ngens(G)
      R[i+ngens(C)+ind-1, i+ngens(C)+ind-1] = G.snf[i]
    end
    ind += ngens(G)
  end
  for i = 1:ngens(H)
    R[ngens(C)+nG+i, ngens(C)+nG+i] = 2
  end

  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  tobeeval = Vector{FacElem{nf_elem, AnticNumberField}}(undef, nU+ngens(C))
  for i = 1:nU
    tobeeval[i] = ctx.units[i]
  end
  for i = 1:ngens(C)
    tobeeval[i+nU] = _preproc(ctx.princ_gens[i]*FacElem(Dict(Kel[i] => C.snf[i])), expon)
  end

  ind = 1
  for i = 1:length(groups_and_maps)
    exp_q = gcd(expon, norm(powers[i][2])- divexact(norm(powers[i][2]), norm(powers[i][1])))
    evals = fac_elems_eval(powers[i][1], powers[i][2], tobeeval, exp_q)
    Qr = quo_rings[i][1]
    mG = groups_and_maps[i][2]
    for j = 1:nU
      a = (mG\Qr(evals[j])).coeff
      for s = 1:ncols(a)
        R[j+nG+ngens(H)+ngens(C), ngens(C)+s+ind-1] = a[1, s]
      end
    end

    for j = 1:ngens(C)
      a = (mG\Qr(evals[j+nU])).coeff
      for s = 1:ncols(a)
        R[j, ngens(C)+ind+s-1] = -a[1, s]
      end
    end
    ind += ngens(groups_and_maps[i][1])
  end

  if !isempty(p)
    for j = 1:nU
      aa = lH(tobeeval[j])::GrpAbFinGenElem
      for s = 1:ngens(H)
        R[j+nG+ngens(C)+ngens(H), ngens(C)+ind-1+s] = aa[s]
      end
    end
    for j = 1:ngens(C)
      aa = lH(tobeeval[j+nU])::GrpAbFinGenElem
      for s = 1:ngens(H)
        R[j, ngens(C)+ind-1+s] = -aa[s]
      end
    end
  end

  X = AbelianGroup(R)

  invd = invmod(fmpz(diffC), expon)
  local disclog
  let X = X, mC = mC, invd = invd, C = C, exp_class = exp_class, powers = powers, groups_and_maps = groups_and_maps, quo_rings = quo_rings, lH = lH, diffC = diffC, n_quo = n_quo, m = m, p = p, expon = expon

    # Discrete logarithm
    function disclog(J::FacElem{NfOrdIdl, NfOrdIdlSet})
      @vprint :RayFacElem 1 "Disc log of element $J \n"
      a1 = id(X)
      for (f, k) in J
        a1 += k*disclog(f)
      end
      return a1
    end

    function disclog(J::NfOrdIdl)
      @hassert :RayFacElem 1 iscoprime(J, m)
      if isone(J)
        @vprint :RayFacElem 1 "J is one \n"
        return id(X)
      end
      coeffs = zero_matrix(FlintZZ, 1, ngens(X))
      if J.is_principal == 1 && isdefined(J, :princ_gen)
        z = FacElem(Dict(J.princ_gen.elem_in_nf => diffC))
      else
        L = mC\J
        for i = 1:ngens(C)
          coeffs[1, i] = L[i]
        end
        @vprint :RayFacElem 1 "Disc log of element J in the Class Group: $(L.coeff) \n"
        eL = exp_class(L)
        inv!(eL)
        add_to_key!(eL.fac, J, 1)
        pow!(eL, diffC)
        @vprint :RayFacElem 1 "This ideal is principal: $eL \n"
        z = principal_gen_fac_elem(eL)
      end
      ii = 1
      z1 = _preproc(z, expon)
      for i = 1:length(powers)
        P, Q = powers[i]
        exponq = gcd(expon, norm(Q)-divexact(norm(Q), norm(P)))
        el = fac_elems_eval(P, Q, FacElem{nf_elem, AnticNumberField}[z1], exponq)
        y = (invd*(groups_and_maps[i][2]\quo_rings[i][1](el[1]))).coeff
        for s = 1:ncols(y)
          coeffs[1, ii-1+ngens(C)+s] = y[1, s]
        end
        ii += ngens(groups_and_maps[i][1])
      end
      if !isempty(p)
        b = lH(z).coeff
        for s = 1:ncols(b)
          coeffs[1, ii-1+s+ngens(C)] = b[1, s]
        end
      end
      return X(coeffs)
    end
  end

  Dgens = Tuple{NfOrdElem, GrpAbFinGenElem}[]
  ind = 1
  #We need generators of the full multiplicative group
  #In particular, we need the idempotents...
  if !isone(length(powers))
    for i = 1:length(powers)
      P, QQ = powers[i]
      mG = groups_and_maps[i][2]
      J = ideal(O, 1)
      minJ = fmpz(1)
      mins = fmpz(1)
      for j = 1:length(powers)
        if minimum(powers[j][1], copy = false) != minimum(P, copy = false)
          mins = lcm(mins, minimum(powers[j][2], copy = false))
          continue
        end
        if j != i
          J *= powers[j][2]
          minJ = lcm(minJ, minimum(powers[j][2], copy = false))
        end
      end
      J.minimum = minJ
      i1, i2 = idempotents(QQ, J)
      if !isone(mins)
        d, u1, v1 = gcdx(minimum(QQ, copy = false), mins)
        i1 = i1*(u1*minimum(QQ, copy = false) + v1*mins) + u1*minimum(QQ, copy = false) *i2
        i2 = v1*mins*i2
      end
      @hassert :RayFacElem 1 isone(i1+i2)
      @hassert :RayFacElem 1 i1 in QQ
      @hassert :RayFacElem 1 i2 in J
      @hassert :RayFacElem 1 i2 in ideal(O, mins)
      gens = mG.generators
      if isempty(p)
        if haskey(mG.tame, P)
          gens_tame = mG.tame[P].generators
          for s = 1:length(gens_tame)
            gens_tame[s] = gens_tame[s]*i2 + i1
            @hassert :RayFacElem 1 isone(mod(gens_tame[s], J*ideal(O, mins)))
          end
          mG.tame[P].generators = gens_tame
        end
        if haskey(mG.wild, P)
          gens_wild = mG.wild[P].generators
          for s = 1:length(gens_wild)
            gens_wild[s] = gens_wild[s]*i2 + i1
            @hassert :RayFacElem 1 isone(mod(gens_wild[s], J*ideal(O, mins)))
          end
          mG.wild[P].generators = gens_wild
        end
        for s = 1:length(gens)
          el_to_push = gens[s].elem*i2+i1
          @hassert :RayFacElem 1 isone(mod(el_to_push, J*ideal(O, mins)))
          push!(Dgens, (el_to_push, X[ngens(C)+ind-1+s]))
        end
      else
        if haskey(mG.tame, P)
          gens_tame = mG.tame[P].generators
          for s = 1:length(gens_tame)
            gens_tame[s] = make_positive(gens_tame[s]*i2 + i1, minimum(m, copy = false))
          end
          mG.tame[P].generators = gens_tame
        end
        if haskey(mG.wild, P)
          gens_wild = mG.wild[P].generators
          for s = 1:length(gens_wild)
            gens_wild[s] = make_positive(gens_wild[s]*i2 + i1, minimum(m, copy = false))
          end
          mG.wild[P].generators = gens_wild
        end
        for s = 1:length(gens)
          elgen = make_positive(gens[s].elem*i2 + i1, minimum(m, copy = false))
          push!(Dgens, (elgen, X[ngens(C)+ind-1+s]))
        end
      end
      ind += length(gens)
    end
  else
    mG = groups_and_maps[1][2]
    gens = mG.generators
    if !isempty(p)
      for s = 1:length(gens)
        elgen = make_positive(gens[s].elem, minimum(m, copy = false))
        push!(Dgens, (elgen, X[ngens(C)+s]))
      end
    else
      for s = 1:length(gens)
        elgen = gens[s].elem
        push!(Dgens, (elgen, X[ngens(C)+s]))
      end
    end
  end

  ind = 1
  for i = 1:length(powers)
    mG = groups_and_maps[i][2]
    for (prim, mprim) in mG.tame
      dis = zero_matrix(FlintZZ, 1, ngens(X))
      to_be_c = mprim.disc_log.coeff
      for i = 1:length(to_be_c)
        dis[1, ind-1+i+ngens(C)] = to_be_c[1, i]
      end
      mprim.disc_log = X(dis)
    end
    ind += ngens(domain(mG))
  end

  disc_log_inf = Dict{InfPlc, GrpAbFinGenElem}()
  for i = 1:length(p)
    eldi = zeros(FlintZZ, ngens(X))
    eldi[ngens(X) - length(inf_plc) + i] = 1
    disc_log_inf[p[i]] = X(eldi)
  end

  mp = MapRayClassGrp()
  mp.header = MapHeader(X, FacElemMon(parent(m)))
  mp.header.preimage = disclog
  mp.fact_mod = lp
  mp.defining_modulus = (m, p)
  mp.powers = powers
  mp.groups_and_maps = groups_and_maps
  mp.disc_log_inf_plc = disc_log_inf
  mp.gens_mult_grp_disc_log = Dgens
  mp.clgrpmap = mC
  return X, mp

end

function log_infinite_primes(O::NfOrd, p::Array{InfPlc,1})
  if isempty(p)
    S = DiagonalGroup(Int[])

    local log1
    let S = S
      function log1(B::T) where T <: Union{nf_elem ,FacElem{nf_elem, AnticNumberField}}
        return id(S)
      end
    end
    return S, log1
  end

  S = DiagonalGroup(Int[2 for i=1:length(p)])
  local log
  let S = S, p = p
    function log(B::T) where T <: Union{nf_elem ,FacElem{nf_elem, AnticNumberField}}
      emb = signs(B, p)
      ar = zero_matrix(FlintZZ, 1, length(p))
      for i = 1:length(p)
        if emb[p[i]] == -1
          ar[1, i] = 1
        end
      end
      return S(ar)
    end
  end
  return S, log
end

function ray_class_group_quo(O::NfOrd, m::Int, wprimes::Dict{NfOrdIdl,Int}, inf_plc::Array{InfPlc,1}, ctx::ctx_rayclassgrp; GRH::Bool = true)

  d1 = Dict{NfOrdIdl, Int}()
  lp = factor(m)
  I = ideal(O, 1)
  minI = fmpz(1)
  for q in keys(lp.fac)
    lq = prime_decomposition(O, q)
    for (P, e) in lq
      I *= P
      d1[P] = 1
    end
    minI = q*minI
  end
  if !isempty(wprimes)
    for (p, v) in wprimes
      qq = p^v
      I *= qq
      minI = lcm(minI, minimum(qq, copy = false))
    end
  end
  I.minimum = minI
  return ray_class_group_quo(I, d1, wprimes, inf_plc, ctx, GRH = GRH)

end

function ray_class_group_quo(O::NfOrd, y::Dict{NfOrdIdl, Int}, inf_plc::Array{InfPlc, 1}, ctx::ctx_rayclassgrp; GRH::Bool = true, check::Bool = true)

  y1=Dict{NfOrdIdl,Int}()
  y2=Dict{NfOrdIdl,Int}()
  n = ctx.n
  for (q,e) in y
    if gcd(norm(q)-1, n) != 1
      y1[q] = Int(1)
      if gcd(norm(q), n) != 1 && e >= 2
        y2[q] = Int(e)
      end
    elseif gcd(norm(q), n) != 1 && e >= 2
      y2[q] = Int(e)
    end
  end
  I=ideal(O,1)
  for (q,vq) in y1
    I*=q
  end
  for (q,vq) in y2
    I*=q^vq
  end
  return ray_class_group_quo(I, y1, y2, inf_plc, ctx, GRH = GRH, check = check)

end

###############################################################################
#
#  Maximal abelian subfield for fields function
#
###############################################################################

function check_abelian_extensions(class_fields::Vector{Tuple{Hecke.ClassField{Hecke.MapRayClassGrp,GrpAbFinGenMap}, Vector{GrpAbFinGenMap}}}, autos::Array{NfToNfMor, 1}, emb_sub::NfToNfMor)
  @vprint :MaxAbExt 3 "Starting checking abelian extension\n"
  K = base_field(class_fields[1][1])
  d = degree(K)
  G = domain(class_fields[1][2][1])
  expo = G.snf[end]
  com, uncom = ppio(Int(expo), d)
  if com == 1
    return trues(length(class_fields))
  end
  #I extract the generators that restricted to the domain of emb_sub are the identity.
  #Notice that I can do this only because I know the way I am constructing the generators of the group.
  expG_arr = Int[]
  act_indices = Int[]
  p = 11
  d1 = discriminant(domain(emb_sub))
  d2 = discriminant(K)
  while iszero(mod(d1, p)) || iszero(mod(d2, p))
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  mp_pol = Rx(emb_sub.prim_img)
  for i = 1:length(autos)
    pol = Rx(autos[i].prim_img)
    if mp_pol ==  Hecke.compose_mod(mp_pol, pol, fmod)
      push!(act_indices, i)
      #I compute the order of the automorphisms. I need the exponent of the relative extension!
      push!(expG_arr, _order(autos[i]))
    end
  end
  expG = lcm(expG_arr)
  expG1 = ppio(expG, com)[1]
  @vprint :MaxAbExt 3 "Context for ray class groups\n"

  OK = maximal_order(K)
  rcg_ctx = Hecke.rayclassgrp_ctx(OK, com*expG1)

  cfields = trues(length(class_fields))
  for i = 1:length(class_fields)
    @vprint :MaxAbExt 3 "Class Field $i\n"
    C, res_act = class_fields[i]
    res_act_new = Vector{GrpAbFinGenMap}(undef, length(act_indices))
    for i = 1:length(act_indices)
      res_act_new[i] = res_act[act_indices[i]]
    end
    cfields[i] = check_abelian_extension(C, res_act_new, emb_sub, rcg_ctx, expG)
  end
  return cfields

end

function check_abelian_extension(C::Hecke.ClassField, res_act::Vector{GrpAbFinGenMap}, emb_sub::NfToNfMor, rcg_ctx::Hecke.ctx_rayclassgrp, exponent_extension::Int)

  #I consider the action on every P-sylow and see if it is trivial
  G = codomain(C.quotientmap)
  expG = Int(exponent(G))
  fac = factor(rcg_ctx.n)
  prime_to_test = Int[]
  new_prime = false
  for (P, v) in fac
    # I check that the action of the P-Sylow has no fixed points.
    PS, mPS = psylow_subgroup(G, P, false)
    s, ms = snf(PS)
    act_sub = induce_action_on_subgroup(ms*mPS, res_act)
    if !isfixed_point_free(act_sub)
      push!(prime_to_test, P)
    end
  end
  if isempty(prime_to_test)
    return true
  end
  n = prod(prime_to_test)
  n1, m = ppio(Int(G.snf[end]), n)
  if m != 1
    # Now, I compute the maximal abelian subextension of the n-part of C
    Q, mQ = quo(G, n1, false)
    C1 = ray_class_field(C.rayclassgroupmap, Hecke.GrpAbFinGenMap(Hecke.compose(C.quotientmap, mQ)))
    #@vtime :MaxAbExt 1
    fl = _maximal_abelian_subfield(C1, emb_sub, rcg_ctx, exponent_extension)
  else
    #@vtime :MaxAbExt 1
    fl = _maximal_abelian_subfield(C, emb_sub, rcg_ctx, exponent_extension)
  end
  return fl

end

function _bound_exp_conductor_wild(O::NfOrd, n::Int, q::Int, bound::fmpz)
  d = degree(O)
  lp = prime_decomposition_type(O, q)
  f_times_r = divexact(d, lp[1][2])
  sq = fmpz(q)^f_times_r
  nbound = n+n*lp[1][2]*valuation(n,q)-div(fmpz(n), q^(valuation(n,q)))
  obound = flog(bound, sq)
  bound_max_ap = min(nbound, obound)  #bound on ap
  return div(q*bound_max_ap, n*(q-1)) #bound on the exponent in the conductor
end

function minimumd(D::Dict{NfOrdIdl, Int}, deg_ext::Int)
  primes_done = Int[]
  res = 1
  for (P, e) in D
    p = Int(minimum(P))
    if p in primes_done
      continue
    else
      push!(primes_done, p)
    end
    ram_index = P.splitting_type[1]
    s, t = divrem(e, ram_index)
    if iszero(t)
      d = min(s, valuation(deg_ext, p)+2)
      res *= p^d
    else
      d = min(s+1, valuation(deg_ext, p)+2)
      res *= p^d
    end
  end
  return res
end

function _maximal_abelian_subfield(A::Hecke.ClassField, mp::Hecke.NfToNfMor, ctx::Hecke.ctx_rayclassgrp, exp_extension::Int)
  K = base_field(A)
  k = domain(mp)
  ZK = maximal_order(K)
  zk = maximal_order(k)
  expected_order = div(degree(K), degree(k))
  if gcd(expected_order, degree(A)) == 1
    return false
  end
  # disc(ZK/Q) = N(disc(ZK/zk)) * disc(zk)^deg
  # we need the disc ZK/k, well a conductor.
  d = abs(div(discriminant(ZK), discriminant(zk)^expected_order))
  mR1 = A.rayclassgroupmap
  expo = Int(exponent(codomain(A.quotientmap)))
  deg = expected_order
  #First, a suitable modulus for A over k
  #I take the discriminant K/k times the norm of the conductor A/K
  fm0 = Dict{NfOrdIdl, Int}()
  for (P, e) in mR1.fact_mod
    p = intersect_prime(mp, P)
    if !haskey(fm0, p)
      if !iscoprime(minimum(P, copy = false), expo)
        if e > 1
          fm0[p] = e
        end
      else
        fm0[p] = 1
      end
    end
  end
  lp = factor(ideal(zk, d))
  for (P, e) in lp
    if haskey(fm0, P)
      if !iscoprime(minimum(P, copy = false), deg*expo)
        fm0[P] = max(e, fm0[P])
      end
    else
      if !iscoprime(minimum(P, copy = false), deg*expo)
        fm0[P] = e
      else
        fm0[P] = 1
      end
    end
  end
  #Now, I extend this modulus to K
  fM0 = Dict{NfOrdIdl, Int}()
  for (p, v) in fm0
    lp = prime_decomposition(mp, p)
    if iscoprime(minimum(p, copy = false), expo*deg)
      for (P, e) in lp
        fM0[P] = 1
      end
    else
      for (P, e) in lp
        fM0[P] = e*v
      end
    end
  end
  ind = 0
  if isdefined(ctx, :computed)
    flinf = isempty(defining_modulus(mR1)[2])
    for i = 1:length(ctx.computed)
      idmr, ifmr, mRRR = ctx.computed[i]
      if flinf != ifmr
        continue
      end
      contained = true
      for (P, v) in fM0
        if !haskey(idmr, P) || idmr[P] < v
          contained = false
        end
      end
      if contained
        ind = i
        break
      end
    end
  end
  if iszero(ind)
    @vtime :MaxAbExt 3 R, mR = Hecke.ray_class_group_quo(ZK, fM0, defining_modulus(mR1)[2], ctx, check = false)
    if isdefined(ctx, :computed)
      push!(ctx.computed, (fM0, isempty(defining_modulus(mR1)[2]), mR))
    else
      ctx.computed = [(fM0, isempty(defining_modulus(mR1)[2]), mR)]
    end
  else
    mR = ctx.computed[ind][3]
    R = domain(mR)
  end
  if degree(zk) != 1
    if istotally_real(K) && isempty(defining_modulus(mR1)[2])
      inf_plc = InfPlc[]
    else
      inf_plc = real_places(k)
    end
    #@vtime :MaxAbExt 1
    r, mr = Hecke.ray_class_group(zk, fm0, inf_plc, n_quo = ctx.n)
  else
    rel_plc = true
    if istotally_real(K) && isempty(defining_modulus(mR1)[2])
      rel_plc = false
    end
    modulo = minimumd(fm0, expo * expected_order)
    #@vtime :MaxAbExt 1
    r, mr = Hecke.ray_class_groupQQ(zk, modulo, rel_plc, ctx.n)
  end
  @vtime :MaxAbExt 1 lP, gS = Hecke.find_gens(mR, coprime_to = minimum(defining_modulus(mR1)[1]))
  listn = NfOrdIdl[norm(mp, x) for x in lP]
  # Create the map between R and r by taking norms
  preimgs = Vector{GrpAbFinGenElem}(undef, length(listn))
  for i = 1:length(preimgs)
    preimgs[i] = mr\listn[i]
  end
  proj = hom(gS, preimgs)
  #compute the norm group of A in R
  prms = Vector{GrpAbFinGenElem}(undef, length(lP))
  for i = 1:length(lP)
    if haskey(mR1.prime_ideal_preimage_cache, lP[i])
      f = mR1.prime_ideal_preimage_cache[lP[i]]
    else
      f = mR1\lP[i]
      mR1.prime_ideal_preimage_cache[lP[i]] = f
    end
    prms[i] = A.quotientmap(f)
  end
  proj1 = hom(gS, prms)
  S, mS = kernel(proj1, false)
  mS1 = mS*proj
  G, mG = Hecke.cokernel(mS1, false)
  return (order(G) == expected_order)::Bool

end


text:
function Base.:-(a::GeoArray, b::GeoArray)
    size(a) == size(b) || throw(DimensionMismatch("The sizes of a and b do not match."))
    a.f == b.f || error("The affine information does not match")
    a.crs == b.crs || error("The crs information does not match")
    GeoArray(a.A .- b.A, a.f, a.crs)
end

function Base.:+(a::GeoArray, b::GeoArray)
    size(a) == size(b) || throw(DimensionMismatch("The sizes of a and b do not match."))
    a.f == b.f || error("The affine information does not match")
    a.crs == b.crs || error("The crs information does not match")
    GeoArray(a.A .+ b.A, a.f, a.crs)
end

function Base.:*(a::GeoArray, b::GeoArray)
    size(a) == size(b) || throw(DimensionMismatch("The sizes of a and b do not match."))
    a.f == b.f || error("The affine information does not match")
    a.crs == b.crs || error("The crs information does not match")
    GeoArray(a.A .* b.A, a.f, a.crs)
end

function Base.:/(a::GeoArray, b::GeoArray)
    size(a) == size(b) || throw(DimensionMismatch("The sizes of a and b do not match."))
    a.f == b.f || error("The affine information does not match")
    a.crs == b.crs || error("The crs information does not match")
    GeoArray(a.A ./ b.A, a.f, a.crs)
end

ext: jl


text:
"""
Bridges `CP.Reification{CP.DifferentFrom}` to indicator constraints, both with equality
and inequalities (CP.DifferentFrom).
"""
struct ReificationDifferentFrom2IndicatorBridge{T <: Real} <: MOIBC.AbstractBridge
    indic_true::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}}
    indic_false::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{ReificationDifferentFrom2IndicatorBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Reification{CP.DifferentFrom{T}},
) where {T}
    return MOIBC.bridge_constraint(
        ReificationDifferentFrom2IndicatorBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{ReificationDifferentFrom2IndicatorBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Reification{CP.DifferentFrom{T}},
) where {T <: Real}
    # Only change with respect to CP.Reification{MOI.EqualTo}: change the
    # sign of the first variable.
    f_scalars = MOIU.scalarize(f)
    f = MOIU.vectorize(
        [
            one(T) - f_scalars[1],
            f_scalars[2]
        ]
    )

    indic_true = MOI.add_constraint(
        model,
        f,
        MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(s.set.value))
    )
    indic_false = MOI.add_constraint(
        model,
        f,
        MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(CP.DifferentFrom(s.set.value))
        # TODO: helper to build CP.\neq from MOI.EqTo, CP.Strictly from inequalities, like `!()`?
    )

    return ReificationDifferentFrom2IndicatorBridge{T}(indic_true, indic_false)
end

function MOI.supports_constraint(
    ::Type{ReificationDifferentFrom2IndicatorBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Reification{CP.DifferentFrom{T}}},
) where {T <: Real}
    return true
end

function MOIB.added_constrained_variable_types(::Type{ReificationDifferentFrom2IndicatorBridge{T}}) where {T <: Real}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{ReificationDifferentFrom2IndicatorBridge{T}}) where {T <: Real}
    return [
        (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}}),
        (MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}}),
    ]
end

function MOI.get(
    ::ReificationDifferentFrom2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}},
    },
) where {T <: Real}
    return 1
end

function MOI.get(
    ::ReificationDifferentFrom2IndicatorBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}},
    },
) where {T <: Real}
    return 1
end

function MOI.get(
    b::ReificationDifferentFrom2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ONE, MOI.EqualTo{T}},
    },
) where {T <: Real}
    return [b.indic_true]
end

function MOI.get(
    b::ReificationDifferentFrom2IndicatorBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, MOI.Indicator{MOI.ACTIVATE_ON_ZERO, CP.DifferentFrom{T}},
    },
) where {T <: Real}
    return [b.indic_false]
end


text:
"""
    perpendicular_vector(vec)

Compute a vector perpendicular to `vec` by switching the two elements with
largest absolute value, flipping the sign of the second largest, and setting the
remaining element to zero.
"""
function perpendicular_vector(vec::SVector{3})
    T = eltype(vec)

    # find indices of the two elements of vec with the largest absolute values:
    absvec = abs.(vec)
    ind1 = argmax(absvec) # index of largest element
    tmin = typemin(T)
    @inbounds absvec2 = @SVector [ifelse(i == ind1, tmin, absvec[i]) for i = 1 : 3] # set largest element to typemin(T)
    ind2 = argmax(absvec2) # index of second-largest element

    # perp[ind1] = -vec[ind2], perp[ind2] = vec[ind1], set remaining element to zero:
    @inbounds perpind1 = -vec[ind2]
    @inbounds perpind2 = vec[ind1]
    tzero = zero(T)
    perp = @SVector [ifelse(i == ind1, perpind1, ifelse(i == ind2, perpind2, tzero)) for i = 1 : 3]
end


text:
using Random, Distributions;

"""
# randomwalk_normal_bounded(n::Array{<:Number}; σ::Number, nmin::Number, nmax::Number, steps::Number=1)

Execute an arbitrary number of random walk steps on the values of the n array,
consisting in the addition of random values from the normal distribution defined by σ,
and bounding the values between a minimum and a maximum.

## Examples
```randomwalk_normal_bounded!([.1,.2,.3], σ = .1, nmin = 0, nmax = 1, steps = 1);
julia> randomwalk_normal_bounded!([.1,.2,.3], σ = .1, nmin = 0, nmax = 1, steps = 1)
3-element Array{Float64,1}:
0.0480345320809181
0.07878951041425464
0.3094678083144717
```
"""
function randomwalk_normal_bounded!(n::Array{<:Number}; σ::Number, nmin::Number, nmax::Number, steps::Number=1)
    for step = 1 : steps
        n = broadcast(
          min,
          broadcast(
            max,
            n + rand(Normal(0.0, σ), length(n)),
            nmin),
          nmax);
    end;
    return n;
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
using GLVisualize, GeometryTypes, Plots, GLWindow
using GLAbstraction
if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end
description = """
Combining GLVisualize with Plots.jl.
"""

plot_area, glv_area = y_partition(window.area, 50) # part parent area at 50%

# create one screen for Plots.jl
plot_screen = Screen(
    window, name=:plot, area=plot_area
)
# make it the current screen for plotting
GLVisualize.add_screen(plot_screen)
# tell Plots.jl the size of the window:
#
glvisualize(size=widths(plot_screen))

# one for GLVisualize
glv_screen = Screen(
    window, name=:glvisualize, area=glv_area
)

N = 10_000
_view(visualize(
    (Circle(Point3f0(0), 0.01f0), rand(Point3f0, N).*4f0),
    color = rand(RGBA{Float32}, N)
), glv_screen, camera=:perspective)

# plot one of the examples. This could be any Plots.plot command!
group = rand(map(i->"group $i",1:4),100)
plot(rand(100), layout=@layout([a b;c]), group=group, linetype=[:bar :scatter :steppre])
gui()

if !isdefined(:runtests)
    renderloop(window)
end


text:
# circle.jl

abstract type AbstractHyperSphere{N,T} <: EuclideanDomain{N,T} end

"""
The unit sphere (of radius 1) in `N` dimensions.
"""
struct UnitHyperSphere{N,T} <: AbstractHyperSphere{N,T} end

UnitHyperSphere{N}() where N = UnitHyperSphere{N,Float64}()

const UnitCircle{T} = UnitHyperSphere{2,T}
const UnitSphere{T} = UnitHyperSphere{3,T}

convert(::Type{Domain{SVector{N,T}}}, d::UnitHyperSphere{N}) where {N,T} =
    UnitHyperSphere{N,T}()

indomain(x, ::UnitHyperSphere) = norm(x) == 1

approx_indomain(x, ::UnitHyperSphere, tolerance) = 1-tolerance <= norm(x) <= 1+tolerance

boundary(::UnitHyperBall{N,T}) where {N,T} = UnitHyperSphere{N,T}()

isempty(::UnitHyperSphere) = false

"Create an ellipse curve with semi-axes lengths `a` and `b` respectively."
ellipse(a::Number, b::Number) = ellipse(promote(a,b)...)
ellipse(a::T, b::T) where {T <: Number} = scaling_map(a, b) * UnitCircle{T}()

"Create an ellipse-shaped domain with semi-axes lengths `a` and `b` respectively."
ellipse_shape(a::Number, b::Number) = ellipse_shape(promote(a,b)...)
ellipse_shape(a::T, b::T) where {T <: Number} = scaling_map(a, b) * UnitDisk{T}()

"""
The map `[cos(2πt), sin(2πt)]` from `[0,1)` to the unit circle in `ℝ^2`.
"""
struct UnitCircleMap{S,T} <: AbstractMap{S,T} end

parameterization(d::UnitCircle) = UnitCircleMap{subeltype(d),eltype(d)}()

domain(d::UnitCircleMap{S}) where S = HalfOpenRightInterval{S}(0, 1)

image(m::UnitCircleMap{S}) where S = UnitCircle{S}()

applymap(m::UnitCircleMap{S}, t) where S = SVector(cos(2*S(pi)*t), sin(2*S(pi)*t))

function gradient(m::UnitCircleMap{S}, t) where S
    a = 2*S(pi)
    SVector(-a*sin(a*t), a*cos(a*t))
end

"""
`AngleMap` is a left inverse of `UnitCircleMap`. A 2D vector `x` is projected onto
the intersection point with the unit circle of the line connecting `x` to the
origin. The angle of this point, scaled to the interval `[0,1)`, is the result.
"""
struct AngleMap{S,T} <: AbstractMap{S,T}
end

domain(d::AngleMap{S,T}) where {S,T} = FullSpace{S}()

range(m::AngleMap{S,T}) where {S,T} = HalfOpenRightInterval{T}(0, 1)

function applymap(m::AngleMap, x)
    twopi = 2*convert(codomaintype(m), pi)
    θ = atan(x[2],x[1])
    if θ < 0
        # atan2 returns an angle in (-π,π], convert to [0,2π) using periodicity.
        θ += twopi
    end
    # And divide by 2π to scale to [0,1)
    θ / twopi
end

left_inverse(m::UnitCircleMap{S,T}) where {S,T} = AngleMap{T,S}()

right_inverse(m::AngleMap{S,T}) where {S,T} = UnitCircleMap{T,S}()


text:
function inpcelltype(::Type{CT}) where {CT}
    if CT === Triangle
        return "CPS3"
    elseif CT === QuadraticTriangle
        return "CPS6"
    elseif CT === Tetrahedron
        return "C3D4"
    elseif CT === QuadraticTetrahedron
        return "C3D10"
    elseif CT === Quadrilateral
        return "CPS4"
    elseif CT === QuadraticQuadrilateral
        return "CPS8"
    elseif CT === Hexahedron
        return "C3D8"
    elseif CT === Ferrite.QuadraticHexahedron
        return "C3D20"
    else
        return ""
    end
end

function import_inp(filepath_with_ext)
    problem = extract_inp(filepath_with_ext)
    return inp_to_ferrite(problem)
end

function inp_to_ferrite(problem::InpContent)
    _celltype = problem.celltype
    if _celltype == "CPS3"
        # Linear triangle
        celltype = Triangle
        geom_order = 1
        refshape = RefTetrahedron
        dim = 2
    elseif _celltype == "CPS6"
        # Quadratic triangle
        celltype = QuadraticTriangle
        geom_order = 2
        refshape = RefTetrahedron
        dim = 2
    elseif _celltype == "C3D4"
        # Linear tetrahedron
        celltype = Tetrahedron
        geom_order = 1
        refshape = RefTetrahedron
        dim = 3
    elseif _celltype == "C3D10"
        # Quadratic tetrahedron
        celltype = QuadraticTetrahedron
        geom_order = 2
        refshape = RefTetrahedron
        dim = 3
    elseif _celltype == "CPS4"
        # Linear quadrilateral
        celltype = Quadrilateral
        geom_order = 1
        refshape = RefCube
        dim = 2
    elseif _celltype == "CPS8" || _celltype == "CPS8R"
        # Quadratic quadrilateral
        celltype = QuadraticQuadrilateral
        geom_order = 2
        refshape = RefCube
        dim = 2
    elseif _celltype == "C3D8" || _celltype == "C3D8R"
        # Linear hexahedron
        celltype = Hexahedron
        geom_order = 1
        refshape = RefCube
        dim = 3
    elseif _celltype == "C3D20" || _celltype == "C3D20R"
        # Quadratic hexahedron
        celltype = QuadraticHexahedron
        geom_order = 2
        refshape = RefCube
        dim = 3
        #elseif _celltype == "C3D6"
        # Linear wedge
        #elseif _celltype == "C3D15"
        # Quadratic wedge
    else
        throw("Unsupported cell type $_celltype.")
    end
    cells = celltype.(problem.cells)
    nodes = Node.(problem.node_coords)
    grid = Grid(cells, nodes)

    for k in keys(problem.cellsets)
        grid.cellsets[k] = Set(problem.cellsets[k])
    end
    for k in keys(problem.nodesets)
        grid.nodesets[k] = Set(problem.nodesets[k])
    end
    for k in keys(problem.facesets)
        grid.facesets[k] = Set(problem.facesets[k])
    end
    # Define boundary faces
    grid.boundary_matrix = extract_boundary_matrix(grid)

    dh = DofHandler(grid)
    # Isoparametric
    field_interpolation = Lagrange{dim,refshape,geom_order}()
    push!(dh, :u, dim, field_interpolation)
    close!(dh)

    ch = ConstraintHandler(dh)
    for k in keys(problem.nodedbcs)
        vec = problem.nodedbcs[k]
        f(x, t) = [vec[i][2] for i in 1:length(vec)]
        components = [vec[i][1] for i in 1:length(vec)]
        dbc = Dirichlet(:u, getnodeset(grid, k), f, components)
        add!(ch, dbc)
        close!(ch)
        update!(ch, 0.0)
    end

    return ch
end

function extract_boundary_matrix(grid::Grid{dim}) where {dim}
    nfaces = length(Ferrite.faces(grid.cells[1]))
    ncells = length(grid.cells)
    countedbefore = Dict{NTuple{dim,Int},Bool}()
    boundary_matrix = ones(Bool, nfaces, ncells) # Assume all are boundary faces
    for (ci, cell) in enumerate(getcells(grid))
        for (fi, face) in enumerate(Ferrite.faces(cell))
            sface = Ferrite.sortface(face) # TODO: faces(cell) may as well just return the sorted list
            token = Base.ht_keyindex2!(countedbefore, sface)
            if token > 0 # haskey(countedbefore, sface)
                boundary_matrix[fi, ci] = 0
            else # distribute new dofs
                Base._setindex!(countedbefore, true, sface, -token)# countedbefore[sface] = true,  mark the face as counted
            end
        end
    end
    return sparse(boundary_matrix)
end


text:
###########################################################
### Manifold Truncated Newton
###########################################################
"""
    ManNewtonCG

Optimizer performing inexact Newton steps on a prescribed manifold. The inexactness in
the approximation of the Newton step is controlled by `νₖ`: one expects a direction `dᴺ`
such that
        ||∇²(f+g)(x)[dᴺ] + ∇(f+g)|| ≤ νₖ ||∇(f+g)||
"""
@with_kw struct ManNewtonCG <: ManifoldUpdate
    linesearch::ManifoldLinesearch = Armijo()
    # ν_reductionfactor::Float64 = 0.5
end
Base.summary(::ManNewtonCG) = "ManNewtonCG"

@with_kw mutable struct ManNewtonCGState <: AbstractUpdateState
    norm_∇fgₘ::Float64 = -1.0
    d_type::Symbol = :Unsolved
    νₖ::Float64 = 1.0
    CG_niter::Int64 = -1
    CG_ε::Float64 = -1
    CG_residual::Float64 = -1
    norm_dᴺ::Float64 = -1
    cosθ::Float64 = -1
    λ_x::Float64 = -1
    ls_niter::Int64 = -1
end

initial_state(::ManNewtonCG, x, reg) = ManNewtonCGState()

# str_updatelog(o::ManifoldGradient, t::ManifoldGradientState) = @sprintf "ls-nit %2i\t\t||gradₘ f+g|| %.3e" t.ls_niter t.norm_∇fgₘ

function str_updatelog(::ManNewtonCG, t::ManNewtonCGState)
    resstyle = (t.CG_residual <= t.CG_ε) ? "0" : "4;1"
    return string(
        (@sprintf "%-17s CG: nit:%3i  ε:%.1e res:" t.d_type t.CG_niter t.CG_ε),
        "\033[$(resstyle)m",
        (@sprintf "%.1e" t.CG_residual),
        "\033[0m",
        (@sprintf "     |dᴺ|:%.3e   cos(θₖ):%.1e  λ(x):%.1e  ls-nit %2i" t.norm_dᴺ t.cosθ t.λ_x t.ls_niter)
    )
end

function update_iterate!(state::PartlySmoothOptimizerState{Tx}, pb, o::ManNewtonCG) where Tx
    ###
    grad_fgₖ = state.tempvec_man
    M = state.M
    x = state.x.man_repr

    @assert state.x.repr == manifold_repr
    @assert is_manifold_point(M, x)

    state_TN = state.update_to_updatestate[o]

    # @unpack ν_reductionfactor = o
    ν_reductionfactor = 1.0
    @unpack νₖ = state_TN

    ## 1. Get first & second order information
    # TODO: remove intermediate alloc from .+= op.
    grad_fgₖ = egrad_to_rgrad(M, x, state.∇f_x) + ∇M_g(pb, M, x); state.ncalls_gradₘF += 1
    function hessfg_x_h(ξ)
        state.ncalls_HessₘF += 1
        return ehess_to_rhess(M, x, state.∇f_x, ∇²f_h(pb, x, ξ), ξ) + ∇²M_g_ξ(pb, M, x, ξ)
    end

    norm_rgrad = norm(M, x, grad_fgₖ)
    state_TN.norm_∇fgₘ = norm_rgrad

    # check_tangent_vector(M, x, grad_fgₖ)

    ## 2. Get Truncated Newton direction
    # TODO: turn this into approximalte solution finding
    dₖ, state_TN.d_type, state_TN.CG_niter = solve_tCG_capped(M, x, grad_fgₖ, hessfg_x_h; ϵ=1e-15, ζ = 1e-15, maxiter=1e5, check_tvectors=false, printlev=0)

    if state_TN.d_type == :NegativeCurvature
        dₖ = - sign(inner(M, x, dₖ, grad_fgₖ)) * abs(inner(M, x, dₖ, hessfg_x_h(dₖ))) * norm(M, x, dₖ)^(-3) * dₖ
    elseif state_TN.d_type == :MaxIter
    else
        @assert state_TN.d_type == :Solution "Unknown d_type" state_TN.d_type
    end

    # check_tangent_vector(M, x, dᴺ)

    ## 3. Execute linesearch
    # TODO: make linesearch inplace for x, return status.
    hist_ls = Dict()
    x_ls = linesearch(o.linesearch, state, nothing, pb, M, x, grad_fgₖ, dₖ, hist=hist_ls)

    x = x_ls
    state_TN.ls_niter = hist_ls[:niter]

    ## 4. Update CG precision
    ν_strat = :ls
    if ν_strat == :ls
        if hist_ls[:niter] == 1
            state_TN.νₖ *= ν_reductionfactor
        end
    elseif ν_strat == :lsnormgrad
        if hist_ls[:niter] == 1 && νₖ ≥ 1e-2 * norm_rgrad^2
            state_TN.νₖ *= ν_reductionfactor
        end
    else
        @error "unknown reduction strategy"
    end

    state.x.man_repr = x_ls

    ### Logging data
    # state_TN.CG_ε = ϵ_residual
    state_TN.CG_residual = norm(M, x, hessfg_x_h(dₖ) + grad_fgₖ)
    state_TN.norm_dᴺ = norm(M, x, dₖ)
    state_TN.cosθ = inner(M, x, -dₖ, grad_fgₖ)/(norm(M, x, dₖ)*norm_rgrad)
    state_TN.λ_x = inner(M, x, -dₖ, grad_fgₖ)
    state_TN.ls_niter = hist_ls[:niter]

    return
end


text:
# # An equation of second order in two dimmensions
#
# Let \\(\Omega = (0, 1)^3\\) the unit cube and denote \\( \Gamma = \partial \Omega\\) its boundary.
# Consider the following problem. Given \\(f \in C(\Omega)\\), find a function \\(u\\) satisfying
# -\Delta u + u = f in \Omega
# u = 0 on \Gamma

using Revise
using SymFEL
using SymPy
using LinearAlgebra
using SparseArrays
using JLD2
using WriteVTK

## discretization parameters
# we use the mesh square.msh (in gmsh format)
# obtained from square.geo using gmsh
# this mesh is formed by quad elements
# we use gmsh module for read the mesh
# there are some problems using Threads and gmsh (in linux)
# the file is square.jld2 is prepared separetely
@load "cube.jld2"

# elementary matrices -
elem_Mx = SymFEL.get_lagrange_em(1, 0, 0)
elem_Kx = SymFEL.get_lagrange_em(1, 1, 1)

nc = ([1, 2, 2, 1, 1, 2, 2, 1],
      [1, 1, 2, 2, 1, 1, 2, 2],
      [1, 1, 1, 1, 2, 2, 2, 2])
nr = nc

elem_Mxyz = SymFEL.get_cube_em(elem_Mx,
                               elem_Mx,
                               elem_Mx,
                               nc,
                               nr)

elem_Kxyz = SymFEL.get_cube_em(elem_Kx,
                               elem_Mx,
                               elem_Mx,
                               nc,
                               nr) +
            SymFEL.get_cube_em(elem_Mx,
                               elem_Kx,
                               elem_Mx,
                               nc,
                               nr) +
            SymFEL.get_cube_em(elem_Mx,
                               elem_Mx,
                               elem_Kx,
                               nc,
                               nr)

dx = norm(nodes_coordinate[:, elements_bound[1,1]] - nodes_coordinate[:, elements_bound[2,1]])

elem_Kxyz_dx = convert(Matrix{Float64}, elem_Kxyz.subs(h, dx))
elem_Mxyz_dx = convert(Matrix{Float64}, elem_Mxyz.subs(h, dx));

# global matrices
K = SymFEL.assemble_cubemesh_FE_matrix(elem_Kxyz_dx, elements_int, order1=1, order2=1)
M = SymFEL.assemble_cubemesh_FE_matrix(elem_Mxyz_dx, elements_int, order1=1, order2=1)

f = (3*pi^2 + 1) * (sin.(pi * nodes_coordinate[1,:])
                    .* sin.(pi * nodes_coordinate[2,:])
                    .* sin.(pi * nodes_coordinate[3,:]))

F = M * f

A = K + M

# boundary condition
tgv = 1e30
A[nodes_boundary, nodes_boundary] += tgv * sparse(Matrix{Float64}(I, nodes_boundary_N, nodes_boundary_N))

u = A \ F
u_exact = sin.(pi*nodes_coordinate[1,:]) .* sin.(pi*nodes_coordinate[2,:]) .* sin.(pi*nodes_coordinate[3,:])
err = u - u_exact

println("L2 error : ", sqrt(err' * M * err))
println("H1 error : ", sqrt(err' * K * err))

# export to vtk
cells = [MeshCell(VTKCellTypes.VTK_HEXAHEDRON, elements_int[1:8, i]) for i = 1:elements_int_N]

points_x = nodes_coordinate[1, :]
points_y = nodes_coordinate[2, :]
points_z = nodes_coordinate[3, :]
vtkfile = vtk_grid("ex5-output", points_x, points_y, points_z, cells)

vtkfile["u", VTKPointData()] = u
outfiles = vtk_save(vtkfile)


text:
using TuringModels

Turing.setadbackend(:reversediff);
# Turing.turnprogress(false);

d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "Kline.csv"), delim=';');
size(d) # Should be 10x5

# New col log_pop, set log() for population data
d[!, :log_pop] = map(x -> log(x), d[:, :population]);

# New col contact_high, set binary values 1/0 if high/low contact
d[!, :contact_high] = map(x -> ifelse(x=="high", 1, 0), d[:, :contact]);

# This is supposed to be a "bad" model since we take non-centered data for the
# predictor log_pop
@model m10_10stan(total_tools, log_pop, contact_high) = begin
    α ~ Normal(0, 100)
    βp ~ Normal(0, 1)
    βc ~ Normal(0, 1)
    βpc ~ Normal(0, 1)

    for i ∈ 1:length(total_tools)
        λ = exp(α + βp*log_pop[i] + βc*contact_high[i] +
            βpc*contact_high[i]*log_pop[i])
        total_tools[i] ~ Poisson(λ)
    end
end;

posterior = sample(m10_10stan(d[:, :total_tools], d[:, :log_pop],
    d[:, :contact_high]), Turing.NUTS(0.65), 1000);

# Rethinking result

m_10_10t_result = "
     mean   sd  5.5% 94.5% n_eff Rhat
 a    0.94 0.37  0.36  1.53  3379    1
 bp   0.26 0.04  0.21  0.32  3347    1
 bc  -0.08 0.84 -1.41  1.23  2950    1
 bpc  0.04 0.09 -0.10  0.19  2907    1
";

# Describe the draws

describe(chns)

# End of m10.10t.jl


text:
module BranchAndBound

using JuMP
using MathOptInterface

const MOI = MathOptInterface

include("node.jl")
include("tree.jl")

# specific implementation
include("impl-JuMP.jl")

end # module


text:
module Blueprints

using Pkg
using TOML
using Dates
using UUIDs
using LibGit2
using Mustache
using GarishPrint
using Configurations

export Blueprint, MustacheBlueprint

include("types.jl")

include("blueprints/git.jl")
include("blueprints/project_file.jl")

end


text:
module ShaderAbstractions

using StaticArrays, ColorTypes, FixedPointNumbers, StructArrays
import GeometryBasics, GeometryTypes, Tables

include("types.jl")
include("uniforms.jl")
include("api.jl")

end # module


text:
#===============================================================================

    Convenience functions for:
        ✔   packaging together auxiliary objects and observations so
            that the shapes of the two match
        ✔   building recordrings that share the observation schemes

===============================================================================#
"""
    package(objs::AbstractArray, all_obs::AllObservations)

    package(objs::AbstractArray, recording::NamedTuple)

    package(obj, recording::NamedTuple)

    package(objs::AbstractArray, old_to_new_idx::Dict, all_obs::AllObservations)

    package(obj, all_obs::AllObservations)

Various packaging functions that zip together objects/functions/structs together
with observations defined in AllObservations/recording
"""
function package end

function package(objs::AbstractArray, all_obs::AllObservations)
    N = length(objs)
    @assert N == length(all_obs.recordings)
    packaged = []
    for i in 1:N
        push!(packaged, package(objs[i], all_obs.recordings[i]))
    end
    packaged
end

function package(objs::AbstractArray, recording::NamedTuple)
    length(objs) == length(recording.obs) && return objs
    @assert length(objs) == 1
    num_segments = length(recording.obs)
    [deepcopy(obj[1]) for _ in 1:num_segments]
end

function package(obj, recording::NamedTuple)
    num_segments = length(recording.obs)
    [deepcopy(obj) for _ in 1:num_segments]
end

function package(objs::AbstractArray, old_to_new_idx::Dict, all_obs::AllObservations)
    N = length(objs)
    @assert N == length(keys(old_to_new_idx))
    packaged = []
    counter = 1
    for i in 1:N
        for j in 1:length(old_to_new_idx[i])
            num_segments = length(all_obs.recordings[counter].obs)
            push!(packaged, [deepcopy(objs[i]) for _ in 1:num_segments])
            counter += 1
        end
    end
    packaged
end

function package(obj, all_obs::AllObservations)
    N = length(all_obs.recordings)
    packaged = []
    for i in 1:N
        num_segments = length(all_obs.recordings[i].obs)
        push!(packaged, [deepcopy(obj) for _ in 1:num_segments])
    end
    packaged
end

"""
    setup_time_grids(
        all_obs::AllObservations,
        dt=0.01,
        τ=identity,
        eltype=Float64,
        already_arranged_tt=nothing
    )

Same as a version for a single recording, but applies recursively  to all
recordings in `all_obs`.
"""
function setup_time_grids(
        all_obs::AllObservations,
        dt=0.01,
        τ=identity,
        eltype=Float64,
        already_arranged_tt=nothing
    )
    N = length(all_obs.recordings)
    already_arranged_tt === nothing || (
        @assert legnth(already_arranged_tt) == length(all_obs.recordings);
        @assert all([
            length(already_arranged_tt[i]) == length(all_obs.recordings[i].obs)
            for i in 1:N
        ]);
        return already_arranged_tt
    )

    arranged_tt = map(1:N) do i
        _dt = _idx_if_array(dt, i)
        _τ = _idx_if_array(τ, i)
        _eltype = _idx_if_array(eltype, i)

        setup_time_grids(all_obs.recordings[i], _dt, _τ, _eltype)
    end
    arranged_tt
end

_idx_if_array(x, i) = typeof(x) <: AbstractArray ? x[i] : x

"""
    setup_time_grids(
        recording::NamedTuple,
        dt=0.01,
        τ=identity,
        eltype=Float64,
        already_arranged_tt=nothing
    )

Same as a version for a single interval, but applies recursively to all
inter-observation intervals in `recording`.
"""
function setup_time_grids(
        recording::NamedTuple,
        dt=0.01,
        τ=identity,
        eltype=Float64,
        already_arranged_tt=nothing
    )
    N = length(recording.obs)
    already_arranged_tt === nothing || (
        @assert length(already_arranged_tt) == length(recordings.obs);
        return already_arranged_tt
    )
    arranged_tt = map(1:N) do i
        _dt = _idx_if_array(dt, i)
        _τ = _idx_if_array(τ, i)
        _eltype = _idx_if_array(eltype, i)

        _t0 = i == 1 ? recording.t0 : recording.obs[i-1].t
        _T = recording.obs[i].t
        setup_time_grids(_t0, _T, _dt, _τ, _eltype)
    end
    arranged_tt
end

"""
    setup_time_grids(
        t0::Number,
        T::Number,
        dt=0.01,
        τ=identity,
        eltype=Float64,
        already_arranged_tt=nothing
    )

Set up a time grid for a single inter-observation interval [`t0`, `T`]. Starts
with an equidistant time grid with mesh witdh `dt` and applies a transformation
`τ` afterwards. Alternatively, if `already_arranged_tt` is passed then the
initial setting up of an equidistant grid is omitted and `already_arrange_tt` is
used in its place for the `τ` transformation.
"""
function setup_time_grids(
        t0::Number,
        T::Number,
        dt=0.01,
        τ=identity,
        eltype=Float64,
        already_arranged_tt=nothing
    )
    already_arranged_tt === nothing || (
        @assert already_arranged_tt[1] == t0 && already_arranged_tt[end] == T;
        return already_arranged_tt
    )
    _tt = collect(t0:dt:T)
    _tt[end] == T || append!(_tt, [T])
    tt = eltype.(τ(_tt))
end

"""
    num_recordings(all_obs::AllObservations)

Return the total number of recordings in `all_obs`.
"""
num_recordings(all_obs::AllObservations) = length(all_obs.recordings)

"""
    num_obs(all_obs::AllObservations)

Return the total number of observations in all recordings of `all_obs`.
"""
function num_obs(all_obs::AllObservations)
    mapreduce(r->length(r.obs), +, all_obs.recordings)
end


text:
export rule

@rule KernelGCV(:x, Marginalisation) (m_y::MvNormalMeanCovariance, q_z::MvNormalMeanCovariance, meta::KernelGCVMetadata) = begin
    kernelfunction = get_kernelfn(meta)
    Λ_out = approximate_kernel_expectation(get_approximation(meta), (s) -> cholinv(kernelfunction(s)), q_z)
    return MvNormalMeanCovariance(mean(m_y), cov(m_y) + cholinv(Λ_out))
end

@rule KernelGCV(:x, Marginalisation) (m_y::MvNormalMeanPrecision, q_z::MvNormalMeanPrecision, meta::KernelGCVMetadata) = begin
    kernelfunction = get_kernelfn(meta)
    Λ_out = approximate_kernel_expectation(get_approximation(meta), (s) -> cholinv(kernelfunction(s)), q_z)
    return MvNormalMeanPrecision(mean(m_y), cholinv(cov(m_y) + cholinv(Λ_out)))
end


text:
"""
Calculates the SVD deomposition of a matrix in a type generic manner.
Can for example be used with a `Matrix{BigFloat}`.
Returns a `SVD` factorization object.
"""
genericsvd(x) = genericsvd!(copy(x))
export genericsvd
"""
Same as `genericsvd` but saves space by overwriting the input matrix.
"""
function genericsvd!(x)
  GenericSVD.svd!(x)
end
export genericsvd!


text:
# this standalone script can be put anywhere
#
# run from shell/command-line:
# julia fusion.jl -f <REF> -b <BED> -l <READ1> -r <READ2>
#
# required arguments:
#   -f, --ref REF      the reference folder, which contains chr1.fa,chr2fa...,
#                      can be downloaded from
#                      http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
#   -b, --bed BED      a gene list with their coordination intervals
#   -l, --read1 READ1  read1 input fastq file, can be gzipped
#   -r, --read2 READ2  read2 input fastq file, can be gzipped

using FusionDirect
using ArgParse

function main(args)

    s = ArgParseSettings(description = "FusionDirect: detect gene fusion directly from fastq, without alignment required")

    @add_arg_table s begin
        "--ref",  "-f"
            help = "the reference folder, which contains chr1.fa, chr2fa..., can be downloaded from http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz"
            required = true
        "--bed",  "-b"
            help = "a gene list with their coordination intervals, see the example bed files in data folder"
            required = true
        "--read1",   "-l"
            help = "read1 input fastq file, can be gzipped"
            required = true
        "--read2",   "-r"
            help = "read2 input fastq file, can be gzipped"
            required = true
    end

    options = parse_args(s) # the result is a Dict{String,Any}
    detect(options["ref"], options["bed"], options["read1"],options["read2"])
end

# REPL is not supported
if !isinteractive()
    main(ARGS)
end


text:
using FastArrayOps
using Base.Test

include("testingmacros.jl")

for nmax in 20:30,
    i in 1:20,
    inc in [-10:1:-1,1:10]
    @test nel2nmax(i, inc, nmax2nel(i, inc, nmax)) <= nmax
    @test nel2nmax(i, inc, nmax2nel(i, inc, nmax)+1) > nmax
    @test nmax2nel(i, inc, nel2nmax(i, inc, nmax)) == nmax
    @test fast_range2args(fast_args2range(i, inc, nmax)) == (i, inc, nmax)
    @test fast_args2range(fast_range2args(i:inc:nmax)...) == i:inc:nmax
end # for

const TYPES = (Float64, Float32, Complex128, Complex64)
const IX = (1, 2)
const INC1 = (1, 2)
const INC2 = ((1, 1), (2, 2), (2, 3), (-1, 2), (-1, -3))
const INC3 = ((1, 1, 1), (2, 2, 2), (2, 3, 1), (-1, 2, 3), (-1, -3, -2))
const ANUM = 1.123

println("fast_scale!(x, ix, incx, a, n) ...")
for dtype in TYPES,
    ix in IX,
    incx in INC1,
    nmax in (FastArrayOps.NLIM_SCALE-1, FastArrayOps.NLIM_SCALE+1)

    n = nmax2nel(ix, incx, nmax)
    @printtest1(dtype, nmax, ix, incx, n)

    x_exp, x0 = makesignals1(dtype, nmax)
    a = convert(dtype,ANUM)
    x_exp[fast_args2range(ix, incx, n)] *= a

    @testfunc_arr1a(fast_scale!, x_exp, x0, ix, incx, a, n)
end # for

println("fast_scale!(x, ix, incx, y, iy, incy, a, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy) in INC2,
    nmax in (FastArrayOps.NLIM_SCALE_OOP1-1, FastArrayOps.NLIM_SCALE_OOP1+1)

    iy = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    n = min(nx, ny)
    @printtest2(dtype, nmax, ix, incx, iy, incy, n)

    x_exp, x0, y0 = makesignals2(dtype, nmax)
    a = convert(dtype,ANUM)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    x_exp[rx] = a*y0[ry]

    @testfunc_arr2a(fast_scale!, x_exp, x0, ix, incx, y0, iy, incy, a, n)
end # for

println("fast_scale!(x, ix, incx, y, iy, incy, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy) in INC2,
    nmax in (15,22)  # NLIM is max Int

    iy = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    n = min(nx, ny)
    @printtest2(dtype, nmax, ix, incx, iy, incy, n)

    x_exp, x0, y0 = makesignals2(dtype, nmax)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    x_exp[rx] .*= y0[ry]

    @testfunc_arr2(fast_scale!, x_exp, x0, ix, incx, y0, iy, incy, n)
end # for

println("fast_scale!(x, ix, incx, y, iy, incy, z, iz, incz, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy, incz) in INC3,
    nmax in (15,22)  # NLIM is max Int

    iy = ix
    iz = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    nz = nmax2nel(iz, incz, nmax)
    n = min(min(nx, ny), nz)
    @printtest3(dtype, nmax, ix, incx, iy, incy, iz, incz, n)

    x_exp, x0, y0, z0 = makesignals3(dtype, nmax)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    rz = fast_args2range(iz, incz, n)
    x_exp[rx] = y0[ry].*z0[rz]

    @testfunc_arr3(fast_scale!, x_exp, x0, ix, incx, y0, iy, incy, z0, iz, incz, n)
end # for

println("fast_add!(x, ix, incx, a, n) ...")
for dtype in TYPES,
    ix in IX,
    incx in INC1,
    nmax in (15,22)  # NLIM is max Int

    n = nmax2nel(ix, incx, nmax)
    @printtest1(dtype, nmax, ix, incx, n)

    x_exp, x0 = makesignals1(dtype, nmax)
    a = convert(dtype,ANUM)
    x_exp[fast_args2range(ix, incx, n)] += a

    @testfunc_arr1a(fast_add!, x_exp, x0, ix, incx, a, n)
end # for

println("fast_add!(x, ix, incx, y, iy, incy, a, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy) in INC2,
    nmax in (15,22)  # NLIM is max Int

    iy = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    n = min(nx, ny)
    @printtest2(dtype, nmax, ix, incx, iy, incy, n)

    x_exp, x0, y0 = makesignals2(dtype, nmax)
    a = convert(dtype,ANUM)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    x_exp[rx] = y0[ry] + a

    @testfunc_arr2a(fast_add!, x_exp, x0, ix, incx, y0, iy, incy, a, n)
end # for

println("fast_add!(x, ix, incx, y, iy, incy, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy) in INC2,
    nmax in (FastArrayOps.NLIM_ADDARR-1, FastArrayOps.NLIM_ADDARR+1)

    iy = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    n = min(nx, ny)
    @printtest2(dtype, nmax, ix, incx, iy, incy, n)

    x_exp, x0, y0 = makesignals2(dtype, nmax)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    x_exp[rx] .+= y0[ry]

    @testfunc_arr2(fast_add!, x_exp, x0, ix, incx, y0, iy, incy, n)
end # for

println("fast_add!(x, ix, incx, y, iy, incy, z, iz, incz, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy, incz) in INC3,
    nmax in (FastArrayOps.NLIM_ADDARR_OOP1-1, FastArrayOps.NLIM_ADDARR_OOP1+1)

    iy = ix
    iz = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    nz = nmax2nel(iz, incz, nmax)
    n = min(min(nx, ny), nz)
    @printtest3(dtype, nmax, ix, incx, iy, incy, iz, incz, n)

    x_exp, x0, y0, z0 = makesignals3(dtype, nmax)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    rz = fast_args2range(iz, incz, n)
    x_exp[rx] = y0[ry] .+ z0[rz]

    @testfunc_arr3(fast_add!, x_exp, x0, ix, incx, y0, iy, incy, z0, iz, incz, n)
end # for

println("fast_addscal!(x, ix, incx, y, iy, incy, a, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy) in INC2,
    nmax in (FastArrayOps.NLIM_ADDARRSCAL-1, FastArrayOps.NLIM_ADDARRSCAL+1)

    iy = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    n = min(nx, ny)
    @printtest2(dtype, nmax, ix, incx, iy, incy, n)

    x_exp, x0, y0 = makesignals2(dtype, nmax)
    a = convert(dtype,ANUM)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    x_exp[rx] .+= a*y0[ry]

    @testfunc_arr2a(fast_addscal!, x_exp, x0, ix, incx, y0, iy, incy, a, n)
end # for

println("fast_addscal!(x, ix, incx, y, iy, incy, z, iz, incz, a, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy, incz) in INC3,
    nmax in (FastArrayOps.NLIM_ADDARRSCAL_OOP1-1, FastArrayOps.NLIM_ADDARRSCAL_OOP1+1)

    iy = ix
    iz = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    nz = nmax2nel(iz, incz, nmax)
    n = min(min(nx, ny), nz)
    @printtest3(dtype, nmax, ix, incx, iy, incy, iz, incz, n)

    x_exp, x0, y0, z0 = makesignals3(dtype, nmax)
    a = convert(dtype,ANUM)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    rz = fast_args2range(iz, incz, n)
    x_exp[rx] = y0[ry] .+ a*z0[rz]

    @testfunc_arr3a(fast_addscal!, x_exp, x0, ix, incx, y0, iy, incy, z0, iz, incz, a, n)
end # for

println("fast_copy!(x, ix, incx, y, iy, incy, n) ...")
for dtype in TYPES,
    ix in IX,
    (incx, incy) in INC2,
    nmax in (FastArrayOps.NLIM_COPY1-1, FastArrayOps.NLIM_COPY1+1)

    iy = ix
    nx = nmax2nel(ix, incx, nmax)
    ny = nmax2nel(iy, incy, nmax)
    n = min(nx, ny)
    @printtest2(dtype, nmax, ix, incx, iy, incy, n)

    x_exp, x0, y0 = makesignals2(dtype, nmax)
    rx = fast_args2range(ix, incx, n)
    ry = fast_args2range(iy, incy, n)
    x_exp[rx] = y0[ry]

    @testfunc_arr2(fast_copy!, x_exp, x0, ix, incx, y0, iy, incy, n)
end # for

println("fast_fill!(x, ix, incx, a, n) ...")
for dtype in TYPES,
    ix in IX,
    incx in INC1,
    nmax in (FastArrayOps.NLIM_FILL-1, FastArrayOps.NLIM_FILL+1),
    aa in (ANUM, 0)

    n = nmax2nel(ix, incx, nmax)
    @printtest1(dtype, nmax, ix, incx, n)

    x_exp, x0 = makesignals1(dtype, nmax)
    a = convert(dtype,aa)
    x_exp[fast_args2range(ix, incx, n)] = a

    @testfunc_arr1a(fast_fill!, x_exp, x0, ix, incx, a, n)
end # for

println("\nTESTING: SUCCESS")


text:
module Circle_Mod

export Circle_const

mutable struct Circle
    radius
end

function Circle_const(r::Real)
    Circle(r)
end

end # end of module

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
using HardSphereDynamics

using LinearAlgebra
using Test

using StaticArrays

import HardSphereDynamics: normal, collision_time, flow!, collide!

@testset "HardSphereDynamics.jl" begin

    b = MovableBall(SVector(0.0, 0.0), SVector(1.0, 0.0), 2.0)
    @test normal(b, SVector(2.0, 0.0)) == SVector(1.0, 0.0)

    Π = FixedPlane(SVector(10, 0), SVector(1, 0))
    @test collision_time(b, Π) == 8

    b1 = MovableBall(SVector(0.0, 0.0), normalize(SVector(1.0, 0.0)), 1.0)
    b2 = MovableBall(SVector(10.0, 0.5), normalize(SVector(-1.0, 0.0)), 1.0)

    @test collision_time(b1, b2) == 4.031754163448146

    flow!(b1, t)
    flow!(b2, t)

    collide!(b1, b2)

    @test b1.v == SVector(-0.8750000000000002, -0.48412291827592724)

end

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
module Coverings

#--- Dependencies

using ..Lenny: AbstractCovering
import ..Lenny: initialise

#using ..ZeroProblems: dimdeficit, zeroproblem

#--- DefaultCovering - automatically selects an appropriate covering depending on the dimensionality deficit
struct DefaultCovering{T} <: AbstractCovering{T}
end

function initialise(cov::DefaultCovering{T}, prob) where T
	# Choose a covering based on the dimensionality deficit
	dimdef = dimdeficit(zeroproblem(prob))
	if dimdef == 0
		return initialise(Covering0d{T}(), prob)
	elseif dimdef == 1
		return initialise(Covering1d{T}(), prob)
	else
		error("DefaultCovering can only compute 0 or 1 dimensional manifolds")
	end
end

#--- 0 dimensional covering code
struct Covering0d{T} <: AbstractCovering{T}
end

function initialise(cov::Covering0d{T}, prob) where T
	# grab any necessary options from the problem structure
	return cov
end

#--- 1 dimensional covering code
struct Covering1d{T} <: AbstractCovering{T}
end

function initialise(cov::Covering1d{T}, prob) where T
	# grab any necessary options from the problem structure
	return cov
end

end # end module


text:
export SSPKnoth

abstract type RosenbrockAlgorithm <: DistributedODEAlgorithm end

struct RosenbrockTableau{N, RT, N²}
    A::SMatrix{N, N, RT, N²}
    C::SMatrix{N, N, RT, N²}
    Γ::SMatrix{N, N, RT, N²}
    m::SMatrix{N, 1, RT, N}
end

struct RosenbrockCache{Nstages, RT, N², A}
    tableau::RosenbrockTableau{Nstages, RT, N²}
    U::A
    fU::A
    k::NTuple{Nstages, A}
    W
    linsolve!
end

function cache(
    prob::DiffEqBase.AbstractODEProblem,
    alg::RosenbrockAlgorithm; kwargs...)

    tab = tableau(alg, eltype(prob.u0))
    Nstages = length(tab.m)
    U = zero(prob.u0)
    fU = zero(prob.u0)
    k = ntuple(n -> similar(prob.u0), Nstages)
    W = prob.f.jac_prototype
    linsolve! = alg.linsolve(Val{:init}, W, prob.u0; kwargs...)

    return RosenbrockCache(tab, U, fU, k, W, linsolve!)
end

function step_u!(int, cache::RosenbrockCache{Nstages, RT}) where {Nstages, RT}
    tab = cache.tableau

    f! = int.prob.f
    Wfact_t! = int.prob.f.Wfact_t

    u = int.u
    p = int.p
    t = int.t
    dt = int.dt
    W = cache.W
    U = cache.U
    fU = cache.fU
    k = cache.k
    linsolve! = cache.linsolve!

    # 1) compute jacobian factorization
    γ = dt * tab.Γ[1,1]
    Wfact_t!(W, u, p, γ, t)
    for i in 1:Nstages
        U .= u
        for j = 1:i-1
            U .+= tab.A[i,j] .* k[j]
        end
        # TODO: there should be a time modification here (t + c * dt)
        # if f does depend on time, would need to add tgrad term as well
        f!(fU, U, p, t)
        for j = 1:i-1
            fU .+= (tab.C[i,j] / dt) .* k[j]
        end
        linsolve!(k[i], W, fU)
    end
    for i = 1:Nstages
        u .+= tab.m[i] .* k[i]
    end
end

struct SSPKnoth{L} <: RosenbrockAlgorithm
    linsolve::L
end
SSPKnoth(;linsolve)=SSPKnoth(linsolve)

function tableau(::SSPKnoth, RT)
  # ROS.transformed=true;
    N = 3
    N² = N*N
    α = @SMatrix RT[
      0 0 0;
      1 0 0;
      1/4 1/4 0]
    # ROS.d=ROS.alpha*ones(ROS.nStage,1);
    b = @SMatrix RT[1/6 1/6 2/3]
    Γ = @SMatrix RT[
        1 0 0;
        0 1 0;
        -3/4 -3/4 1]
    A = α / Γ
    C = -inv(Γ)
    m = b / Γ
    return RosenbrockTableau{N, RT, N²}(A, C, Γ, m)
#   ROS.SSP.alpha=[1 0 0
#                  3/4 1/4 0
#                  1/3 0 2/3];

end


text:
@testset "IO" begin
  for (l, r) in [ (Dense(20, 30, relu), rand(Float32, 20, 2))
                , (Pointwise(tanh), rand(Float32, 15, 24, 2, 5))
                , (Conv(3, 5), rand(Float32, 7, 7, 3, 12))
                , (Deconv(5, 3), rand(Float32, 8, 8, 5, 4))
                , (Pool(), rand(Float32, 5, 5, 5, 5))
                , (Dropout(0.7, elu), rand(Float32, 7,7,7,7))
                , (Batchnorm(10), rand(Float32, 10, 5)) ]

    ld = l |> Jtac.decompose |> Jtac.compose
    @test typeof(l) == typeof(ld)
    @test all(l(r) .== ld(r))

  end

  for G in [TicTacToe, MetaTac]

    feats = [ConstantFeature([0.7, 0.5, 1.]), ConstantFeature([-5.])]
    chain = @chain G Conv(10) Deconv(7) Batchnorm() Pool() Dense(25) Dropout(0.7)
    model = NeuralModel(G, chain, feats)
    name  = tempname()

    save_model(name, model)
    model2 = load_model(name)
    @test typeof(model) == typeof(model2)

    game = random_turn!(G())
    @test model(game) == model2(game)

    player = MCTSPlayer(model, power = 10)
    dataset = record_self(player, 1, augment = true)

    name = tempname()
    save_dataset(name, dataset)
    ds = load_dataset(name)
    @test all(representation(ds.games) .== representation(dataset.games))
    @test all(ds.label .== dataset.label)
    @test all(ds.flabel .== dataset.flabel)

  end

end

max_stars_repo_name: roSievers/jtac


text:
using SplitApplyCombine
using Dictionaries
using StaticArrays
using Test

include("only.jl")
include("splitdims.jl")
include("combinedims.jl")
include("map.jl")
include("product.jl")
include("group.jl")
include("innerjoin.jl")
include("leftgroupjoin.jl")
include("invert.jl")


text:
using Random

export RandomSource

struct RandomSource
  fs::Float32
  blksize::Int
  nchannels::Int
  streaming::Ref{Bool}
end

function RandomSource(fs; blksize=4096, nchannels=1)
  RandomSource(fs, blksize, nchannels, Ref(false))
end

connect(::RandomSource) = true
Base.close(src::RandomSource) = stopinputstream(src)

inputframerate(src::RandomSource) = src.fs
inputblocksize(src::RandomSource) = src.blksize
inputchannels(src::RandomSource) = src.nchannels

function startinputstream(src::RandomSource, callback)
  src.streaming[] && return
  src.streaming[] = true
  @async begin
    t0 = time()
    n = 0
    data = Matrix{Float32}(undef, src.blksize, src.nchannels)
    while src.streaming[]
      t = n / src.fs
      randn!(data)
      Δt = t0 + t - time()
      Δt > 0 && sleep(Δt)
      callback(round(Int, 1000000t), data)
      n += src.blksize
    end
  end
end

function stopinputstream(src::RandomSource)
  src.streaming[] = false
end


text:
# test scalar field

using Test
using ImageCore, ImageIO
using TensorCast
using Interpolations
using RoME

import IncrementalInference: LevelSetGridNormal

##

@testset "Basic low-res ScalarField localization" begin
##

# # load dem (18x18km span, ~17m/px)
x_min, x_max = -9000, 9000
y_min, y_max = -9000, 9000
# north is regular map image up
global img
x, y, img = RoME.generateField_CanyonDEM(1, 100, x_is_north=false, x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max)

## modify to generate elevation measurements (data/smallData as in Boxy) and priors

dem = Interpolations.LinearInterpolation((x,y), img) # interpolated DEM
elevation(p) = dem(getPPE(fg, p, :simulated).suggested[1:2]'...)
sigma_e = 0.01 # elevation measurement uncertainty

## test buildDEMSimulated

im = (j->((i->dem(i,j)).(x))).(y);
@cast im_[i,j] := im[j][i];
@test norm(im_ - img) < 1e-10

##

function cb(fg_, lastpose)
  global dem, img

  # query DEM at ground truth
  z_e = elevation(lastpose)

  # generate noisy measurement
  @info "Callback for DEM LevelSet priors" lastpose ls(fg_, lastpose) z_e

  # create prior
  hmd = LevelSetGridNormal(img, (x,y), z_e, sigma_e, N=10000, sigma_scale=1)
  pr = PartialPriorPassThrough(hmd, (1,2))
  addFactor!(fg_, [lastpose], pr, tags=[:DEM;], graphinit=false, nullhypo=0.1)
  nothing
end

## Testing

# 0. init empty FG w/ datastore
fg = initfg()
storeDir = joinLogPath(fg,"data")
mkpath(storeDir)
datastore = FolderStore{Vector{UInt8}}(:default_folder_store, storeDir)
addBlobStore!(fg, datastore)

# new feature, going to temporarily disable as WIP
getSolverParams(fg).attemptGradients = false

##

# 1. load DEM into the factor graph
# point uncertainty - 2.5m horizontal, 1m vertical
# horizontal uncertainty chosen so that 3sigma is approx half the resolution
if false
  sigma = diagm([2.5, 2.5, 1.0])
  @time loadDEM!(fg, img, (x), (y), meshEdgeSigma=sigma);
end

##

# 2. generate trajectory

μ0 = [-7000;-2000.0;pi/2]
@time generateCanonicalFG_Helix2DSlew!(10, posesperturn=30, radius=1500, dfg=fg, μ0=μ0, graphinit=false, postpose_cb=cb) #, slew_x=1/20)
deleteFactor!(fg, :x0f1)

##

# ensure specific solve settings
getSolverParams(fg).useMsgLikelihoods = true
getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true

## optional prior at start

mu0 = getPPE(fg, :x0, :simulated).suggested
pr0 = PriorPose2(MvNormal(mu0, 0.01.*[1;1;1;]))
addFactor!(fg, [:x0], pr0)

##

tree = solveTree!(fg);

## check at least the first five poses

for lb in sortDFG(ls(fg,r"x\d+"))[1:4]
  sim = getPPE(fg, lb, :simulated).suggested
  ppe = getPPE(fg, lb).suggested
  @test isapprox(sim[1:2], ppe[1:2], atol=400)
  @test isapprox(sim[3], ppe[3], atol=0.75)
end

##

@error "Skipping latter part of testScalarTest.jl, see #518"
# try

# for lb in sortDFG(ls(fg,r"x\d+"))[5:end]
#   sim = getPPE(fg, lb, :simulated).suggested
#   ppe = getPPE(fg, lb).suggested
#   @test isapprox(sim[1:2], ppe[1:2], atol=400)
#   @test isapprox(sim[3], ppe[3], atol=0.5)
# end

# catch
#   @error "ScalarField test failure on latter half poses"
# end

##
end

##

# using Cairo, RoMEPlotting
# Gadfly.set_default_plot_size(35cm,20cm)

# plotSLAM2D_KeyAndRef(fg)

##


text:
import Base: *

"""
    *(f::Real, g::GP)

Multiplication of a GP `g` from the left by a scalar `f`.
"""
*(f::Real, g::GP) = GP(*, ConstantMean(f), g)
*(f::Function, g::GP) = GP(*, CustomMean(f), g)
μ_p′(::typeof(*), f::MeanFunction, g::GP) = f * mean(g)
k_p′(::typeof(*), f::MeanFunction, g::GP) = OuterKernel(f, kernel(g))
k_pp′(f_p::GP, ::typeof(*), f::MeanFunction, g::GP) = kernel(f_p, g) * f
k_p′p(::typeof(*), f::MeanFunction, g::GP, f_p::GP) = f * kernel(g, f_p)

"""
    *(g::GP, f::Real)

Multiplication of a GP `g` from the right by a scalar `f`.
"""
*(g::GP, f::Real) = GP(*, g, ConstantMean(f))
*(g::GP, f::Function) = GP(*, g, CustomMean(f))
μ_p′(::typeof(*), g::GP, f::MeanFunction) = mean(g) * f
k_p′(::typeof(*), g::GP, f::MeanFunction) = OuterKernel(f, kernel(g))
k_pp′(f_p::GP, ::typeof(*), g::GP, f::MeanFunction) = kernel(f_p, g) * f
k_p′p(::typeof(*), g::GP, f::MeanFunction, f_p::GP) = f * kernel(g, f_p)


text:
# Reference: https://opensees.berkeley.edu/wiki/index.php/Simply_supported_beam_modeled_with_two_dimensional_solid_elements
using UncertaintyQuantification, DelimitedFiles, Formatting

# To run the model distributed add the desired workers and load the required packages with @everywhere
# using Distributed, Formatting
# addprocs(6; exeflags="--project")
# @everywhere using UncertaintyQuantification, DelimitedFiles

E = RandomVariable(Normal(1000, 5), :E)

# Source/Extra files are expected to be in this folder
sourcedir = joinpath(pwd(), "demo/models")

# These files will be rendere through Mustach.jl and have values injected
sourcefiles = ["supported-beam.tcl"]

# These files will be copied to the working directory without injecting values
extrafiles = []

# Dictionary to map FormatSpecs (Formatting.jl) to variables
numberformats = Dict(:E => FormatSpec(".8e"))

# UQ will create subfolders in here to run the solver and store the results
workdir = joinpath(pwd(), "supported-beam")

# Read output file and compute maximum (absolute) displacement
# An extractor is based the working directory for the current sample
disp = Extractor(base -> begin
    file = joinpath(base, "displacement.out")
    data = readdlm(file, ' ')

    return maximum(abs.(data[:, 2]))
end, :disp)

opensees = Solver(
    "OpenSees", # path to OpenSees binary
    "", # extra args passed to the solver binary
    "supported-beam.tcl",
)

ext = ExternalModel(
    sourcedir, sourcefiles, extrafiles, numberformats, workdir, [disp], opensees
)

pf, samples = probability_of_failure(ext, df -> 0.35 .- df.disp, [E], MonteCarlo(1000))

println("Probability of failure: $pf")

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using Plots
gr()

################################
### Algoritmo Newton-Raphson ###
################################

function NewtonRaphson(f, fd, x; tol = 1e-8, max_iter=1000, max_tempo=10.0)
    fx = f(x)               # primeira iteracao
    i = 0                   # contador de iteracoes
    t0 = time()             # tempo de execucao
    t1 = time() - t0
    estado = "Desconhecido" # estado do algoritmo

    # define erro absoluto e erro relativo e calcula o erro
    atol = rtol = tol
    erro = atol + rtol*abs(fx)

    # regras de parada
    rp1 = (abs(fx) ≤ erro)
    rp2 = (t1 ≥ max_tempo || i ≥ max_iter)

    while !(rp1 || rp2)
        fdx = fd(x)

        x = x - fx/fdx
        if abs(fdx) ≤ erro
            estado = "Falha na convergência. Derivada nula encontrada."
            break
        end
        fx = f(x)

        i += 1
        t1 = time() - t0

        rp1 = (abs(fx) ≤ erro)
        rp2 = (t1 ≥ max_tempo || i ≥ max_iter)
    end

    if rp1
        estado = "Convergência Obtida."
    elseif rp2
        i ≥ max_iter ? estado = "Falha na convergência. O Algoritmo atingiu iteração limite" : estado = "Falha na convergência. O algoritmo atingiu tempo limite"
    end
    println("\n",estado,"\n")
    println("Raiz em $x com $i iterações.\n")
    return x, fx, estado
end

Plots.plot(x -> x^2 - 2, -5, 5, leg=false)
f(x) = x^2 - 2
fd(x) = 2x
NewtonRaphson(f, fd, 500.0, tol = 1e-20)
NewtonRaphson(f, fd, -100.0)

f(x) = x * exp(x) - 1
fd(x) = exp(x) + x * exp(x)
Plots.plot(f, -2, 5, leg=false)
NewtonRaphson(f, fd, 0.0)

f(x) = x^2 + 3
fd(x) = 2x
Plots.plot(f, -2, 2, leg=false)
NewtonRaphson(f, fd, 10)

##################################################################
### Newton-Raphson aproximando derivada com diferenças finitas ###
##################################################################

function newton_fin_dif(f, x; h = 0.1, tol = 1e-10, N = 1000)
    c = 0
    for i = 1:N
        c += 1
        xnew = x - (2*h*f(x)) / (f(x+h) - f(x-h))
        if abs(xnew - x) <= tol
            print("root at $x with $(c) iterations.\nError: $(round(xnew-x, sigdigits=5))\n")
            return xnew, c
        end
        x = xnew
    end
    print("WARNING: Function did NOT converge.\nRoot at $x with $(c) iterations.\nError: $(round(xnew-x, sigdigits=5))\n")
    return xnew, c
end

f(x) = 2x^2 - 5x
plot(f, -2, 5, leg=false)
newton_fin_dif(f, 20)

newton_fin_dif(x -> 2x^2 - 5x, 100);
newton_fin_dif(x -> 2x^2 - 5x, -0.5)
newton_fin_dif(x -> 2x^2 - 5x, -500);

plot(x -> x^3 - 2x - 5, -2, 5, leg=false)
newton_fin_dif(x -> x^3 - 2x - 5, 5; h = 0.01, tol=0.0001);

##############################################
### Newton-Raphson Com multiplas dimensões ###
##############################################

function newton_mult_dim(f, G, x; alpha = 0.1, tol = 1e-10, N = 1000)
    c = 0
    for i = 1:N
        c += 1
        xnew = x - alpha*(f(x)./G(x))
        if abs(f(xnew)) <= tol
            print("roots at $xnew with $(c) iterations.\nError: $(round.(xnew-x, sigdigits=5))\n")
            return xnew, c
        end
        x = xnew
    end
    print("WARNING: Function did NOT converge.\nRoots at $x with $(c) iterations.\nError: $(round.(xnew-x, sigdigits=5))\n")
    return xnew, c
end

a1 = range(-20,20,step=0.1)
a2 = range(-20,20,step=0.1)
zplot(x,y) = x^2 + y^2 - 5^2
Plots.surface(a1, a2, zplot)

z(x) = x[1]^2 + x[2]^2 - 5^2
G(x) = [2x[1], 2x[2]]
x0 = [3.0; 5.0]
newton_mult_dim(z, G, x0)

#############################################
### Newton-Raphson usando pacote Roots.jl ###
#############################################

using Roots
f(x) = x^3 - 2x - 5
fp(x) = 3x^2 - 2
newton_fin_dif(f, 2);
NewtonRaphson(f, fp, 10)
x = Roots.newton(f, fp, 2, verbose=true)
x = find_zero((f, fp), 2, Roots.Newton(), verbose=true)
x, f(x)


text:
###############

function loss(model, S, c, u, λ, cTrue)
    return Flux.Losses.mse(S*vec(model(c)), u) +  λ*Flux.Losses.mse(vec(model(c)), vec(cTrue))
  end

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

  function reconstructionDIP(S, N, u; model = make_model_unet(n),
                                      opt = ADAM(0.001, (0.7, 0.999)),
                                      cTrue = ones(Float32,N...),
                                      noiseData = vec(rand(Float32,N...).*Float32(0.7)),
                                      iterations = 100,
                                      λ = Float32(0),
                                      iters = nothing)

    parameters = params(model)
    data = [(model, S, noiseData, u, λ, device(cTrue))]; # λ, device(cTrue)

    @info "training"
    iter = 1:iterations
    #@showprogress 1.0 "Computing..."
    residual = Float64[]
    error = Float64[]
    ssims = Float64[]
    for i=iter
      Flux.train!(loss, parameters, data, opt)
      currentC = model(noiseData) |> cpu
      push!(residual, loss(data[1]...))
      push!(error, norm(vec(currentC)-vec(cTrue))/norm(vec(cTrue)) )
      push!(ssims, ssim(squeeze(currentC),squeeze(cTrue)) )
      if iters != nothing
        push!(iters, currentC)
      end
      if mod(i,20) == 1
        @info "res=$(residual[end]) err=$(error[end]) ssim=$(ssims[end]) time=$(i/iter[end])"
      end
    end
    cTrained = model(noiseData) |> cpu

    return cTrained, residual, error, ssims
  end

  function lossPre(model, c, u)
    return Flux.Losses.mae(vec(model(c)), vec(u))
  end

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

  function reconstructionDIPPre(N; model = make_model_unet(n),
                                      opt = ADAM(0.001, (0.7, 0.999)),
                                      cTrue = ones(Float32,N...),
                                      noiseData = vec(rand(Float32,N...).*Float32(0.7)),
                                      iterations = 100)

    parameters = params(model)
    data = [(model, noiseData, cTrue)];

    @info "training"
    iter = 1:iterations
    #@showprogress 1.0 "Computing..."
    residual = Float64[]

    for i=iter
      Flux.train!(lossPre, parameters, data, opt)

      push!(residual, lossPre(data[1]...))

      if mod(i,20) == 1
        @info "res=$(residual[end])  time=$(i/iter[end])"
      end
    end
    cTrained = model(noiseData) |> cpu

    return cTrained, residual
  end


text:
module Polyorder

using CxxWrap
using YAML
using Optim

if Sys.isapple()
    @wrapmodule(joinpath(@__DIR__, "../lib/libpolyorderjl.dylib"))
elseif Sys.islinux()
    # @wrapmodule(joinpath(@__DIR__, "../lib/libpolyorderjl.so"))
    error("Linux is currently not supported by Polyorder. Loading library failed!")
elseif Sys.iswindows()
    # @wrapmodule(joinpath(@__DIR__, "../lib/libpolyorderjl.dll"))
    error("Windows is currently not supported by Polyorder. Loading library failed!")
else
    error("OS is not supported by Polyorder. Loading library failed!")
end

include("common.jl")
include("chain2component2.jl")

function __init__()
    @initcxx
end

end


text:
# ---
# title: 331. Verify Preorder Serialization of a Binary Tree
# id: problem331
# author: Qling
# date: 2021-03-12
# difficulty: Medium
# categories: Stack
# link: <https://leetcode.com/problems/verify-preorder-serialization-of-a-binary-tree/description/>
# hidden: true
# ---
#
# One way to serialize a binary tree is to use pre-order traversal. When we
# encounter a non-null node, we record the node's value. If it is a null node,
# we record using a sentinel value such as `#`.
#
#
#
#          _9_
#         /   \
#        3     2
#       / \   / \
#      4   1  #  6
#     / \ / \   / \
#     # # # #   # #
#
#
# For example, the above binary tree can be serialized to the string
# `"9,3,4,#,#,1,#,#,2,#,6,#,#"`, where `#` represents a null node.
#
# Given a string of comma separated values, verify whether it is a correct
# preorder traversal serialization of a binary tree. Find an algorithm without
# reconstructing the tree.
#
# Each comma separated value in the string must be either an integer or a
# character `'#'` representing `null` pointer.
#
# You may assume that the input format is always valid, for example it could
# never contain two consecutive commas such as `"1,,3"`.
#
# **Example 1:**
#
#
#
#     Input:"9,3,4,#,#,1,#,#,2,#,6,#,#"
#     Output:true
#
# **Example 2:**
#
#
#
#     Input:"1,#"
#     Output:false
#
#
# **Example 3:**
#
#
#
#     Input:"9,#,#,1"
#     Output:false
#
#
## @lc code=start
using LeetCode

function is_valid_serialization(preorder::String)::Bool
    p = split(preorder, ",")
    slot = 1

    for node in p
        slot == 0 && return false

        slot = (node == "#") ? (slot - 1) : (slot + 1)
    end

    return slot == 0
end

## @lc code=end


text:
using JLBoost
using Test
using DataFrames

@testset "smoke test" begin
    df = DataFrame(x=[1,1,1,0,0], y = [1,1,1,0,0])
    jlboost(df, :y; nrounds=4)
end


text:
include("lib_for_robot.jl")
function task10(r)
    side=Ost
    n=0
    k=0
    while isborder(r,Nord)==false
        a,b=moveToSide(r,side)
        k+=a
        n+=b
        move!(r,Nord)
        side=reverse(side)
    end
    a,b=moveToSide(r,side)
    k+=a
    n+=b
    print(k/n)
end


text:
using Plots

#=
Another OEIS function. This one calculates the length of binary as a number continues upwards!

Plots wasn't working correctly here for some reason...
  [Update]...a restart of Julia worked apparently.
=#

println("Enter Maximum Number Please")

(function ()
  maxNum = parse(UInt64, readline())
  mod = 0
  len = 1
  arr = Array{Int64}(undef, maxNum, 1)

  Threads.@threads for i in 1:maxNum
    if (i == len)
      mod += 1
      len *= 2
    end

    arr[i] = mod
    println(i," ",mod)
  end

  println("Would you like to see the plot? [y or yes]")
  if (readline() == "y" || readline() == "yes")
    gui(plot(1:maxNum,arr))

  println("Would you like to save the plot? [y or yes]")
    if (readline() == "y" || readline() == "yes")

      println("Directory? [Use only backslashes]")
      png(readline())
    end
  end

end )()

exit()


text:
import Base.(*)

dict_types(d::Dict{S,T}) where {S,T} = (S,T)

"""
For dictionaries `f` and `g`, the composition `f*g` is a new
dictionary with the property that `(f*g)[x]` yields `f[g[x]]`.
"""
function *(f::Dict, g::Dict)
  A,B = dict_types(f)
  C,D = dict_types(g)

  KK = keys(g)
  faults::Int = 0   # Number of bad keys
  fg = Dict{C,B}()

  if A != D
    warn("Dictionary type mismatch")
  end

  for x in KK
    y = g[x]
    if haskey(f,y)
      fg[x] = f[y]
    else
      faults += 1
    end
  end

  if faults > 0
    warn("$faults keys were not mapped")
  end

  return fg
end


text:
@testset "Constructors" begin
    block = DCTBlock{Float32}(1:64)
    @test block[3, 7] == 42.0f0
    @test block[42] == 42.0f0
end

@testset "Traits" begin
    block = DCTBlock{Float32}(1:64)
    @test sizeof(DCTBlock{Float32}) == 64 * 4
    @test sizeof(block) == 64 * 4
    @test length(block) == 64
    @test size(block) == (8, 8)
end

@testset "Functions" begin
    block = DCTBlock{Float32}(1:64)
    mat = collect(block)
    @test mat isa Matrix{Float32}
    @test mat[3, 7] == 42.0f0
end


text:
using Test, GlobalSearchRegression, DataFrames, Distributions

data = DataFrame(Array{Union{Missing,Float64}}(randn(100,6)))
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data,2) - 1 ] ]
names!(data, headers )

rename!(data, :x5 => Symbol("weird_name"))

data[1, 2] = missing
data[4, 5] = missing

############# EQUATION #################

# unsorted (x2 -the firstone- is the depvar)
gsreg("x2 x1 y", data)

# Stata like
gsreg("y x1 x2 x3", data)

# Stata like with comma
gsreg("y,x1,x2,x3", data)

# Unconventional varname
gsreg("y x1 weird_name", data)

# R like
gsreg("y ~ x1 + x2 + x3", data)
gsreg("y ~ x1 + x2 + x3", data = data)

# Array of strings
gsreg(["y", "x1", "x2", "x3"], data)

# Also, with wildcard
gsreg("y *", data)

gsreg("y x*", data)

gsreg("y ~ x*", data)
gsreg("y ~ .", data)

################ DATA #################

# dataframe with implicit datanames
gsreg("x2 x1 y", data)

headers = map(e->string(e),headers)

# array with explicit datanames
gsreg("x2 x1 y", convert(Array,data); datanames=headers)

# tuple with explcit datanames
gsreg("x2 x1 y", (convert(Array,data), headers) )

############### ARGUMENTS ##############

gsreg("x2 x1 y", data; intercept=false) # without constant

gsreg("x2 x1 y", data; outsample=10) # expect 88 obs

gsreg("x2 x1 y", data; outsample=10, criteria=[:r2adj, :bic, :aic, :aicc, :cp, :rmse, :rmseout, :sse])

gsreg("x2 x1 y", data; ttest=true)

gsreg("x2 x1 y", data; vectoroperation=true)

gsreg("x2 x1 y", data; modelavg=true)

gsreg("x2 x1 y", data; residualtest=true)

gsreg("x2 x1 y", data; time=:x4)
gsreg("x2 x1 y", convert(Array,data); datanames=headers, time=:x4)
gsreg("x2 x1 y", (convert(Array,data), headers); time=:x4)

gsreg("x2 x1 y", data; summary="summary.txt")

gsreg("x2 x1 y", data; orderresults=true)

gsreg("x2 x1 y", data; onmessage= message -> (println(message)))

gsreg("x2 x1 y", data; method="fast")

gsreg("x2 x1 y", data; method="precise")

gsreg("x2 x1 y", data; csv="results.csv")

res = gsreg("x2 x1 y", data; resultscsv="results.csv")

@test size(res.results,1) == 3

data[:x2] = data[:x3]*2
@test_throws ErrorException gsreg("y x1 x2 x3", data)


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Random
using ElasticArrays
using ArraysOfArrays
using FillArrays
import TypedTables

@testset "abstract_value_shape" begin
    @testset "default_datatype" begin
        @test @inferred(ValueShapes.default_datatype(Integer)) == Int
        @test @inferred(ValueShapes.default_datatype(AbstractFloat)) == Float64
        @test @inferred(ValueShapes.default_datatype(Real)) == Float64
        @test @inferred(ValueShapes.default_datatype(Complex)) == Complex{Float64}

        @test @inferred(elshape(Complex(1.0, 2.0))) == ScalarShape{Complex{Float64}}()
        @test @inferred(elshape([[3, 5], [3, 2]])) == ArrayShape{Int,1}((2,))

        @test ValueShapes.stripscalar(Ref(ScalarShape{Real})) == ScalarShape{Real}

        # @test Vector{Real}(undef, ArrayShape{Real}((2,1))) ==  # weird typing going on with default_unshaped_eltype

        arrshape = ArrayShape{Real, 2}((2,3))
        vec = Vector{Real}(undef, arrshape)
        @test @inferred(length(vec)) == 6
        @test @inferred(size(vec)) == (6,)

        data1 = [1;2;3;4;7;8;9]
        scalarshape = ScalarShape{Real}()
        ntshape = NamedTupleShape(a=arrshape, b=scalarshape)
        shapedasnt = ntshape(data1)
        @test stripscalar(Ref(shapedasnt)) == Ref(shapedasnt)[]

        @test_throws ArgumentError Broadcast.broadcastable(ntshape)

        named_shapes = (
            a = ArrayShape{Real}(2, 3),
            b = ScalarShape{Real}(),
            c = ConstValueShape(4.2),
            x = ConstValueShape([11 21; 12 22]),
            y = ArrayShape{Real}(4)
        )
        shape = NamedTupleShape(;named_shapes...)
        data2 = VectorOfSimilarVectors(reshape(collect(1:22), 11, 2))
        @test_throws ArgumentError ValueShapes._checkcompat_inner(ntshape, data2)
        @test ValueShapes._checkcompat_inner(shape, data2) == nothing

        @test @inferred(unshaped(4.2)) isa Fill{Float64,1}
        @test unshaped(4.2) == [4.2]
        @test @inferred(unshaped(view([4.2], 1))) isa SubArray{Float64,1,Vector{Float64}}
        @test unshaped(view([4.2], 1)) == [4.2]
        @test @inferred(unshaped(Array(view([4.2], 1)))) isa SubArray{Float64,1,Vector{Float64}}
        @test unshaped(Array(view([4.2], 1))) == [4.2]
        let x = rand(15)
            @test @inferred(unshaped(x)) === x
            @test @inferred(unshaped(Base.ReshapedArray(x, (3, 5), ()))) === x

            @test @inferred(broadcast(unshaped, x)) isa ArrayOfSimilarArrays{Float64,1,1,2,<:Base.ReshapedArray}
            @test broadcast(unshaped, x) == nestedview(reshape(x, 1, 15))
        end

        let A = rand(1,15)
            @test @inferred(broadcast(unshaped, view(A, 1, :))) isa ArrayOfSimilarArrays{Float64,1,1,2,<:SubArray}
            @test broadcast(unshaped, view(A, 1, :)) == nestedview(A)
        end
    end
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to estimate Caliendo and Parro (2015) trade elasticities
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, RData, XLSX, LinearAlgebra, Statistics, CSV

dir = "X:/VIVES/1-Personal/Florian/git/IOTables/src/"
include(dir * "model/transform_data.jl") # Script with functions to import and transform raw data
# include(dir * "elasticity/elasticities_functions.jl") # Script with functions to compute statistics for estimating sectoral trade elasticities (method of Caliendo and Parro, 2015)
# include(dir * "elasticity/wto_tariffs_2018.jl") # Script to prepare WTO tariff data from 2018

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Data specification
dir_raw = "C:/Users/u0148308/data/raw/" # location of raw data

# WIOD rev. 2013
source = "WIOD"
revision = "2013"
year = 2011 # specified year
N = 41 # number of countries
S = 35 # number of industries

# # WIOD rev. 2016
# source = "WIOD"
# revision = "2016"
# year = 2014 # specified year
# N = 44 # number of countries
# S = 56 # number of industries

# # OECD rev. 2021
# source = "OECD"
# revision = "2021"
# year = 1995 # specified year
# N = 71 # number of countries
# S = 45 # number of industries

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F = transform_data(dir_raw, source, revision, year, N, S)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

ctry_names = ["AUS", "AUT", "BEL", "BGR", "BRA", "CAN", "CHN", "CYP", "CZE", "DEU", "DNK", "ESP", "EST", "FIN", "FRA", "GBR",
    "GRC", "HUN", "IDN", "IND", "IRL", "ITA", "JPN", "KOR", "LTU", "LUX", "LVA", "MEX", "MLT", "NLD", "POL", "PRT", "ROM", "RUS",
    "SVK", "SVN", "SWE", "TUR", "TWN", "USA"]

ctry_EU = ["AUT", "BEL", "BGR", "CYP", "CZE", "DNK", "EST", "FIN", "FRA", "DEU", "GBR", "GRC", "HUN", "IRL", "ITA", "LVA",
    "LTU", "LUX", "MLT", "NLD", "POL", "PRT", "ROM", "SVK", "SVN", "ESP", "SWE"] # with GBR but without HVR

countries = [ctry_names; "ZoW"]

industries = lpad.(1:S, 2, '0') # left pad with leading zeros so we can sort later on

# -------

function index_subset(N::Integer, S::Integer,
    ctry_all::Vector{String}, ctry_subset::Vector{String}, ind_all::Vector{String}, ind_subset::Vector{String})

    # subsetting rows
    N_subset = length(ctry_subset)
    S_subset = length(ind_subset)

    all = repeat(ctry_all, inner=S) .* "_" .* repeat(ind_all, outer=N) # all row country-industry pairs
    subset = repeat(ctry_subset, inner=S_subset) .* "_" .* repeat(ind_subset, outer=N_subset) # subset of country-industry pairs
    index_subset = findall(in(subset), all)

    return index_subset
end

function ctry_aggregation(Z::Matrix, F::Matrix, N::Integer, S::Integer,
    ctry_all::Vector{String}, ctry_subset::Vector{String}, ind_all::Vector{String}, ind_subset::Vector{String})

    # obtain indices for subsetting Z and F
    index_row_subset = index_subset(N, S, ctry_all, ctry_subset, ind_all, ind_subset)
    index_col_subset = index_subset(N, S, ctry_all, ctry_subset, ind_all, ind_subset)
    index_col_subset_F = findall(in(ctry_subset), ctry_all) # only have countries in columns

    # subset columns, construct aggregate and rebuild matrix
    N_col_subset = length(ctry_subset)
    S_col_subset = length(ind_subset)

    Z_col_subset = Z[:, index_col_subset]
    Z_col_subset = [sum(Z_col_subset[i,j:S_col_subset:(N_col_subset-1)*S_col_subset+j]) for i in 1:N*S, j in 1:S_col_subset] # sum over industries to create aggregate
    Z_col_other = Z[:, Not(index_col_subset)] # other countries
    Z_new = hcat(Z_col_other, Z_col_subset) # rebuild Z with new aggregate as last columns

    F_col_subset = F[:, index_col_subset_F]
    F_col_subset = [sum(F_col_subset[i,:]) for i in 1:N*S]
    F_col_other = F[:, Not(index_col_subset_F)] # other countries
    F_new = hcat(F_col_other, F_col_subset) # rebuild F with new aggregate as last column

    # subset rows from new matrix, construct aggregate and rebuild matrix
    N_row_subset = length(ctry_subset)
    S_row_subset = length(ind_subset)
    N_new = size(F_new, 2) # new number of countries

    Z_row_subset = Z_new[index_row_subset, :]
    Z_row_subset = [sum(Z_row_subset[i:S_row_subset:(N_row_subset-1)*S_row_subset+i,j]) for i in 1:S_row_subset, j in 1:N_new*S_row_subset]
    Z_row_other = Z_new[Not(index_row_subset), :]
    Z_new = vcat(Z_row_other, Z_row_subset)

    F_row_subset = F_new[index_row_subset, :]
    F_row_subset = [sum(F_row_subset[i:S_row_subset:(N_row_subset-1)*S_row_subset+i,j]) for i in 1:S_row_subset, j in 1:N_new]
    F_row_other = F_new[Not(index_row_subset), :]
    F_new = vcat(F_row_other, F_row_subset)

    return Z_new, F_new
end

a, b = ctry_aggregation(Z, F, N, S, countries, ctry_EU, industries, industries)

function ctry_aggregation(Z::Matrix, F::Matrix, N::Integer, S::Integer, ctry_all::Vector{String}, ctry_subset::Vector{String})

    # obtain indices for subsetting Z and F
    N_subset = length(ctry_subset)
    industries = lpad.(1:S, 2, '0') # vector of industries

    all = repeat(ctry_all, inner=S) .* "_" .* repeat(industries, outer=N) # all row country-industry pairs
    subset = repeat(ctry_subset, inner=S) .* "_" .* repeat(industries, outer=N_subset) # subset of country-industry pairs

    index_subset = findall(in(subset), all) # can be used for both rows and columns since IO table is a square matrix
    index_subset_F = findall(in(ctry_subset), ctry_all) # only have countries in columns

    # subset columns, construct aggregate and rebuild matrix
    Z_col_subset = Z[:, index_subset]
    Z_col_subset = [sum(Z_col_subset[i, j:S:(N_subset-1)*S+j]) for i in 1:N*S, j in 1:S] # sum over industries to create aggregate
    Z_col_other = Z[:, Not(index_subset)] # other countries
    Z_new = hcat(Z_col_other, Z_col_subset) # rebuild Z with new aggregate as last columns

    F_col_subset = F[:, index_subset_F]
    F_col_subset = [sum(F_col_subset[i,:]) for i in 1:N*S]
    F_col_other = F[:, Not(index_subset_F)] # other countries
    F_new = hcat(F_col_other, F_col_subset) # rebuild F with new aggregate as last column

    # subset rows from new matrix, construct aggregate and rebuild matrix
    N_new = size(F_new, 2) # new number of countries

    Z_row_subset = Z_new[index_subset, :]
    Z_row_subset = [sum(Z_row_subset[i:S:(N_subset-1)*S+i, j]) for i in 1:S, j in 1:N_new*S]
    Z_row_other = Z_new[Not(index_subset), :]
    Z_new = vcat(Z_row_other, Z_row_subset)

    F_row_subset = F_new[index_subset, :]
    F_row_subset = [sum(F_row_subset[i:S:(N_subset-1)*S+i, j]) for i in 1:S, j in 1:N_new]
    F_row_other = F_new[Not(index_subset), :]
    F_new = vcat(F_row_other, F_row_subset)

    return Z_new, F_new
end

c, d = ctry_aggregation(Z, F, N, S, countries, ctry_EU)

println(a == c)
println(b == d)


text:
module CUTENSOR

using ..CuArrays
using ..CuArrays: libcutensor, @libcutensor, active_context

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CEnum
const cudaDataType_t = cudaDataType

include("libcutensor_common.jl")
include("error.jl")

function version()
    ver = cutensorGetVersion()
    major, ver = divrem(ver, 10000)
    minor, patch = divrem(ver, 100)

    VersionNumber(major, minor, patch)
end

include("libcutensor.jl")
include("highlevel.jl")
include("wrappers.jl")

const _handles = Dict{CuContext,cutensorHandle_t}()
const _handle = Ref{cutensorHandle_t}(C_NULL)

function handle()
    if _handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cutensorCreate()
            atexit(()->CUDAdrv.isvalid(context) && cutensorDestroy(handle))
            handle
        end
    end
    return _handle[]
end

end