text:
add_verbose_scope(:BrauerObst)

###############################################################################
#
#  Function to check the cocycle condition
#
###############################################################################

function check_cocycle(G::Array{NfToNfMor, 1}, Coc::Function, d::Int)

  K = domain(G[1])
  O = maximal_order(K)
  T, mT = torsion_unit_group(O)
  zeta = mT(divexact(order(T), d)*T[1])
  @assert zeta^d == 1
  for g1 in G
    zeta1 = O(g1(K(zeta)))
    # I need the discrete logarithm with respect to zeta
    i = 1
    a = deepcopy(zeta)
    while a != zeta1
      a *= zeta
      i += 1
    end
    expzeta1 = i
    @assert zeta^i == zeta1
    for g2 in G
      for g3 in G
        lhs = mod(Coc(g1, g2*g3) + expzeta1*Coc(g2, g3), d)
        rhs = mod(Coc(g1*g2, g3) + Coc(g1, g2), d)
        @assert lhs == rhs
      end
    end
  end
  return true

end

###############################################################################
#
#  P-Sylow subgroup
#
###############################################################################

function pSylow(Gperm::Main.ForeignGAP.MPtr, permGAP::Vector{Main.ForeignGAP.MPtr}, G::Vector{NfToNfMor}, p::Int)
  p2 = ispower(length(G))[2]
  if p == p2
    return G
  end
  H = GAP.Globals.SylowSubgroup(Gperm, p)
  lGp = GAP.Globals.Size(H)
  Gp = Array{Hecke.NfToNfMor,1}(undef, lGp)
  j = 1
  for ind = 1:length(G)
    if j > lGp
      break
    end
    if GAP.Globals.IN(permGAP[ind], H)
      Gp[j] = G[ind]
      j += 1
    end
  end
  return Gp
end

function pSylow(G::Vector{NfToNfMor}, p::Int)
  com, uncom = ppio(length(G), p)
  if uncom == 1
    return G
  end
  permGC = _from_autos_to_perm(G)
  Gperm = _perm_to_gap_grp(permGC)
  PermGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(permGC))
  for w = 1:length(permGC)
    PermGAP[w] = _perm_to_gap_perm(permGC[w])
  end
  return pSylow(Gperm, PermGAP, G, p)
end

###############################################################################
#
#  Cocycle computation: prime case
#
###############################################################################

function _find_exp(x::Main.ForeignGAP.MPtr, y::Main.ForeignGAP.MPtr)
  # I want to find i such that x^i = y
  i = 1
  z = GAP.Globals.ShallowCopy(x)
  while z != y
    z = z*x
    i += 1
  end
  #@assert x^i == y
  return i
end

function _autos_to_check(G::Main.ForeignGAP.MPtr, E::Main.ForeignGAP.MPtr, mG::Main.ForeignGAP.MPtr)
  AutG = GAP.Globals.AutomorphismGroup(G)
  AutE = GAP.Globals.AutomorphismGroup(E)
  #I want to construct the map between the automorphism groups. The kernel is charachteristic!
  gens = GAP.Globals.GeneratorsOfGroup(AutE)
  ind_auts = Array{Main.ForeignGAP.MPtr, 1}(undef, length(gens))
  for s = 1:length(gens)
    ind_auts[s] = GAP.Globals.InducedAutomorphism(mG, gens[s])
  end
  S = GAP.Globals.Subgroup(AutG, GAP.julia_to_gap(ind_auts))
  T = GAP.Globals.List(GAP.Globals.RightTransversal(AutG, S))
  return GAP.gap_to_julia(T)
end

function cocycle_computation(L::Main.ForeignGAP.MPtr, i::Int)

  proj = GAP.Globals.NaturalHomomorphismByNormalSubgroup(L[1], L[i+1])
  target_grp = GAP.Globals.ImagesSource(proj)
  mH1 = GAP.Globals.NaturalHomomorphismByNormalSubgroup(target_grp, GAP.Globals.Image(proj, L[i]))
  H1 = GAP.Globals.ImagesSource(mH1)
  K = GAP.Globals.Kernel(mH1)
  #@assert (GAP.Globals.IsCyclic(K)).ptr == GAP.Globals.ReturnTrue().ptr
  # The kernel is cyclic by hypothesis, but Gap represent it as a pc group with more than one generator.
  # Brute force at the moment.
  gen = GAP.Globals.MinimalGeneratingSet(K)[1]

  Elems = GAP.Globals.Elements(H1)
  MatCoc = Array{Int, 2}(undef, length(Elems), length(Elems))
  Preimags = Vector{Main.ForeignGAP.MPtr}(undef, length(Elems))
  for i = 1:length(Elems)
    Preimags[i] = GAP.Globals.PreImagesRepresentative(mH1, Elems[i])
  end
  for i = 1:length(Elems)
    x1 = Preimags[i]
    for j = 1:length(Elems)
      y1 = Preimags[j]
      xy = Elems[i] * Elems[j]
      k = 1
      while Elems[k] != xy
        k += 1
      end
      xy1 = Preimags[k]
      MatCoc[i,j] = _find_exp(gen, x1*y1*GAP.Globals.Inverse(xy1))
    end
  end

  function cocycle(x::Main.ForeignGAP.MPtr, y::Main.ForeignGAP.MPtr)
    i = 1
    while Elems[i] != x
      i += 1
    end
    j = 1
    while Elems[j] != y
      j += 1
    end
    return MatCoc[i, j]
  end

  autos = _autos_to_check(H1, target_grp, mH1)
  return mH1, autos, cocycle

end

function _reduce_to_prime(L::Main.ForeignGAP.MPtr, i::Int, d::Int)

  p = ispower(d)[2]
  G = GAP.Globals.FactorGroup(L[1], L[i+1])
  L1 = GAP.Globals.DerivedSeries(G)
  mH1 = GAP.Globals.NaturalHomomorphismByNormalSubgroup(G, L1[end-1])
  K = GAP.Globals.Kernel(mH1)
  gens = GAP.Globals.MinimalGeneratingSet(K)
  @assert length(gens) == 1
  gens[1] = gens[1]^divexact(d, p)
  S = GAP.Globals.Subgroup(G, gens)
  Q = GAP.Globals.FactorGroup(G, S)
  return GAP.Globals.DerivedSeries(Q)

end

###############################################################################
#
#  Brauer obstruction using Crossed Product Algebras: degree 2 case
#
###############################################################################

function _Brauer_at_two(list::Vector{FieldsTower}, L::Main.ForeignGAP.MPtr, i::Int)
  mH, autos, _cocycle_values = cocycle_computation(L, i)
  domcoc = GAP.Globals.ImagesSource(mH)
  admit_ext = falses(length(list))
  for t = 1:length(list)
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Fields to test: $(length(list)-t+1)"
    K = list[t].field
    GC = automorphisms(K, copy = false)
    DGCvect = Vector{Tuple{NfToNfMor, Int}}(undef, length(GC))
    for i = 1:length(GC)
      DGCvect[i] = (GC[i], i)
    end
    DGC = Dict{NfToNfMor, Int}(DGCvect)
    permGC = _from_autos_to_perm(GC) #TODO: Improve a little.
    Gperm = _perm_to_gap_grp(permGC)
    PermGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(permGC))
    for w = 1:length(permGC)
      PermGAP[w] = _perm_to_gap_perm(permGC[w])
    end
    #Restrict to the p-Sylow
    Gp = pSylow(Gperm, PermGAP, GC, 2)
    iso = GAP.Globals.IsomorphismGroups(Gperm, domcoc)
    ElemGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(permGC))
    for w = 1:length(permGC)
      ElemGAP[w] = GAP.Globals.Image(iso, PermGAP[w])
    end
    list[t].imgs_autos = ElemGAP
    obstruction = false
    auts_for_conductors = Vector{Main.ForeignGAP.MPtr}()
    for g in autos
      #I create the cocycle
      local cocycle
      let DGC = DGC, g = g, ElemGAP = ElemGAP
        function cocycle(aut1::NfToNfMor, aut2::NfToNfMor)
          s1 = DGC[aut1]
          a1 = GAP.Globals.Image(g, ElemGAP[s1])
          s2 = DGC[aut2]
          b1 = GAP.Globals.Image(g, ElemGAP[s2])
          return mod(_cocycle_values(a1, b1), 2)
        end
      end

      #Now, we have to see if it splits
      fl = cpa_issplit(K, Gp, cocycle, 2)
      if fl
        obstruction = true
        push!(auts_for_conductors, g)
      end

    end
    if obstruction
      list[t].proj_ext = mH
      list[t].auts_for_conductors = auts_for_conductors
      admit_ext[t] = true
    end
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  return list[findall(admit_ext)]

end

###############################################################################
#
#  Brauer Obstruction: prime case
#
###############################################################################

function _Brauer_prime_case(list::Vector{FieldsTower}, L::Main.ForeignGAP.MPtr, i::Int, p::Int)

  list1 = falses(length(list))
  mH, autos, coc = cocycle_computation(L, i)
  H = GAP.Globals.ImagesSource(mH)
  for t = 1:length(list)
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Fields to test: $(length(list)-t+1)"
    Kc = _ext_cycl(list[t].generators_of_automorphisms, p)
    K = list[t].field
    K1 = Kc.Ka
    autsK = automorphisms(K, copy = false)
    autsK1 = automorphisms(K1, copy = false)
    restr = restriction(autsK1, autsK, Kc.mp[2])
    #I construct the group and the isomorphisms between the automorphisms and the gap group.
    permGC = _from_autos_to_perm(autsK)
    Gperm = _perm_to_gap_grp(permGC)
    dautsK1 = Dict{NfToNfMor, Int}()
    for w = 1:length(autsK1)
      dautsK1[autsK1[w]] = w
    end
    iso = GAP.Globals.IsomorphismGroups(Gperm, H)
    ElemGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(permGC))
    for w = 1:length(permGC)
      ElemGAP[w] = GAP.Globals.Image(iso, _perm_to_gap_perm(permGC[w]))
    end
    list[t].imgs_autos = ElemGAP
    #Restrict to the p-Sylow
    #Unfortunately, I need to compute the group structure.
    Gp = pSylow(autsK1, p)
    obstruction = false
    auts_for_conductor = Vector{Main.ForeignGAP.MPtr}()
    for g in autos
      #I create the cocycle

      function cocycle(aut1::NfToNfMor, aut2::NfToNfMor)
        s1 = restr[dautsK1[aut1]]
        a1 = GAP.Globals.Image(g, ElemGAP[s1])
        s2 = restr[dautsK1[aut2]]
        b1 = GAP.Globals.Image(g, ElemGAP[s2])
        rescoc = coc(a1, b1)
        return mod(rescoc, p)::Int
      end

      if cpa_issplit(K1, Gp, cocycle, p)
        obstruction = true
        push!(auts_for_conductor, g)
      end

    end
    #If I am here, all the algebras don't split. I return false
    if obstruction
      list[t].auts_for_conductors = auts_for_conductor
      list[t].proj_ext = mH
      list1[t] = true
    end
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  return list[findall(list1)]

end

###############################################################################
#
#  Cocycle computation: prime power case
#
###############################################################################

function _autos_to_check(G::Main.ForeignGAP.MPtr, K::Main.ForeignGAP.MPtr, E::Main.ForeignGAP.MPtr, mG::Main.ForeignGAP.MPtr)
  AutG = GAP.Globals.AutomorphismGroup(G)
  AutK = GAP.Globals.AutomorphismGroup(K)
  AutE = GAP.Globals.AutomorphismGroup(E)
  #I want to construct the map between the automorphism groups. The kernel is characteristic!
  gens = GAP.Globals.GeneratorsOfGroup(AutE)
  ind_auts_quo = Array{Main.ForeignGAP.MPtr, 1}(undef, length(gens))
  ind_auts_sub = Array{Main.ForeignGAP.MPtr, 1}(undef, length(gens))
  for s = 1:length(gens)
    ind_auts_quo[s] = GAP.Globals.InducedAutomorphism(mG, gens[s])
    ind_auts_sub[s] = GAP.Globals.RestrictedMapping(gens[s], K)
  end
  GProd = GAP.Globals.DirectProduct(AutG, AutK)
  EmbAutG = GAP.Globals.Embedding(GProd, 1)
  EmbAutK = GAP.Globals.Embedding(GProd, 2)
  gensubs = Vector{Main.ForeignGAP.MPtr}(undef, length(gens))
  for s = 1:length(gens)
    gensubs[s] = GAP.Globals.Image(EmbAutG, ind_auts_quo[s]) * GAP.Globals.Image(EmbAutK, ind_auts_sub[s])
  end
  S = GAP.Globals.Subgroup(GProd, GAP.julia_to_gap(gensubs))
  T = GAP.Globals.List(GAP.Globals.RightTransversal(GProd, S))
  #Now, I have to recover the automorphisms in AutG and AutK
  ProjAutG = GAP.Globals.Projection(GProd, 1)
  ProjAutK = GAP.Globals.Projection(GProd, 2)
  res = Vector{Tuple{Main.ForeignGAP.MPtr, Main.ForeignGAP.MPtr}}(undef, length(T))
  for i = 1:length(T)
    res[i] = (GAP.Globals.Image(ProjAutG, T[i]), GAP.Globals.Image(ProjAutK, T[i]))
  end
  return res
end

function cocycle_computation_pp(L::Main.ForeignGAP.MPtr, i::Int)

  proj = GAP.Globals.NaturalHomomorphismByNormalSubgroup(L[1], L[i+1])
  target_grp = GAP.Globals.ImagesSource(proj)
  mH1 = GAP.Globals.NaturalHomomorphismByNormalSubgroup(target_grp, GAP.Globals.Image(proj, L[i]))
  H1 = GAP.Globals.ImagesSource(mH1)
  K = GAP.Globals.Kernel(mH1)
  # The kernel is cyclic by hypothesis, but Gap represent it as a pc group with more than one generator.
  # Brute force at the moment.
  gen = GAP.Globals.MinimalGeneratingSet(K)[1]

  Elems = GAP.Globals.Elements(H1)
  MatCoc = Array{Int, 2}(undef, length(Elems), length(Elems))
  Preimags = Vector{Main.ForeignGAP.MPtr}(undef, length(Elems))
  for i = 1:length(Elems)
    Preimags[i] = GAP.Globals.PreImagesRepresentative(mH1, Elems[i])
  end
  for i = 1:length(Elems)
    x1 = Preimags[i]
    for j = 1:length(Elems)
      y1 = Preimags[j]
      xy = Elems[i] * Elems[j]
      k = 1
      while Elems[k] != xy
        k += 1
      end
      xy1 = Preimags[k]
      MatCoc[i, j] = _find_exp(gen, x1*y1*GAP.Globals.Inverse(xy1))
    end
  end
  local cocycle
  let Elems = Elems, MatCoc = MatCoc
    function cocycle(x::Main.ForeignGAP.MPtr, y::Main.ForeignGAP.MPtr)
      i = 1
      while Elems[i] != x
        i += 1
      end
      j = 1
      while Elems[j] != y
        j += 1
      end
      return MatCoc[i, j]
    end
  end

  autos = _autos_to_check(H1, K, target_grp, mH1)
  #For the automorphisms of K, since K is cyclic, I only need the image of
  # 1 under the automorphism.
  autos_new = Vector{Tuple{Main.ForeignGAP.MPtr, Int}}(undef, length(autos))
  for i = 1:length(autos)
    ex = _find_exp(gen, GAP.Globals.Image(autos[i][2], gen))
    autos_new[i] = (autos[i][1], ex)
  end

  action = Vector{Int}(undef, length(Elems))
  for i = 1:length(Elems)
    x1 = Preimags[i]
    y = x1 * gen * GAP.Globals.Inverse(x1)
    action[i] = _find_exp(gen, y)
  end

  return mH1, autos_new, cocycle, action, Elems

end

function _ext_cycl(G::Array{Hecke.NfToNfMor, 1}, d::Int)

  K = domain(G[1])
  Kc = cyclotomic_extension(K, d, compute_maximal_order = true, compute_LLL_basis = false, cached = false)
  automorphisms(Kc; gens = G, copy = false)
  return Kc
end

###############################################################################
#
#  Brauer obstruction: prime power case
#
###############################################################################

function action_on_roots(G::Vector{NfToNfMor}, zeta::nf_elem, pv::Int)
  act_on_roots = Array{Int, 1}(undef, length(G))
  #@assert isone(zeta^pv)
  p = 11
  K = domain(G[1])
  Qx = parent(K.pol)
  d = discriminant(K.pol)
  while iszero(mod(numerator(d), p)) || iszero(mod(pv, p))
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  polsG = gfp_poly[Rx(g.prim_img) for g in G]
  zetaP = Rx(zeta)
  units = Vector{gfp_poly}(undef, pv-1)
  units[1] = zetaP
  for i = 2:pv-1
    units[i] = mod(units[i-1]*zetaP, fmod)
  end
  for w = 1:length(G)
    act = Hecke.compose_mod(zetaP, polsG[w], fmod)
    s = 1
    while act != units[s]
      s += 1
    end
    act_on_roots[w] = s
    #@assert G[w](zeta) == zeta^s
  end
  return act_on_roots
end

function restriction(autsK1::Vector{NfToNfMor}, autsK::Vector{NfToNfMor}, mp::NfToNfMor)
  K = domain(mp)
  K1 = codomain(mp)
  p = 11
  d1 = numerator(discriminant(K.pol))
  d2 = numerator(discriminant(K1.pol))
  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(K1.pol)
  mp_pol = Rx(mp.prim_img)
  imgs = Vector{gfp_poly}(undef, length(autsK))
  for i = 1:length(imgs)
    imgs[i] = Hecke.compose_mod(Rx(autsK[i].prim_img), mp_pol, fmod)
  end
  res = Vector{Int}(undef, length(autsK1))
  for i = 1:length(res)
    img = Hecke.compose_mod(mp_pol, Rx(autsK1[i].prim_img), fmod)
    for j = 1:length(autsK)
      if img == imgs[j]
        res[i] = j
        break
      end
    end
    #@assert mp(autsK[res[i]].prim_img) == autsK1[i](mp.prim_img)
  end
  return res
end

function pSylow(G::Vector{NfToNfMor}, p::Int)
  com, uncom = ppio(length(G), p)
  if uncom == 1
    return G
  end
  permGC = _from_autos_to_perm(G)
  Gperm = _perm_to_gap_grp(permGC)
  PermGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(permGC))
  for w = 1:length(permGC)
    PermGAP[w] = _perm_to_gap_perm(permGC[w])
  end
  return pSylow(Gperm, PermGAP, G, p)
end

function check_Brauer_obstruction_pp(list::Vector{FieldsTower}, L::Main.ForeignGAP.MPtr, i::Int, p::Int, v::Int)
  list1 = falses(length(list))
  pv = p^v
  mH, autos, coc, action_grp, Elems = cocycle_computation_pp(L, i)
  for t = 1:length(list)
    @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Fields to test: $(length(list)-t+1)"
    Kc = _ext_cycl(list[t].generators_of_automorphisms, pv)
    K = list[t].field
    K1 = Kc.Ka
    autsK = automorphisms(K, copy = false)
    autsK1 = automorphisms(K1, copy = false)
    restr = restriction(autsK1, autsK, Kc.mp[2])
    #I construct the group and the isomorphisms between the automorphisms and the gap group.
    permGC = _from_autos_to_perm(autsK)
    Gperm = _perm_to_gap_grp(permGC)
    dautsK1 = Dict{NfToNfMor, Int}()
    for w = 1:length(autsK1)
      dautsK1[autsK1[w]] = w
    end
    H = GAP.Globals.ImagesSource(mH)
    iso = GAP.Globals.IsomorphismGroups(Gperm, H)
    ElemGAP = Vector{Main.ForeignGAP.MPtr}(undef, length(permGC))
    for w = 1:length(permGC)
      ElemGAP[w] = GAP.Globals.Image(iso, _perm_to_gap_perm(permGC[w]))
    end
    list[t].imgs_autos = ElemGAP
    #Restrict to the p-Sylow
    #Unfortunately, I need to compute the group structure.
    Gp = pSylow(autsK1, p)
    act_on_roots = action_on_roots(Gp, Kc.mp[1]\(gen(Kc.Kr)), pv)
    obstruction = false
    auts_for_conductors = Vector{Main.ForeignGAP.MPtr}()
    for (g1, g2) in autos
      #I create the cocycle
      local cocycle
      let restr = restr, ElemGAP = ElemGAP, dautsK1 = dautsK1, coc = coc, pv = pv
        function cocycle(aut1::NfToNfMor, aut2::NfToNfMor)
          s1 = restr[dautsK1[aut1]]
          a1 = GAP.Globals.PreImagesRepresentative(g1, ElemGAP[s1])
          s2 = restr[dautsK1[aut2]]
          b1 = GAP.Globals.PreImagesRepresentative(g1, ElemGAP[s2])
          rescoc = coc(a1, b1)*g2
          return mod(rescoc, pv)::Int
        end
      end
      #I have to find the subgroup of Gp such that the action of Gp on the roots of unity
      #coincides with the action
      Stab = NfToNfMor[]
      for w = 1:length(Gp)
        if Gp[w].prim_img == gen(K1)
          push!(Stab, Gp[w])
          continue
        end
        s1 = restr[dautsK1[Gp[w]]]
        img_el = GAP.Globals.PreImagesRepresentative(g1, ElemGAP[s1])
        ind = 1
        for s = 1:length(Elems)
          if img_el == Elems[s]
            ind = s
            break
          end
        end
        if action_grp[ind] == act_on_roots[w]
          push!(Stab, Gp[w])
        end
      end
      #@assert check_cocycle(Stab, cocycle, pv)

      #If the group acts as the roots of unity, we have a Brauer embedding problem, so
      # we only have to check one algebra.
      if length(Stab) == length(Gp)
        fl = cpa_issplit(K1, Gp, cocycle, pv)
        if fl
          obstruction = true
          push!(auts_for_conductors, g1)
          continue
        end
      end
      #Otherwise, I reduce the embedding problem to 2 different Brauer embedding problems.
      #One corresponds to the subextension of prime degree, the other one corresponds
      # the the embedding problem given by the subgroup of G having the same action on G and
      # on the roots of unity.
      fl = cpa_issplit(K1, Gp, cocycle, p) && cpa_issplit(K1, Gp, Stab, cocycle, p, v)

      if fl
        obstruction = true
        push!(auts_for_conductors, g1)
        continue
      end

    end
    #If I am here, all the algebras don't split. I return false
    if obstruction
      list[t].auts_for_conductors = auts_for_conductors
      list[t].proj_ext = mH
      list1[t] = true
    end
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  it = findall(list1)
  return FieldsTower[list[t] for t in it]

end

###############################################################################
#
#  Brauer obstruction: interface
#
###############################################################################

function check_Brauer_obstruction(list::Vector{FieldsTower}, L::Main.ForeignGAP.MPtr, i::Int, l::Int)
  d = degree(list[1])
  if gcd(d, l) == 1
    return list
  end
  fac = factor(l)
  for (p, v) in fac
    if d % Int(p) != 0
      continue
    end
    if p == 2
      list = _Brauer_at_two(list, L, i)
      if v > 1
        @vprint :Fields 1 "Fields after first step: $(length(list))\n"
        list = check_Brauer_obstruction_pp(list, L, i, 2, v)
      end
    else
      if v > 1
        list = check_Brauer_obstruction_pp(list, L, i, Int(p), v)
      else
        list = _Brauer_prime_case(list, L, i, Int(p))
      end
    end
  end
  return list
end

###############################################################################
#
#  IsSplit for crossed product algebras: the action of G on the extension
#  correspond to the action on the roots of unity
#
###############################################################################

#p must be a prime power
function cpa_issplit(K::AnticNumberField, G::Vector{NfToNfMor}, Coc::Function, pv::Int)

  p = ispower(pv)[2]
  @vtime :BrauerObst 1 if p == 2 && !istotally_real(K) && !is_split_at_infinity(K, Coc)
    return false
  end
  # Now, the finite primes.
  # The crossed product algebra is ramified at the ramified primes of the field (and at the primes dividing the values
  # of the cocycle, but our cocycle has values in the roots of unity...
  # I only need to check the tame ramification.
  # The exact sequence on Brauer groups and completion tells me that I have one degree of freedom! :)
  O = maximal_order(K)
  lp = Hecke.ramified_primes(O)
  if p in lp
    for q in lp
      if q != p
        @vtime :BrauerObst 1 fl = is_split_at_p(O, G, Coc, Int(q), pv)
        if !fl
          return false
        end
      end
    end
  else
    for i = 1:length(lp)-1
      @vtime :BrauerObst 1 fl = is_split_at_p(O, G, Coc, Int(lp[i]), pv)
      if !fl
        return false
      end
    end
  end
  return true

end

function is_split_at_infinity(K::AnticNumberField, Coc::Function)

  aut = complex_conjugation(K)
  if Coc(aut, aut) == 1
    @vprint :Fields 3 "Real places!\n"
    return false
  else
    return true
  end

end

function inertia_subgroup(F, mF, G::Array{NfToNfMor, 1})

  @assert !isempty(G)
  K = domain(G[1])
  O = maximal_order(K)
  gF = gen(F)
  igF = K(mF\gF)
  inertia_grp = NfToNfMor[]
  for g in G
    if g.prim_img == gen(K)
      push!(inertia_grp, g)
      continue
    end
    y = mF(O(g(igF), false))
    if y == gF
      push!(inertia_grp, g)
    end
  end
  return inertia_grp

end

function decomposition_group(G::Array{NfToNfMor, 1}, P::NfOrdIdl, orderG::Int = length(G))

  @assert !isempty(G)
  O = order(P)
  K = nf(O)
  dec_group = NfToNfMor[]
  F, mF = Hecke.ResidueFieldSmall(O, P)
  for g in G
    if g.prim_img == gen(K)
      push!(dec_group, g)
      continue
    end
    y = O(g(K(P.gen_two)), false)
    if iszero(image(mF, y))
      push!(dec_group, g)
    end
    if length(dec_group) == orderG
      break
    end
  end
  return dec_group

end

function _find_theta(G::Vector{NfToNfMor}, F::FqNmodFiniteField, mF::Hecke.NfOrdToFqNmodMor, e::Int)
  K = domain(G[1])
  O = maximal_order(K)
  p = ispower(e)[2]
  t = div(e, p)
  gF = gen(F)
  igF = K(mF\gF)
  gK = gen(K)
  theta = G[1]
  for i = 1:length(G)
    theta = G[i]
    #I make sure that the element is a generator of the inertia subgroup
    if mF(O(theta(igF), false)) == gF && (theta^t).prim_img != gK
      break
    end
  end
  return theta
end

function _find_frob(G::Vector{NfToNfMor}, F::FqNmodFiniteField, mF::Hecke.NfOrdToFqNmodMor, e::Int, f::Int, q::Int, theta::NfToNfMor)
  K = domain(G[1])
  O = maximal_order(K)
  gK = gen(K)
  p = ispower(e)[2]
  t = div(e, p)
  frob = G[1]
  expo = mod(q, e)
  gF = gen(F)
  igF = K(mF\gF)
  rhs = gF^q
  for i = 1:length(G)
    frob = G[i]
    if frob == theta
      continue
    end
    #I make sure that g1 is a Frobenius
    if mF(O(frob(igF), false)) != rhs
      continue
    end
    #Now, I check the relation
    gc = frob * theta
    gq = (theta^expo) * frob
    if gc == gq
      break
    end
  end
  return frob
end

#See Gerald J. Janusz (1980) Crossed product orders and the schur index,
#Communications in Algebra, 8:7, 697-706

function is_split_at_p(O::NfOrd, GC::Array{NfToNfMor, 1}, Coc::Function, p::Int, n::Int)

  lp = prime_decomposition(O, p, cached = true)
  e = gcd(length(GC), lp[1][2])
  if e == 1
    return true
  end
  P = lp[1][1]
  f = gcd(length(GC), degree(P))
  if f == 1 && iszero(mod(norm(P)-1, e*n))
    return true
  end
  if length(lp) != 1
    @vtime :BrauerObst 1  Gp = decomposition_group(GC, P, e*f)
    #I don't really want the decomposition group of the p-sylow, but the p-sylow of the decomposition group.
    #Therefore, if the p-sylow is not normal, I have to try different primes.
    i = 1
    while length(Gp) != e*f
      i += 1
      P = lp[i][1]
      Gp = decomposition_group(GC, P, e*f)
    end
  else
    Gp = GC
  end
  f1 = divexact(degree(P), f)
  q = p^f1 #Cardinality of the residue field

  F, mF = Hecke.ResidueFieldSmall(O, P)
  theta = _find_theta(Gp, F, mF, e)

  #I have found my generators. Now, we find the elements to test if it splits.
  @assert divisible(norm(P)-1, e)
  c = divexact(norm(P)-1, e)
  if !iszero(mod(c, n))
    powtheta = theta
    zeta = 0
    for k = 1:e-1
      zeta += Coc(powtheta, theta)
      powtheta *= theta
    end
    zeta = mod(zeta * c, n)
    if !iszero(zeta)
      return false
    end
  end

  #The second element is a little more complicated.
  if f == 1
    return true
  end
  frob = _find_frob(Gp, F, mF, e, f, q, theta)
  if iszero(mod(q-1, e*n))
    lambda = mod(Coc(frob, theta)- Coc(theta, frob), n)
    return iszero(lambda)
  end
  lambda = Coc(frob, theta)
  powtheta = theta
  s, t = divrem(q-1, e)
  if !iszero(mod(s+1, n))
    for k = 1:t
      lambda -= (s+1) * Coc(powtheta, theta)
      powtheta *= theta
    end
  else
    powtheta = theta^(t+1)
  end
  if !iszero(mod(s, n))
    for k = t+1:(e-1)
      lambda -= s * Coc(powtheta, theta)
      powtheta *= theta
    end
  end
  powtheta = theta^mod(q, e)
  lambda = mod(lambda - Coc(powtheta, frob), n)
  return iszero(lambda)

end

################################################################################
#
#  IsSplit for crossed product algebras: prime power case
#
################################################################################

function cpa_issplit(K::AnticNumberField, G::Vector{NfToNfMor}, Stab::Vector{NfToNfMor}, Coc::Function, p::Int, v::Int)

  @vtime :BrauerObst 1 if p == 2 && !istotally_real(K) && !is_split_at_infinity(K, Coc)
    return false
  end
  # Now, the finite primes.
  # The crossed product algebra is ramified at the ramified primes of the field (and at the primes dividing the values
  # of the cocycle, but our cocycle has values in the roots of unity...
  # I only need to check the tame ramification.
  # The exact sequence on Brauer groups and completion tells me that I have one degree of freedom! :)
  O = maximal_order(K)
  lp = Hecke.ramified_primes(O)
  if p in lp
    for q in lp
      if q != p
        @vtime :BrauerObst 1 fl = is_split_at_p(O, G, Stab, Coc, Int(q), p^v)
        if !fl
          return false
        end
      end
    end
  else
    for i = 1:length(lp)-1
      @vtime :BrauerObst 1 fl = is_split_at_p(O, G, Stab, Coc, Int(lp[i]), p^v)
      if !fl
        return false
      end
    end
  end
  return true

end

function is_split_at_p(O::NfOrd, GC::Vector{NfToNfMor}, Stab::Vector{NfToNfMor}, Coc::Function, p::Int, n::Int)

  lp = prime_decomposition(O, p, cached = true)
  e = gcd(length(GC), lp[1][2])
  if e == 1
    return true
  end

  P = lp[1][1]
  @assert divisible(norm(P)-1, e)
  c = divexact(norm(P)-1, e)
  f = gcd(length(GC), degree(P))
  if f == 1 && iszero(mod(c, n))
    return true
  end
  if length(lp) != 1
    @vtime :BrauerObst 1  Gp = decomposition_group(GC, P, e*f)
    #I don't really want the decomposition group of the p-sylow, but the p-sylow of the decomposition group.
    #Therefore, if the p-sylow is not normal, I have to try different primes.
    i = 1
    while length(Gp) != e*f
      i += 1
      P = lp[i][1]
      Gp = decomposition_group(GC, P, e*f)
    end
  else
    Gp = GC
  end
  GpStab = NfToNfMor[]
  for i = 1:length(Gp)
    if Gp[i] in Stab
      push!(GpStab, Gp[i])
    end
  end
  #I need the inertia subgroup to determine e and f
  F, mF = Hecke.ResidueFieldSmall(O, P)
  InGrp = inertia_subgroup(F, mF, GpStab)
  e = length(InGrp)
  if e == 1
    return true
  end
  f = divexact(length(GpStab), e)
  c = divexact(norm(P)-1, e)
  if f == 1 && iszero(mod(c, n))
    return true
  end
  f1 = divexact(degree(P), f)
  q = p^f1 #Cardinality of the residue field

  theta = _find_theta(GpStab, F, mF, e)
  #I have found my generators. Now, we find the elements to test if it splits.
  if !iszero(mod(c, n))
    powtheta = theta
    zeta = 0
    for k = 1:e-1
      zeta += Coc(powtheta, theta)
      powtheta *= theta
    end
    zeta = mod(zeta * c, n)
    if !iszero(zeta)
      return false
    end
  end

  if f == 1
    return true
  end
  frob = _find_frob(GpStab, F, mF, e, f, q, theta)
  if iszero(mod(q-1, e*n))
    lambda = mod(Coc(frob, theta)- Coc(theta, frob), n)
    return iszero(lambda)
  end
  lambda = Coc(frob, theta)
  powtheta = theta
  s, t = divrem(q-1, e)
  if !iszero(mod(s+1, n))
    for k = 1:t
      lambda -= (s+1) * Coc(powtheta, theta)
      powtheta *= theta
    end
  else
    powtheta = theta^(t+1)
  end
  lambda -= Coc(powtheta, frob)
  if !iszero(mod(s, n))
    for k = t+1:(e-1)
      lambda -= s * Coc(powtheta, theta)
      powtheta *= theta
    end
  end
  powtheta = theta^mod(q, e)
  lambda = mod(lambda - Coc(powtheta, frob), n)
  return iszero(lambda)

end


text:
using KernelAbstractions: @index, @kernel
using KernelAbstractions.Extras.LoopInfo: @unroll
using Oceananigans.Architectures: device_event, architecture
using Oceananigans.Utils: launch!
using Oceananigans.Grids

@kernel function _compute_mixed_layer_depth!(h, grid, b, Δb)
    i, j = @index(Global, NTuple)

    # Use this to set mld to surface
    z_surface = znode(Center(), Center(), Face(), i, j, grid.Nz, grid)
    b_surface = @inbounds b[i, j, grid.Nz] # buoyancy at surface

    # at the first iteration, mld is at the surface
    mld_ij = z_surface
    searching_mld = true

    @unroll for k in grid.Nz -1 : -1 : 1 # scroll from point just below surface
        # buoyancy and buoyancy decrease
        b_ijk = @inbounds b[i, j, k]
        Δb_ijk = b_surface-b_ijk

        # if below mixed layer and we havent yet found it
        if (Δb_ijk>Δb)&(searching_mld)
            # buoyancy and buoyancy decrease for the cell above
            b_ijk_above = @inbounds b[i, j, k+1]
            Δb_ijk_above = b_surface-b_ijk_above

            # height for both cells
            z_face = znode(Center(), Center(), Face(), i, j, k, grid)
            z_face_above = znode(Center(), Center(), Face(), i, j, k+1, grid)

            # get Δz and Δb
            Δz_ijk = z_face_above-z_face

            # interpolate linearly to obtain the mld
            mld_ij = max(-((Δz_ijk/(Δb_ijk_above-Δb_ijk))*(Δb-Δb_ijk)+z_face), 0)

            # stop searching (skip this calculation for other levels)
            searching_mld=false
        end
    end

    # Note "-" since `h` is supposed to be "depth" rather than "height"
    @inbounds h[i, j, 1] = mld_ij

end

"""
    MixedLayerDepth!(h, b, Δb)

Compute and update the mixed layer depth field (h) from the buoyancy field (b) using the buoyancy decrease criterium (Δb).
"""
function MixedLayerDepth!(h, b, Δb)
    grid = h.grid
    arch = architecture(grid)

    event = launch!(arch, grid, :xy,
                    _compute_mixed_layer_depth!, h, grid, b, Δb,
                    dependencies = device_event(arch))

    wait(device_event(arch), event)

    return nothing
end

ext: jl

max_stars_repo_name: iuryt/biocenanigans

max_issues_repo_name: iuryt/biocenanigans


text:
import Base.==

"""
Determine if `p, q` is a handle for `braid0`.
"""
function dehornoyishandle(braid0::BraidElement, p::Int, q::Int)::Bool
    braid = simplifybraid(braid0)
    if (braid.crossings[p][2] + braid.crossings[q][2] == 0) && (braid.crossings[p][1] == braid.crossings[q][1])
        v = braid.crossings[(p+1):(q-1)]
        j = braid.crossings[p][1]
        jp1 = filter(x -> x[1] == (j+1), v)
        jp1pos = filter(x -> x[2] > 0, jp1)
        jp1neg = filter(x -> x[2] < 0, jp1)
        if length(jp1pos) > 0 && length(jp1neg) > 0
            return false
        end
        for k in v
            if !(k[1] < (j-1) || k[1] > j) # Original
            #if !(k[1] < (j-1) || k[1] > (j+1)) # Ambiguity in the paper
                return false
            end
        end
        true
    else
        false
    end
end

"""
Find the first handle in `braid0` if one exists.
"""
function dehornoyfirsthandle(braid0::BraidElement)::Union{Tuple{Int,Int},Void}
    braid = simplifybraid(braid0)
    #for q = length(braid.crossings):-1:2
    for q = 2:length(braid.crossings)
    #for p = 1:(length(braid.crossings) - 1)
        #for q = (p+1):length(braid.crossings)
        #for p = 1:(q-1)
        for p = (q-1):-1:1
            if dehornoyishandle(braid0, p, q)
                return (p, q)
            end
        end
    end
    return nothing
end

"""
Reduce according to the alphabetic homomorphism of Dehornoy.
"""
function dehornoyalphareduce(braid0::BraidElement, handlep::Int, handleq::Int)::BraidElement
    nh = []
    j, e = braid0.crossings[handlep]
    for s in braid0.crossings[handlep:handleq]
        idx, ex = s
        if (idx != j) && (idx != (j+1))
            push!(nh, s)
        elseif idx == (j+1)
            push!(nh, (idx, -1*e))
            push!(nh, (j, ex))
            push!(nh, (idx, e))
        end
    end
    nw = cat(1, braid0.crossings[1:handlep-1], nh, braid0.crossings[handleq+1:end])
    b = BraidElement(braid0.group, nw)
    simplifybraid(b)
end

"""
Reduce according to the algorithm of Dehornoy.
"""
function dehornoyreduce(braid0::BraidElement)::BraidElement
    braid = unprettybraid(simplifybraid(braid0))
    while true
        handle = dehornoyfirsthandle(braid)
        if handle == nothing
            return prettybraid(simplifybraid(braid))
        else
            p, q = handle
            braid = dehornoyalphareduce(braid, p, q)
        end
    end
end

"""
Is a braid empty?
"""
function isempty(b::BraidElement)::Bool
    length(b.crossings) == 0
end

"""
Are two braids equal (in terms of ambient isotopy)?
"""
function ==(b1::BraidElement, b2::BraidElement)
    isempty(dehornoyreduce(b1 * inv(b2)))
end


text:
using ModelTES, DifferentialEquations

tes = ModelTES.highEpix()

@time out = rk8(12000,1e-7, tes, 1000, 2000);

@time out2 = ModelTES.pulse(12000,1e-7, tes, 1000, 2000);

nsample=12000
dt = 1e-7
bt = tes
E=1000
npresamples = 2000
p = bt.p
# u = [T,I]
u0 = [bt.T0+E*ModelTES.J_per_eV/p.C, bt.I0]
function du(t,u)
    T,I = u[1],u[2]
    r = ModelTES.R(I,T,p)
    [ModelTES.dT(I, T, p.k, p.n, p.Tbath, p.C, r),
             ModelTES.dI(I,T, bt.V, p.Rl, p.L, r)]
end
prob = ODEProblem(tes, u0, (0.0,nsample*dt))
@time sol = DifferentialEquations.solve(prob,DifferentialEquations.Vern8(),dt=1e-9,abstol=1e-11,reltol=1e-11);
@time out =  sol(0:dt:(nsample-1)*dt);
@time sol2 = DifferentialEquations.solve(prob,DifferentialEquations.Vern8(),dt=1e-9,abstol=1e-11,reltol=1e-11, saveat=0:dt:(nsample-1)*dt, dense=false);

T=[o[1] for o in out];
I=[o[2] for o in out];
ModelTES.TESRecord(T,I, ModelTES.R(I,T,bt.p),dt)

function lorenz(t,u,du)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
lu0 = [1.0;0.0;0.0]
tspan = (0.0,1.0)
dt = 0.0001
times = first(tspan):dt:last(tspan)
method = DifferentialEquations.Tsit5()
lprob = ODEProblem(lorenz,lu0,tspan)
lt1 = @elapsed lsol = solve(lprob, method);
lt2 = @elapsed lsol(times);
lt3 = @elapsed lsol2 = solve(lprob, method, saveat=times, dense=false, save_timeseries=false);
lsol3 = solve(lprob, method, saveat=collect(times), dense=false, save_timeseries=false);
@show length(lsol2.t)
@show length(lsol3.t)

dt=1e-5
nsample=1000
prob = ODEProblem(tes, u0, (0.0,nsample*dt))
times2 = 0:dt:(nsample-1)*dt
t1 = @elapsed sol = solve(prob,method,dt=1e-9,reltol=1e-9);
t2=@elapsed out =  sol(0:dt:(nsample-1)*dt);
t3=@elapsed sol2 = solve(prob,method,dt=1e-9,reltol=1e-9, saveat=times2, dense=false);
t3=@elapsed sol2_abstol = solve(prob,method,dt=1e-9,abstol=1e-9,reltol=1e-9, saveat=times2, dense=false);

@printf("lorenz %s, solve %0.6f, sol interp %0.6f, solve saveat %0.6f\n", repr(typeof(method)), lt1, lt2, lt3);@printf("%s, solve %0.6f, sol interp %0.6f, solve saveat %0.6f", repr(typeof(method)), t1, t2, t3)

@time sol_saveat = solve(prob,method,dt=1e-9,reltol=1e-9, dense=false, saveat=times2, save_timeseries=false);

using PyPlot
figure()
plot(sol2_abstol.t,sol2_abstol[:,2],label="abstol,reltol=1e-9")
plot(sol2.t,sol2[:,2],label="reltol=1e-9")
xlabel("time (s)")
ylabel("current (A)")


text:
using CxxWrap

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, package EigenCpp did not build properly")
end
include(depsfile)

wrap_modules(_l_eigen_wrap)


text:
#Expo Atmoshpere parameters
 ρ0expo(::Type{Earth}) = 1.225 #kg/m^3
 Rexpo(::Type{Earth}) = 287 #J/Kg K For Earth Air
 γexpo(::Type{Earth}) = 1.4 #For Earth Air
 g0expo(::Type{Earth}) = mu(Earth)/equatorial_radius(Earth)^2*1000
 H_sexpo(::Type{Earth}) = 7125 #m
 Texpo(::Type{Earth}) = H_s*g0/R
 aexpo(::Type{Earth}) = sqrt(γ*R*T)

function ExpoAtmos(Body::Type{<:CelestialBody}, PosRsph::Vector{Float64})
  h = Altitude(Body, PosRsph)
  ρ = ρ0expo(Body)*exp(-h/H_sexpo(Body))
  p = ρ*Rexpo(Body)*Texpo(Body)
  return [ρ, p, Texpo(Body), aexpo(Body)]
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct EXAC1A <: AVR
        Tr::Float64
        Tb::Float64
        Tc::Float64
        Ka::Float64
        Ta::Float64
        Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Te::Float64
        Kf::Float64
        Tf::Float64
        Kc::Float64
        Kd::Float64
        Ke::Float64
        E_sat::Tuple{Float64, Float64}
        Se::Tuple{Float64, Float64}
        Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        V_ref::Float64
        saturation_coeffs::Tuple{Float64, Float64}
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes.StateType}
        internal::InfrastructureSystemsInternal
    end

Modified ESAC1A. This excitation systems consists of an alternator main exciter feeding its output via non-controlled rectifiers.
The exciter does not employ self-excitation, and the voltage regulator power is taken from a source that is not affected by external transients.
Parameters of IEEE Std 421.5 Type AC1A Excitacion System. EXAC1A in PSSE and PSLF

# Arguments
- `Tr::Float64`: Regulator input filter time constant in s, validation range: `(0, 0.5)`, action if invalid: `warn`
- `Tb::Float64`: Regulator denominator (lag) time constant in s, validation range: `(0, 20)`, action if invalid: `warn`
- `Tc::Float64`: Regulator numerator (lead) time constant in s, validation range: `(0, 20)`, action if invalid: `warn`
- `Ka::Float64`: Regulator output gain, validation range: `(0, 1000)`, action if invalid: `warn`
- `Ta::Float64`: Regulator output time constant in s, validation range: `(0, 10)`
- `Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for regulator output `(Va_min, Va_max)`
- `Te::Float64`: Exciter field time constant in s, validation range: `(eps(), 2)`, action if invalid: `error`
- `Kf::Float64`: Rate feedback excitation system stabilizer gain, validation range: `(0, 0.3)`, action if invalid: `warn`
- `Tf::Float64`: Rate feedback time constant, validation range: `(eps(), 1.5)`, action if invalid: `error`
- `Kc::Float64`: Rectifier loading factor proportional to commutating reactance, validation range: `(0, 1)`, action if invalid: `warn`
- `Kd::Float64`: Demagnetizing factor, function of exciter alternator reactances, validation range: `(0, 1)`, action if invalid: `warn`
- `Ke::Float64`: Exciter field proportional constant, validation range: `(0, 1)`
- `E_sat::Tuple{Float64, Float64}`: Exciter output voltage for saturation factor: (E1, E2)
- `Se::Tuple{Float64, Float64}`: Exciter saturation factor at exciter output voltage: (Se(E1), Se(E2))
- `Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limits for exciter field voltage: `(Vr_min, Vr_max)`
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `saturation_coeffs::Tuple{Float64, Float64}`: Coefficients (A,B) of the function: Se(x) = B(x - A)^2/x
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vm: Sensed terminal voltage,
	Vr1: Lead-lag state,
	Vr2: Regulator output state,
	Ve: Integrator output state,
	Vr3: Feedback output state
- `n_states::Int`: EXAC1A has 5 states
- `states_types::Vector{StateTypes.StateType}`: EXAC1A has 5 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct EXAC1A <: AVR
    "Regulator input filter time constant in s"
    Tr::Float64
    "Regulator denominator (lag) time constant in s"
    Tb::Float64
    "Regulator numerator (lead) time constant in s"
    Tc::Float64
    "Regulator output gain"
    Ka::Float64
    "Regulator output time constant in s"
    Ta::Float64
    "Limits for regulator output `(Va_min, Va_max)`"
    Va_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Exciter field time constant in s"
    Te::Float64
    "Rate feedback excitation system stabilizer gain"
    Kf::Float64
    "Rate feedback time constant"
    Tf::Float64
    "Rectifier loading factor proportional to commutating reactance"
    Kc::Float64
    "Demagnetizing factor, function of exciter alternator reactances"
    Kd::Float64
    "Exciter field proportional constant"
    Ke::Float64
    "Exciter output voltage for saturation factor: (E1, E2)"
    E_sat::Tuple{Float64, Float64}
    "Exciter saturation factor at exciter output voltage: (Se(E1), Se(E2))"
    Se::Tuple{Float64, Float64}
    "Limits for exciter field voltage: `(Vr_min, Vr_max)`"
    Vr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Reference Voltage Set-point"
    V_ref::Float64
    "Coefficients (A,B) of the function: Se(x) = B(x - A)^2/x"
    saturation_coeffs::Tuple{Float64, Float64}
    ext::Dict{String, Any}
    "The states are:
	Vm: Sensed terminal voltage,
	Vr1: Lead-lag state,
	Vr2: Regulator output state,
	Ve: Integrator output state,
	Vr3: Feedback output state"
    states::Vector{Symbol}
    "EXAC1A has 5 states"
    n_states::Int
    "EXAC1A has 5 states"
    states_types::Vector{StateTypes.StateType}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function EXAC1A(Tr, Tb, Tc, Ka, Ta, Va_lim, Te, Kf, Tf, Kc, Kd, Ke, E_sat, Se, Vr_lim, V_ref=1.0, saturation_coeffs=PowerSystems.get_avr_saturation(E_sat, Se), ext=Dict{String, Any}(), )
    EXAC1A(Tr, Tb, Tc, Ka, Ta, Va_lim, Te, Kf, Tf, Kc, Kd, Ke, E_sat, Se, Vr_lim, V_ref, saturation_coeffs, ext, [:Vm, :Vr1, :Vr2, :Ve, :Vr3], 5, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function EXAC1A(; Tr, Tb, Tc, Ka, Ta, Va_lim, Te, Kf, Tf, Kc, Kd, Ke, E_sat, Se, Vr_lim, V_ref=1.0, saturation_coeffs=PowerSystems.get_avr_saturation(E_sat, Se), ext=Dict{String, Any}(), states=[:Vm, :Vr1, :Vr2, :Ve, :Vr3], n_states=5, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    EXAC1A(Tr, Tb, Tc, Ka, Ta, Va_lim, Te, Kf, Tf, Kc, Kd, Ke, E_sat, Se, Vr_lim, V_ref, saturation_coeffs, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function EXAC1A(::Nothing)
    EXAC1A(;
        Tr=0,
        Tb=0,
        Tc=0,
        Ka=0,
        Ta=0,
        Va_lim=(min=0.0, max=0.0),
        Te=0,
        Kf=0,
        Tf=0,
        Kc=0,
        Kd=0,
        Ke=0,
        E_sat=(0.0, 0.0),
        Se=(0.0, 0.0),
        Vr_lim=(min=0.0, max=0.0),
        V_ref=0,
        saturation_coeffs=(0.0, 0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get [`EXAC1A`](@ref) `Tr`."""
get_Tr(value::EXAC1A) = value.Tr
"""Get [`EXAC1A`](@ref) `Tb`."""
get_Tb(value::EXAC1A) = value.Tb
"""Get [`EXAC1A`](@ref) `Tc`."""
get_Tc(value::EXAC1A) = value.Tc
"""Get [`EXAC1A`](@ref) `Ka`."""
get_Ka(value::EXAC1A) = value.Ka
"""Get [`EXAC1A`](@ref) `Ta`."""
get_Ta(value::EXAC1A) = value.Ta
"""Get [`EXAC1A`](@ref) `Va_lim`."""
get_Va_lim(value::EXAC1A) = value.Va_lim
"""Get [`EXAC1A`](@ref) `Te`."""
get_Te(value::EXAC1A) = value.Te
"""Get [`EXAC1A`](@ref) `Kf`."""
get_Kf(value::EXAC1A) = value.Kf
"""Get [`EXAC1A`](@ref) `Tf`."""
get_Tf(value::EXAC1A) = value.Tf
"""Get [`EXAC1A`](@ref) `Kc`."""
get_Kc(value::EXAC1A) = value.Kc
"""Get [`EXAC1A`](@ref) `Kd`."""
get_Kd(value::EXAC1A) = value.Kd
"""Get [`EXAC1A`](@ref) `Ke`."""
get_Ke(value::EXAC1A) = value.Ke
"""Get [`EXAC1A`](@ref) `E_sat`."""
get_E_sat(value::EXAC1A) = value.E_sat
"""Get [`EXAC1A`](@ref) `Se`."""
get_Se(value::EXAC1A) = value.Se
"""Get [`EXAC1A`](@ref) `Vr_lim`."""
get_Vr_lim(value::EXAC1A) = value.Vr_lim
"""Get [`EXAC1A`](@ref) `V_ref`."""
get_V_ref(value::EXAC1A) = value.V_ref
"""Get [`EXAC1A`](@ref) `saturation_coeffs`."""
get_saturation_coeffs(value::EXAC1A) = value.saturation_coeffs
"""Get [`EXAC1A`](@ref) `ext`."""
get_ext(value::EXAC1A) = value.ext
"""Get [`EXAC1A`](@ref) `states`."""
get_states(value::EXAC1A) = value.states
"""Get [`EXAC1A`](@ref) `n_states`."""
get_n_states(value::EXAC1A) = value.n_states
"""Get [`EXAC1A`](@ref) `states_types`."""
get_states_types(value::EXAC1A) = value.states_types
"""Get [`EXAC1A`](@ref) `internal`."""
get_internal(value::EXAC1A) = value.internal

"""Set [`EXAC1A`](@ref) `Tr`."""
set_Tr!(value::EXAC1A, val) = value.Tr = val
"""Set [`EXAC1A`](@ref) `Tb`."""
set_Tb!(value::EXAC1A, val) = value.Tb = val
"""Set [`EXAC1A`](@ref) `Tc`."""
set_Tc!(value::EXAC1A, val) = value.Tc = val
"""Set [`EXAC1A`](@ref) `Ka`."""
set_Ka!(value::EXAC1A, val) = value.Ka = val
"""Set [`EXAC1A`](@ref) `Ta`."""
set_Ta!(value::EXAC1A, val) = value.Ta = val
"""Set [`EXAC1A`](@ref) `Va_lim`."""
set_Va_lim!(value::EXAC1A, val) = value.Va_lim = val
"""Set [`EXAC1A`](@ref) `Te`."""
set_Te!(value::EXAC1A, val) = value.Te = val
"""Set [`EXAC1A`](@ref) `Kf`."""
set_Kf!(value::EXAC1A, val) = value.Kf = val
"""Set [`EXAC1A`](@ref) `Tf`."""
set_Tf!(value::EXAC1A, val) = value.Tf = val
"""Set [`EXAC1A`](@ref) `Kc`."""
set_Kc!(value::EXAC1A, val) = value.Kc = val
"""Set [`EXAC1A`](@ref) `Kd`."""
set_Kd!(value::EXAC1A, val) = value.Kd = val
"""Set [`EXAC1A`](@ref) `Ke`."""
set_Ke!(value::EXAC1A, val) = value.Ke = val
"""Set [`EXAC1A`](@ref) `E_sat`."""
set_E_sat!(value::EXAC1A, val) = value.E_sat = val
"""Set [`EXAC1A`](@ref) `Se`."""
set_Se!(value::EXAC1A, val) = value.Se = val
"""Set [`EXAC1A`](@ref) `Vr_lim`."""
set_Vr_lim!(value::EXAC1A, val) = value.Vr_lim = val
"""Set [`EXAC1A`](@ref) `V_ref`."""
set_V_ref!(value::EXAC1A, val) = value.V_ref = val
"""Set [`EXAC1A`](@ref) `saturation_coeffs`."""
set_saturation_coeffs!(value::EXAC1A, val) = value.saturation_coeffs = val
"""Set [`EXAC1A`](@ref) `ext`."""
set_ext!(value::EXAC1A, val) = value.ext = val
"""Set [`EXAC1A`](@ref) `states_types`."""
set_states_types!(value::EXAC1A, val) = value.states_types = val


text:
function LCG(model; args...)
    Linesearch(model, CGlin, profil=true)
end


text:
#File treating all the prediction functions

const pred_nodes,pred_weights = gausshermite(100) |> x->(x[1].*sqrt2,x[2]./sqrtπ)

"""
Compute the mean of the predicted latent distribution of `f` on `X_test` for the variational GP `model`

Return also the diagonal variance if `covf=true` and the full covariance if `fullcov=true`
"""
function _predict_f(model::AbstractGP{T},X_test::AbstractMatrix{<:Real};covf::Bool=true,fullcov::Bool=false) where {T}
    k_star = get_σ_k(model).*kernelmatrix.(get_kernel(model),[X_test],get_Z(model),obsdim=1)
    μf = k_star.*(get_K(model).\get_μ(model))
    if !covf
        return (μf,)
    end
    A = get_K(model).\([I].-get_K(model).\get_Σ(model))
    if fullcov
        k_starstar = get_σ_k(model).*(kernelmatrix.(get_kernel(model),[X_test],obsdim=1).+T(jitter)*[I])
        Σf = Symmetric.(k_starstar .- k_star.*A.*transpose.(k_star))
        return μf,Σf
    else
        k_starstar = get_σ_k(model).*(kerneldiagmatrix.(get_kernel(model),[X_test],obsdim=1).+[T(jitter)*ones(T,size(X_test,1))])
        σ²f = k_starstar .- opt_diag.(k_star.*A,k_star)
        return μf,σ²f
    end
end

function _predict_f(model::GP{T},X_test::AbstractMatrix{<:Real};covf::Bool=true,fullcov::Bool=false) where {T}
    k_star = get_σ_k(model).*kernelmatrix.(get_kernel(model),[X_test],get_Z(model),obsdim=1)
    μf = k_star.*mean_f(model)
    if !covf
        return (μf,)
    end
    A = [inv(model.f[1].K+model.likelihood.σ²*I)]
    if fullcov
        k_starstar = get_σ_k(model).*(kernelmatrix.(get_kernel(model),[X_test],obsdim=1).+T(jitter)*[I])
        Σf = Symmetric.(k_starstar .- k_star.*A.*transpose.(k_star))
        return μf,Σf
    else
        k_starstar = get_σ_k(model).*(kerneldiagmatrix.(get_kernel(model),[X_test],obsdim=1).+[T(jitter)*ones(T,size(X_test,1))])
        σ²f = k_starstar .- opt_diag.(k_star.*A,k_star)
        return μf,σ²f
    end
end

function _predict_f(model::MOSVGP{T},X_test::AbstractMatrix{<:Real};covf::Bool=true,fullcov::Bool=false) where {T}
    k_star = get_σ_k(model).*kernelmatrix.(get_kernel(model),[X_test],get_Z(model),obsdim=1)
    μf = k_star.*(get_K(model).\get_μ(model))
    μf = [[sum(vec(model.A[i,j,:]).*μf) for j in 1:model.nf_per_task[i]] for i in 1:model.nTask]
    if !covf
        return (μf,)
    end
    A = get_K(model).\([I].-get_K(model).\get_Σ(model))
    if fullcov
        k_starstar = get_σ_k(model).*(kernelmatrix.(get_kernel(model),[X_test],obsdim=1).+T(jitter)*[I])
        Σf = k_starstar .-  k_star.*A.*transpose.(k_star)
        Σf = [[sum(vec(model.A[i,j,:]).^2 .*Σf) for j in 1:model.nf_per_task[i]] for i in 1:model.nTask]
        return μf,Σf
    else
        k_starstar = get_σ_k(model).*(kerneldiagmatrix.(get_kernel(model),[X_test],obsdim=1).+[T(jitter)*ones(T,size(X_test,1))])
        σ²f = k_starstar .- opt_diag.(k_star.*A,k_star)
        σ²f = [[sum(vec(model.A[i,j,:]).^2 .*σ²f) for j in 1:model.nf_per_task[i]] for i in 1:model.nTask]
        return μf,σ²f
    end
end

function _predict_f(model::MCGP{T,<:Likelihood,<:GibbsSampling},X_test::AbstractMatrix{T};covf::Bool=true,fullcov::Bool=false) where {T}
    k_star = get_σ_k(model).*kernelmatrix.(get_kernel(model),[X_test],get_Z(model),obsdim=1)
    f = _sample_f(model,X_test,k_star)
    μf = Tuple(vec(mean(f[k],dims=2)) for k in 1:model.nLatent)
    if !covf
        return (μf,)
    end
    if fullcov
        k_starstar = get_σ_k(model).*kernelmatrix.(get_kernel(model),[X_test],obsdim=1)
        Σf = Symmetric.(k_starstar .- invquad.(get_K(model),k_star) .+  cov.(f))
        return μf, Σf
    else
        global k_starstar = get_σ_k(model).*kerneldiagmatrix.(get_kernel(model),[X_test],obsdim=1) .+ [T(jitter)*ones(T,size(X_test,1))]
        σ²f = k_starstar .- opt_diag.(k_star,get_K(model).\k_star) .+  diag.(cov.(f,dims=2))
        return μf,σ²f
    end
end

function _sample_f(model::MCGP{T,<:Likelihood,<:GibbsSampling},X_test::AbstractMatrix{T},k_star=get_σ_k(model).*kernelmatrix.(get_kernel(model),[X_test],get_Z(model),obsdim=1)) where {T}
    return f = [k_star[k]*(model.f[k].K\model.inference.sample_store[:,:,k]') for k in 1:model.nLatent]
end

predict_f(model::AbstractGP{T},X_test::AbstractVector{T};covf::Bool=false,fullcov::Bool=false) where {T} = predict_f(model,reshape(X_test,length(X_test),1),covf=covf,fullcov=fullcov)

predict_f(model::AbstractGP{T,L,I,TGP,1},X_test::AbstractMatrix{T};covf::Bool=false,fullcov::Bool=false) where {T,L,I,TGP} = first.(_predict_f(model,X_test;covf=covf,fullcov=fullcov))

predict_f(model::AbstractGP,X_test::AbstractMatrix{T};covf::Bool=false,fullcov::Bool=false) where {T,L,I,TGP} = _predict_f(model,X_test;covf=covf,fullcov=fullcov)

## Wrapper to predict vectors ##
function predict_y(model::AbstractGP{T},X_test::AbstractVector{T}) where {T}
    return predict_y(model,reshape(X_test,:,1))
end

"""
`predict_y(model::AbstractGP,X_test::AbstractMatrix)`

Return
    - the predictive mean of `X_test` for regression
    - the sign of `X_test` for classification
    - the most likely class for multi-class classification
    - the expected number of events for an event likelihood
"""
function predict_y(model::AbstractGP{T},X_test::AbstractMatrix{T}) where {T}
    return predict_y(model.likelihood,_predict_f(model,X_test,covf=false)[1])
end

predict_y(model::MOSVGP,X_test::AbstractMatrix) = predict_y.(model.likelihood,_predict_f(model,X_test,covf=false))

predict_y(l::RegressionLikelihood,μ::AbstractVector{<:Real}) = μ
predict_y(l::RegressionLikelihood,μ::AbstractVector{<:AbstractVector}) = first(μ)
predict_y(l::ClassificationLikelihood,μ::AbstractVector{<:Real}) = sign.(μ)
predict_y(l::ClassificationLikelihood,μ::AbstractVector{<:AbstractVector}) = sign.(first(μ))
predict_y(l::MultiClassLikelihood,μs::AbstractVector{<:AbstractVector{<:Real}}) = [l.class_mapping[argmax([μ[i] for μ in μs])] for i in 1:length(μs[1])]
predict_y(l::MultiClassLikelihood,μs::AbstractVector{<:AbstractVector{<:AbstractVector{<:Real}}}) = predict_y(l,first(μs))
predict_y(l::EventLikelihood,μ::AbstractVector{<:Real}) = expec_count(l,μ)
predict_y(l::EventLikelihood,μ::AbstractVector{<:AbstractVector}) = expec_count(l,first(μ))

## Wrapper to return proba on vectors ##
proba_y(model::AbstractGP{T},X_test::AbstractVector{T}) where {T} = proba_y(model,reshape(X_test,:,1))

"""
`proba_y(model::AbstractGP,X_test::AbstractMatrix)`

Return the probability distribution p(y_test|model,X_test) :

    - Tuple of vectors of mean and variance for regression
    - Vector of probabilities of y_test = 1 for binary classification
    - Dataframe with columns and probability per class for multi-class classification
"""
function proba_y(model::AbstractGP,X_test::AbstractMatrix)
    μ_f,Σ_f = _predict_f(model,X_test,covf=true)
    pred = compute_proba(model.likelihood,μ_f,Σ_f)
end

function proba_y(model::MOSVGP,X_test::AbstractMatrix)
    μ_f,Σ_f = _predict_f(model,X_test,covf=true)
    preds = compute_proba.(model.likelihood,μ_f,Σ_f)
end

function proba_y(model::AbstractGP{T,<:MultiClassLikelihood},X_test::AbstractMatrix) where {T}
    μ_f,Σ_f = _predict_f(model,X_test,covf=true)
    μ_p = compute_proba(model.likelihood,μ_f,Σ_f)
end

compute_proba(l::Likelihood,μ::AbstractVector{<:AbstractVector},σ²::AbstractVector{<:AbstractVector}) = compute_proba(l,first(μ),first(σ²))

function proba_y(model::MCGP{T,<:Union{<:RegressionLikelihood{T},<:ClassificationLikelihood{T}},<:GibbsSampling},X_test::AbstractMatrix{T};nSamples::Int=200) where {T<:Real}
    N_test = size(X_test,1)
    f = _sample_f(model,X_test)
    k_starstar = kerneldiagmatrix.([X_test],model.kernel)
    K̃ = k_starstar .- opt_diag.(k_star.*model.invKnn,k_star) .+ [zeros(size(X_test,1)) for i in 1:model.nLatent]
    nf = length(model.inference.sample_store[1])
    proba = [zeros(size(X_test,1)) for i in 1:model.nLatent]
    sig_proba = [zeros(size(X_test,1)) for i in 1:model.nLatent]
    for i in 1:nf
        for k in 1:model.nLatent
            proba[k], sig_proba[k] = (proba[k],sig_proba[k]) .+ compute_proba(model.likelihood, getindex.(f,[i])[k],K̃[k])
        end
    end
    if model.nLatent == 1
        return (proba[1]/nf, sig_proba[1]/nf)
    else
        return (proba./nf, sig_proba./nf)
    end
end

function proba_y(model::VGP{T,<:MultiClassLikelihood{T},<:GibbsSampling{T}},X_test::AbstractMatrix{T};nSamples::Int=200) where {T}
    k_star = kernelmatrix.([X_test],[model.inference.x],model.kernel)
    f = [[k_star[min(k,model.nPrior)]*model.invKnn[min(k,model.nPrior)]].*model.inference.sample_store[k] for k in 1:model.nLatent]
    k_starstar = kerneldiagmatrix.([X_test],model.kernel)
    K̃ = k_starstar .- opt_diag.(k_star.*model.invKnn,k_star) .+ [zeros(size(X_test,1)) for i in 1:model.nLatent]
    nf = length(model.inference.sample_store[1])
    proba = zeros(size(X_test,1),model.nLatent)
    labels = Array{Symbol}(undef,model.nLatent)
    for i in 1:nf
        res = compute_proba(model.likelihood,getindex.(f,[i]),K̃,nSamples)
        if i ==  1
            labels = names(res)
        end
        proba .+= Matrix(res)
    end
    return DataFrame(proba/nf,labels)
end
#
function compute_proba(l::Likelihood{T},μ::AbstractVector{T},σ²::AbstractVector{}) where {T<:Real}
    @error "Non implemented for the likelihood $l"
end


text:
module TestHelpers

using Base.Test
using JuMP
using MathProgBase

using MIPVerify
using MIPVerify: find_adversarial_example
using MIPVerify: NeuralNet
using MIPVerify: PerturbationFamily

const TEST_DEFAULT_TIGHTENING_ALGORITHM = lp

if Pkg.installed("Gurobi") == nothing
    using Cbc
    main_solver = CbcSolver(logLevel=0)
    tightening_solver = CbcSolver(logLevel=0, seconds=20)
else
    using Gurobi
    main_solver = GurobiSolver(Gurobi.Env())
    tightening_solver = GurobiSolver(Gurobi.Env(), OutputFlag=0, TimeLimit=20)
end

function get_main_solver()::MathProgBase.SolverInterface.AbstractMathProgSolver
    main_solver
end

function get_tightening_solver()::MathProgBase.SolverInterface.AbstractMathProgSolver
    tightening_solver
end

function get_new_model()::Model
    solver = get_main_solver()
    MathProgBase.setparameters!(solver, Silent = true)
    return Model(solver=solver)
end

"""
Tests the `find_adversarial_example` function.
  - If `x0` is already classified in the target label, `expected_objective_value`
    should be 0.
  - If there is no adversarial example for the specified parameters,
    `expected_objective_value` should be NaN.
  - If there is an adversarial example, checks that the objective value is as expected,
    and that the perturbed output for the target label exceeds the perturbed
    output for all other labels by `tolerance`.
"""
function test_find_adversarial_example(
    nn::NeuralNet,
    x0::Array{<:Real, N},
    target_selection::Union{Integer, Array{<:Integer, 1}},
    pp::PerturbationFamily,
    norm_order::Real,
    tolerance::Real,
    expected_objective_value::Real,
    ) where {N}
    d = find_adversarial_example(
        nn, x0, target_selection, get_main_solver(),
        pp = pp, norm_order = norm_order, tolerance = tolerance, rebuild=false,
        tightening_solver=get_tightening_solver(), tightening_algorithm=TEST_DEFAULT_TIGHTENING_ALGORITHM)
    println(d[:SolveStatus])
    if d[:SolveStatus] == :Infeasible || d[:SolveStatus] == :InfeasibleOrUnbounded
        @test isnan(expected_objective_value)
    else
        actual_objective_value = getobjectivevalue(d[:Model])
        if expected_objective_value == 0
            @test isapprox(actual_objective_value, expected_objective_value; atol=1e-4)
        else
            @test isapprox(actual_objective_value, expected_objective_value; rtol=5e-5)

            perturbed_output = getvalue(d[:PerturbedInput]) |> nn
            perturbed_target_output = maximum(perturbed_output[Bool[i∈d[:TargetIndexes] for i = 1:length(d[:Output])]])
            maximum_perturbed_other_output = maximum(perturbed_output[Bool[i∉d[:TargetIndexes] for i = 1:length(d[:Output])]])
            @test perturbed_target_output/(maximum_perturbed_other_output+tolerance)≈1 atol=5e-5
        end
    end
end

"""
Runs tests on the neural net described by `nn` for input `x0` and the objective values
indicated in `expected objective values`.

# Arguments
- `expected_objective_values::Dict`:
   d[(target_selection, perturbation_parameter, norm_order, tolerance)] = expected_objective_value
   `expected_objective_value` is `NaN` if there is no perturbation that brings the image
   into the target category.
"""
function batch_test_adversarial_example(
    nn::NeuralNet,
    x0::Array{<:Real, N},
    expected_objective_values::Dict
) where {N}
    for (test_params, expected_objective_value) in expected_objective_values
        (target_selection, pp, norm_order, tolerance) = test_params
        @testset "target label = $target_selection, $(string(pp)) perturbation, norm order = $norm_order, tolerance = $tolerance" begin
            test_find_adversarial_example(
                nn, x0,
                target_selection, pp, norm_order, tolerance, expected_objective_value)
            end
        end
    end
end


text:
# Some fft routine like fftfreq
include("fft_util.jl")

"""
    function get_points(N_pix::Int, sky::Array{Float64, 2}; direction::Bool)
# Brief Description
This function return the pointing in accord to the sky scanning-strategy.
Now return the pointing directions as pixels. You can choose between two
different scanning strategy mode: if `direction = true` the function return a
left-right scanning. If `direction = false` the return is a up-down scanning.

Future improvements: I ought code a new `get_points` function in order to
get true sky-direction in Right-ascension and Declination coordinates

# Development
The function as is, has not presented weird behaviours.

# Args description
N_pix:      Number of pixel, in linear coordinate
sky  :      The sky that you want observe, in Array{Float64, 2} format
direction:  If true = Left-Right scanning else Up-Down.
"""
function get_points(N_pix::Int; direction::Bool)

    points = Array{Float64, 1}(undef, 0)

    if direction == true
        x_pix = Array{Int64, 1}(undef, 0)
        r1 = Int.(round.(range(1, N_pix, length=N_pix)))
        r2 = Int.(round.(range(N_pix, 1, length=N_pix)))

        append!(x_pix, r1)
        append!(x_pix, r2)

        x_pix = repeat(x_pix, Int.(round.(N_pix/2)))

        y_pix = Int.(round.(reshape(reshape(repeat(range(
                                   1, N_pix, length = N_pix),
                             N_pix),
                    N_pix, N_pix)',
            1, N_pix*N_pix)))

        append!(points, x_pix[:,1])
        append!(points, y_pix[1,:])
    else
        x_pix = Int.(round.(reshape(reshape(repeat(range(
                                   1, N_pix, length = N_pix),
                             N_pix),
                    N_pix, N_pix)',
            1, N_pix*N_pix)))

        y_pix = Array{Int64, 1}(undef, 0)

        r1 = Int.(round.(range(1, N_pix, length=N_pix)))
        r2 = Int.(round.(range(N_pix, 1, length=N_pix)))
        append!(y_pix, r1)
        append!(y_pix, r2)

        y_pix = repeat(y_pix, Int.(round.(N_pix/2)))

        append!(points, x_pix[1,:])
        append!(points, y_pix[:,1])

    end

    points = reshape(points, Int(round(N_pix^2)), 2)

    return points

end

"""
    function get_tod(sky::Array{Float64, 2}, points::Array{Int64, 2})

# Brief description
Given a set of points, genereted either L-R or U-D scanning, return the
time order data (TOD) series detected by a sensors.

# Development
The function is very simple and has not shown weird behaviours.

# Args description
sky  :      The sky that you want observe, in Array{Float64, 2} format
points:     Time-stream of pointings of your sky scanning-strategy.

"""
function get_tod(sky::Array{Float64, 2}, points::Array{Int64, 2})

    Num_of_pixels = size(points)[1]
    tod = [sky[points[i,1], points[i,2]] for i = 1:Num_of_pixels]

    return tod

end

"""
    function observe_sky(N_pix::Int, sky::Array{Float64, 2};
                         noise::Bool, direction_l_r::Bool)
# Brief description
This function return a noisy tod. The noise is a superimposition of a
one-over-f like: include the instrumental fluctuation and atmospheric
noise, and by a second part of white noise

# Args description
N_pix:      Number of pixel, in linear coordinate
sky  :      The sky that you want observe, in Array{Float64, 2} format
direction:  If true = Left-Right scanning else Up-Down.
noise :     Is a useless args ... I ought clean this mess!

"""
function observe_sky(N_pix::Int, sky::Array{Float64, 2};
                     noise::Bool, direction_l_r::Bool)

    points  = Int.(round.(get_points(N_pix, direction=direction_l_r)))
    tod = get_tod(sky, points)

    if noise==true
        rng  = MersenneTwister()
        Frand = fft(randn(rng, length(tod)))
        #Frand2 = fft(randn(rng, length(tod)))
        noise_spec = tod_noise(length(tod), 0.0083, 0.1, 3.0, 20.76)
        #noise_spec2 = tod_noise(length(tod), 0.0083, 1.0, 2.0, 20.76)
        fnoise = Frand .* (noise_spec .^0.5)
        #fnoise2 = Frand2 .* (noise_spec2 .^0.5)
        tod += real.(ifft(fnoise))
        #tod += real.(ifft(fnoise2))
        dataset = TOD_fake_pointing(tod, points, noise_spec)
    else
        noise_spec = tod_noise(length(tod), 0.0083, 0.1, 3.0, 20.76)
        dataset = TOD_fake_pointing(tod, points, noise_spec)
    end

    return dataset

end

function tod_noise(Nsamp::Int, dt::Float64, fknee::Float64, alpha::Float64,
                   sigma::Float64)
    Freq = abs.(fftfreq(Nsamp, dt))
    return (1 .+ ( max.(Freq, Freq[2]) ./ fknee).^(-alpha)) .* (sigma^2)
end

"""
    function binned_map(N::Int, tod::Array{Float64, 1},
                        pointing::Array{Int64, 2})
# Brief description
The function returns the baselines sum, named Sky Binned Map. If P is the
pointing matrix, returns: PᵀP

# Args description
N_pix:      Number of pixel, in linear coordinate
tod  :      Time Order Data stream, returned by get_tod function
points:     Time-stream of pointings of your sky scanning-strategy.
"""
function binned_map(N_pix::Int, tod::Array{Float64, 1},
                    points::Array{Int64, 2})
    map = zeros(N_pix, N_pix)
    hit = zeros(N_pix, N_pix)

    for i = 1:length(tod)
        map[points[i, 1], points[i, 2]] += tod[i]
        hit[points[i, 1], points[i, 2]] += 1.0
    end

    return map, hit

end


text:
# trivial conversion from PostgreSQL->SQLite

using LibPQ
using SQLite
using PostgresCatalog

#conn = LibPQ.Connection("postgresql:///synpuf-5pct?host=/var/run/postgresql")
conn = LibPQ.Connection("postgresql:///synpuf-10p?host=/var/run/postgresql")
db = SQLite.DB("synpuf-10p.sq3")

cat = PostgresCatalog.introspect(conn)

function visit(tbl)
    columns = String[]
    for col in tbl.columns
       if col.type.name in ("int4", "varchar", "int8",
               "numeric", "text")
           push!(columns, col.name)
       elseif col.type.name in ("date",)
           push!(columns, "$(col.name)::TEXT")
       else
           throw("unknown type $(col.type.name)")
       end
    end
    columns = join(columns, ",")
    res = execute(conn, "SELECT $columns FROM $(tbl.name)")
    SQLite.load!(res, db, tbl.name; analyze=false)
    SQLite.execute(db, "ANALYZE $(tbl.name);")
end

for tbl in cat["public"]
    println("tbl ", tbl.name)
    visit(tbl)
end


text:
module Lib

using Libdl: Libdl

using CEnum

const PKGDIR = dirname(dirname(@__DIR__))
const USRDIR = joinpath(PKGDIR, "deps", "dnnl")
const LIBDIR = joinpath(USRDIR, "lib")

const _flags = Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL

const libdnnl = joinpath(LIBDIR, "libdnnl.so")
Libdl.dlopen(libdnnl, _flags)

#include("types.jl")
include("dnnl.jl")
#include("dnnl_threadpool.jl")

end # module


text:
"""
    mirror(N::T) where {T <: AbstractUnipartiteNetwork}

Returns a mirrored version of the adjacency matrix, required for spectral
radius.
"""
function mirror(N::T) where {T <: AbstractUnipartiteNetwork}
    M = copy(N)
    for i in N
        M[i.to, i.from] = M[i.from, i.to]
    end
    return M
end

ρ_phillips(N::T, v::Float64) where {T <: UnipartiteNetwork}  = v/((links(N)*(richness(N)-1))/richness(N))^0.5
ρ_ska(N::T, v::Float64) where {T <: UnipartiteNetwork} = v/sqrt(links(N))
ρ_raw(N::T, v::Float64) where {T <: UnipartiteNetwork} = v

"""
    ρ(N::T; range=EcologicalNetworks.ρ_ska) where {T <: UnipartiteNetwork}

Returns the spectral radius (the absolute value of the largest real part of the
eigenvalues of the adjacency matrix) of any unipartite network whose
interactions are positive or null. The spectral radius is measured on the
mirrored version of the matrix, so that an interaction from i to j also implies
an interaction from j to i.

Note that the spectral radius has been suggested as a measure of nestedness by
Staniczenko *et al.* (2013). Phillips (2011) uses it as a measure of the ability
of a system to dampen or absorb perturbations.

#### Maximal values

The spectral radius is sensitive to network size, and to a certain extent to the
number of links. The keyword argument `range` will divide return a ranged
version of the spectral radius, so that it is expressed relatively to its
maximal value. The `range` argument takes a *function*, which requires two
arguments: the network (which must be unipartite), and the value of the spectral
radius.

Options that come with `EcologicalNetworks.jl` (where L is the number of links
and S the number of nodes) are:

1. `EcologicalNetworks.ρ_phillips`: divides by the square root of (2L(S-1))/S,
   as in Phillips (20110)
1. `EcologicalNetworks.ρ_ska`: divides by the square root of L, as in
   Staniczenko *et al.* (2013) - this is the **default**
1. `EcologicalNetworks.ρ_raw`: returns the raw value

#### References

- Phillips, J.D., 2011. The structure of ecological state transitions:
  Amplification, synchronization, and constraints in responses to environmental
  change. Ecological Complexity, 8, 336–346.
  https://doi.org/10.1016/j.ecocom.2011.07.004

- Staniczenko, P.P.A., Kopp, J.C., Allesina, S., 2013. The ghost of nestedness
  in ecological networks. Nat Commun 4, 1391. https://doi.org/10.1038/ncomms2422
"""
function ρ(N::T; range=EcologicalNetworks.ρ_ska) where {T <: UnipartiteNetwork}
    links(N) == 0 && return NaN
    M = mirror(N)
    @assert minimum(M.edges) ≥ zero(_interaction_type(M))
    absolute_real_part = abs.(real.(eigvals(Array(M.edges))))
    return range(M, maximum(absolute_real_part))
end

"""
    ρ(N::T; varargs...) where {T <: BipartiteNetwork}

Bipartite version of the spectral radius. In practice, this casts the network
into its unipartite representation, since the spectral radius only makes sense
for square matrices.
"""
function ρ(N::T; varargs...) where {T <: BipartiteNetwork}
    return ρ(convert(UnipartiteNetwork, N); varargs...)
end


text:
# This file was generated, do not modify it. # hide
tree = machine(tree_model, X, y)


text:
import Pkg
Pkg.activate(@__DIR__)

using TestEnv

Pkg.activate(joinpath(@__DIR__, ".."))
TestEnv.activate()

using FrankWolfe
using ProgressMeter
using Arpack
using Plots
using DoubleFloats
using ReverseDiff


text:
include("mpl_colormaps.jl")

function mpl_conv(cm_data)
  N = size(cm_data,1)
  colormap = Array(RGB{Float64},N)
  for i in 1:N
    colormap[i] = RGB{Float64}(cm_data[i,:]...)
  end
  colormap
end

# macro might help here
mpl_inferno() = mpl_conv(inferno_data)
access_dict["inferno"] = mpl_inferno

mpl_viridis() = mpl_conv(viridis_data)
access_dict["viridis"] = mpl_viridis

mpl_magma() = mpl_conv(magma_data)
access_dict["magma"] = mpl_magma

mpl_plasma() = mpl_conv(plasma_data)
access_dict["plasma"] = mpl_plasma

max_stars_repo_name: lobingera/ColorMaps


text:
export kaczmarzUpdated

mutable struct KaczmarzUpdated{matT,T,U,Tsparse} <: AbstractLinearSolver
  S::matT
  u::Vector{T}
  reg::Vector{Regularization}
  denom::Vector{U}
  rowindex::Vector{Int64}
  rowIndexCycle::Vector{Int64}
  cl::Vector{T}
  vl::Vector{T}
  εw::Vector{T}
  τl::T
  αl::T
  weights::Vector{U}
  enforceReal::Bool
  enforcePositive::Bool
  Randomized::Bool
  SubMatrixSize::Int64
  Probabilities::Weights
  shuffleRows::Bool
  seed::Int64
  sparseTrafo::Tsparse
  iterations::Int64
  constraintMask::Union{Nothing,Vector{Bool}}
end

"""
KaczmarzUpdated(S, b=nothing; λ::Real=0.0, reg = Regularization("L2", λ)
            , weights::Vector{R}=ones(Float64,size(S,1))
            , sparseTrafo=nothing
            , enforceReal::Bool=false
            , enforcePositive::Bool=false
            , Randomized::Bool=false
            , SubMatrixSize::Int64=150
            , shuffleRows::Bool=false
            , seed::Int=1234
            , iterations::Int64=10
            , kargs...) where R <: Real

creates a KaczmarzUpdated object

# Arguments
* `S`                                             - system matrix
* `b=nothing`                                     - measurement
* (`λ::Real=0.0`)                                 - regularization parameter
* (`reg=Regularization("L2", λ)`)                 - Regularization object
* (`weights::Vector{R}=ones(Float64,size(S,1))`) - weights for the data term
* (`sparseTrafo=nothing`)                         - sparsifying transformation
* (`enforceReal::Bool=false`)                     - constrain the solution to be real
* (`enforcePositive::Bool=false`)                 - constrain the solution to have positive real part
* (`Randomized::Bool=false`)                      - randomize Kacmarz algorithm
* (`SubMatrixSize::Int64=150`)                    - number of rows used in Randomized Kaczmarz algorithm
* (`shuffleRows::Bool=false`)                     - randomize Kacmarz algorithm
* (`seed::Int=1234`)                              - seed for randomized algorithm
* (iterations::Int64=10)                          - number of iterations
"""

function KaczmarzUpdated(S; b=nothing, λ=[0.0], reg = nothing
    , regName = ["L2"]
    , weights=nothing
    , sparseTrafo=nothing
    , enforceReal::Bool=false
    , enforcePositive::Bool=false
    , Randomized::Bool=false
    , SubMatrixSize::Int64=150
    , shuffleRows::Bool=false
    , seed::Int=1234
    , iterations::Int64=10
    , constraintMask=nothing
    , kargs...)

if typeof(λ) <: Number
λ = [λ]
end

if reg == nothing
reg = Regularization(regName, λ; kargs...)
end

if regName[1] != "L2"
error("Kaczmarz only supports L2 regularizer")
end

T = real(eltype(S))

# make sure weights are not empty
w = (weights!=nothing ? weights : ones(T,size(S,1)))

# setup denom and rowindex
denom, rowindex = initkaczmarz(S, λ[1], w)
rowIndexCycle=collect(1:length(rowindex))

Probabilities = Find_Rows_Probaboloties(S)

M,N = size(S)
if b != nothing
u = b
else
u = zeros(eltype(S),M)
end
cl = zeros(eltype(S),N)
vl = zeros(eltype(S),M)
εw = zeros(eltype(S),length(rowindex))
τl = zero(eltype(S))
αl = zero(eltype(S))

return KaczmarzUpdated(S,u,reg,denom,rowindex,rowIndexCycle,cl,vl,εw,τl,αl
        ,T.(w),enforceReal,enforcePositive,Randomized,SubMatrixSize,Probabilities
        ,shuffleRows,Int64(seed),sparseTrafo,iterations, constraintMask)
end

"""
  init!(solver::KaczmarzUpdated
              ; S::matT=solver.S
              , u::Vector{T}=T[]
              , cl::Vector{T}=T[]
              , shuffleRows=solver.shuffleRows) where {T,matT}

(re-) initializes the CGNR iterator
"""
function init!(solver::KaczmarzUpdated
              ; S::matT=solver.S
              , u::Vector{T}=T[]
              , cl::Vector{T}=T[]
              , weights::Vector{R}=solver.weights
              , shuffleRows=solver.shuffleRows) where {T,matT,R}

  if S != solver.S
    solver.denom, solver.rowindex = initkaczmarz(S,solver.reg.λ,weights)
    solver.rowIndexCycle = collect(1:length(solver.rowindex))
  end

  if shuffleRows
    Random.seed!(solver.seed)
    shuffle!(solver.rowIndexCycle)
  end
  solver.u[:] .= u
  solver.weights=weights

  # start vector
  if isempty(cl)
    solver.cl[:] .= zero(T)
  else
    solver.cl[:] .= cl
  end
  solver.vl[:] .= zero(T)

  for i=1:length(solver.rowindex)
    j = solver.rowindex[i]
    solver.ɛw[i] = sqrt(solver.reg[1].λ) / weights[j]
  end
end

"""
  solve(solver::KaczmarzUpdated, u::Vector{T};
                S::matT=solver.S, startVector::Vector{T}=eltype(S)[]
                , weights::Vector=solver.weights, shuffleRows::Bool=false
                , solverInfo=nothing, kargs...) where {T,matT}

solves Thikhonov-regularized inverse problem using Kaczmarz algorithm.

# Arguments
* `solver::KaczmarzUpdated  - the solver containing both system matrix and regularizer
* `u::Vector`        - data vector
* (`S::matT=solver.S`)                  - operator for the data-term of the problem
* (`startVector::Vector{T}=T[]`)        - initial guess for the solution
* (`weights::Vector{T}=solver.weights`) - weights for the data term
* (`shuffleRows::Bool=false`)           - randomize Kacmarz algorithm
* (`solverInfo=nothing`)                - solverInfo for logging

when a `SolverInfo` objects is passed, the residuals are stored in `solverInfo.convMeas`.
"""
function solve(solver::KaczmarzUpdated, u::Vector{T};
                S::matT=solver.S, startVector::Vector{T}=eltype(S)[]
                , weights::Vector=solver.weights, shuffleRows::Bool=false
                , solverInfo=nothing, kargs...) where {T,matT}

  # initialize solver parameters
  init!(solver; S=S, u=u, cl=startVector, weights=weights, shuffleRows=shuffleRows)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.cl,solver.iterations)

  # perform Kaczmarz iterations
  for item in solver
    solverInfo != nothing && storeInfo(solverInfo,solver.cl,solver.iterations)
  end

  return solver.cl
end

function iterate(solver::KaczmarzUpdated, iteration::Int=0)
    if done(solver,iteration) return nothing end

    if solver.Randomized
        UsedIndex=Int.(StatsBase.sample!(Random.GLOBAL_RNG, solver.rowIndexCycle, solver.Probabilities, zeros(solver.SubMatrixSize), replace=false))
    else
        UsedIndex=solver.rowIndexCycle
    end

    for i in UsedIndex
      j = i #solver.rowindex[i]
      solver.τl = dot_with_matrix_row(solver.S,solver.cl,j)
      solver.αl = solver.denom[i]*(solver.u[j]-solver.τl-solver.ɛw[i]*solver.vl[j])
      kaczmarz_update!(solver.S,solver.cl,j,solver.αl)
      solver.vl[j] += solver.αl*solver.ɛw[i]
    end

    # invoke constraints
    applyConstraints(solver.cl, solver.sparseTrafo,
                                solver.enforceReal,
                                solver.enforcePositive,
                                solver.constraintMask)

    if length(solver.reg) > 1
      # We skip the L2 regularizer, since it has already been applied
      prox!(solver.cl, solver.reg[2:end])
    end

    return solver.vl, iteration+1

  end

  @inline done(solver::KaczmarzUpdated,iteration::Int) = iteration>=solver.iterations

  """
  This function finds the Probabilities of the rows of the system matrix
  """

  function Find_Rows_Probaboloties(S::AbstractMatrix)
    M,N = size(S)
    Norm = norm(S)
    P = Vector{Float64}()
    for i=1:M
      prob = (norm(S[i,:]))^2/(Norm)^2
      push!(P,prob)
    end
    return weights(P)
  end

  ### initKaczmar ###

  """
    initkaczmarz(S::AbstractMatrix,λ,weights::Vector)

This funtion saves the denominators to compute αl in denom and the rowindices,
which lead to an update of cl in rowindex.
"""
function initkaczmarz(S::AbstractMatrix,λ,weights::Vector)
  T = typeof(real(S[1]))
  denom = T[]
  rowindex = Int64[]

  for i=1:size(S,1)
    s² = rownorm²(S,i)*weights[i]^2
    if s²>0
      push!(denom,weights[i]^2/(s²+λ))
      push!(rowindex,i)
    end
  end
  denom, rowindex
end

### kaczmarz_update! ###

"""
    kaczmarz_update!(A::DenseMatrix{T}, x::Vector, k::Integer, beta) where T

This funtion updates x during the kaczmarz algorithm for dense matrices.
"""
function kaczmarz_update!(A::DenseMatrix{T}, x::Vector, k::Integer, beta) where T
  @simd for n=1:size(A,2)
    @inbounds x[n] += beta*conj(A[k,n])
  end
end

"""
    kaczmarz_update!(B::Transpose{T,S}, x::Vector,
                     k::Integer, beta) where {T,S<:DenseMatrix}

This funtion updates x during the kaczmarz algorithm for dense matrices.
"""
function kaczmarz_update!(B::Transpose{T,S}, x::Vector,
			  k::Integer, beta) where {T,S<:DenseMatrix}
  A = B.parent
  @inbounds @simd for n=1:size(A,1)
      x[n] += beta*conj(A[n,k])
  end
end

# kaczmarz_update! with manual simd optimization
for (T,W, WS,shufflevectorMask,vσ) in [(Float32,:WF32,:WF32S,:shufflevectorMaskF32,:vσF32),(Float64,:WF64,:WF64S,:shufflevectorMaskF64,:vσF64)]
    eval(quote
        const $WS = VectorizationBase.pick_vector_width($T)
        const $W = Int(VectorizationBase.pick_vector_width($T))
        const $shufflevectorMask = Val(ntuple(k -> iseven(k-1) ? k : k-2, $W))
        const $vσ = Vec(ntuple(k -> (-1f0)^(k+1),$W)...)
        function kaczmarz_update!(A::Transpose{Complex{$T},S}, b::Vector{Complex{$T}}, k::Integer, beta::Complex{$T}) where {S<:DenseMatrix}
            b = reinterpret($T,b)
            A = reinterpret($T,A.parent)

            N = length(b)
            Nrep, Nrem = divrem(N,4*$W) # main loop
            Mrep, Mrem = divrem(Nrem,$W) # last iterations
            idx = MM{$W}(1)
            iOffset = 4*$W

            vβr = vbroadcast($WS, beta.re) * $vσ # vector containing (βᵣ,-βᵣ,βᵣ,-βᵣ,...)
            vβi = vbroadcast($WS, beta.im) # vector containing (βᵢ,βᵢ,βᵢ,βᵢ,...)

            GC.@preserve b A begin # protect A and y from GC
                vptrA = stridedpointer(A)
                vptrb = stridedpointer(b)
                for _ = 1:Nrep
                    Base.Cartesian.@nexprs 4 i -> vb_i = vload(vptrb, ($W*(i-1) + idx,))
                    Base.Cartesian.@nexprs 4 i -> va_i = vload(vptrA, ($W*(i-1) + idx,k))
                    Base.Cartesian.@nexprs 4 i -> begin
                        vb_i = muladd(va_i, vβr, vb_i)
                        va_i = shufflevector(va_i, $shufflevectorMask)
                        vb_i = muladd(va_i, vβi, vb_i)
                    	vstore!(vptrb, vb_i, ($W*(i-1) + idx,))
                    end
                    idx += iOffset
                end

                for _ = 1:Mrep
	            vb = vload(vptrb, (idx,))
	            va = vload(vptrA, (idx,k))
                    vb = muladd(va, vβr, vb)
                    va = shufflevector(va, $shufflevectorMask)
                    vb = muladd(va, vβi, vb)
		            vstore!(vptrb, vb, (idx,))
                    idx += $W
                end

                if Mrem!=0
                    vloadMask = VectorizationBase.mask($T, Mrem)
                    vb = vload(vptrb, (idx,), vloadMask)
                    va = vload(vptrA, (idx,k), vloadMask)
                    vb = muladd(va, vβr, vb)
                    va = shufflevector(va, $shufflevectorMask)
                    vb = muladd(va, vβi, vb)
                    vstore!(vptrb, vb, (idx,), vloadMask)
                end
            end # GC.@preserve
        end
    end)
end

#=
@doc "This funtion updates x during the kaczmarz algorithm for dense matrices." ->
function kaczmarz_update!{T}(A::Matrix{T}, x::Vector{T}, k::Integer, beta::T)
  BLAS.axpy!(length(x), beta, pointer(A,sub2ind(size(A),1,k)), 1, pointer(x), 1)
end
=#

"""
    kaczmarz_update!(B::Transpose{T,S}, x::Vector,
                          k::Integer, beta) where {T,S<:SparseMatrixCSC}

This funtion updates x during the kaczmarz algorithm for sparse matrices.
"""
function kaczmarz_update!(B::Transpose{T,S}, x::Vector,
                          k::Integer, beta) where {T,S<:SparseMatrixCSC}
  A = B.parent
  N = A.colptr[k+1]-A.colptr[k]
  for n=A.colptr[k]:N-1+A.colptr[k]
    @inbounds x[A.rowval[n]] += beta*conj(A.nzval[n])
  end
end

### kaczmarz ###

"""
    kaczmarz(S, u::Vector; kargs...)

This funtion implements the kaczmarz algorithm.

# Keyword/Optional Arguments
* `λ::Float64`: The regularization parameter, relative to the matrix trace
* `iterations::Int`: Number of iterations of the iterative solver
* `solver::AbstractString`: Algorithm used to solve the imaging equation (currently "kaczmarz" or "cgnr")
* `normWeights::Bool`: Enable row normalization (true/false)
* `sparseTrafo::AbstractString`: Enable sparseTrafo if set to "DCT-IV" or "FFT".
* `shuff::Bool`: Enable shuffeling of rows during iterations in the kaczmarz algorithm.
* `enforceReal::Bool`: Enable projection of solution on real plane during iteration.
* `enforcePositive::Bool`: Enable projection of solution onto positive halfplane during iteration.
"""


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AuthenticationV1beta1Api <: SwaggerApi
    client::Swagger.Client
end

"""

create a TokenReview
Param: body::IoK8sApiAuthenticationV1beta1TokenReview (required)
Param: dryRun::String
Param: fieldManager::String
Param: pretty::String
Return: IoK8sApiAuthenticationV1beta1TokenReview
"""
function createAuthenticationV1beta1TokenReview(_api::AuthenticationV1beta1Api, body; dryRun=nothing, fieldManager=nothing, pretty=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", IoK8sApiAuthenticationV1beta1TokenReview, "/apis/authentication.k8s.io/v1beta1/tokenreviews", ["BearerToken"], body)
    Swagger.set_param(_ctx.query, "dryRun", dryRun)  # type String
    Swagger.set_param(_ctx.query, "fieldManager", fieldManager)  # type String
    Swagger.set_param(_ctx.query, "pretty", pretty)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

get available resources
Return: IoK8sApimachineryPkgApisMetaV1APIResourceList
"""
function getAuthenticationV1beta1APIResources(_api::AuthenticationV1beta1Api; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApimachineryPkgApisMetaV1APIResourceList, "/apis/authentication.k8s.io/v1beta1/", ["BearerToken"])
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"] : [_mediaType])
    Swagger.exec(_ctx)
end

export createAuthenticationV1beta1TokenReview, getAuthenticationV1beta1APIResources


text:
"""
A store that wraps any other AbstractStore but has access to the consolidated metadata
stored in the .zmetadata key. Whenever data attributes or metadata are accessed, the
data will be read from the dictionary instead.
"""
struct ConsolidatedStore{P} <: AbstractStore
  parent::P
  cons::Dict{String}
end
function ConsolidatedStore(s::AbstractStore)
  d = s[".zmetadata"]
  if d === nothing
    throw(ArgumentError("Could not find consolidated metadata for store $s"))
  end
  ConsolidatedStore(s,JSON.parse(String(copy(d)))["metadata"])
end

function Base.show(io::IO,d::ConsolidatedStore)
    b = IOBuffer()
    show(b,d.parent)
    print(io, "Consolidated ", String(take!(b)))
end

storagesize(d::ConsolidatedStore) = storagesize(d.parent)
zname(s::ConsolidatedStore) = zname(s.parent)
Base.getindex(d::ConsolidatedStore,i::String) = d.parent[i]
getmetadata(d::ConsolidatedStore) = Metadata(d.cons[".zarray"])
getattrs(d::ConsolidatedStore) = get(d.cons,".zattrs", Dict{String,Any}())
is_zarray(d::ConsolidatedStore) = haskey(d.cons,".zarray")
is_zgroup(d::ConsolidatedStore) = haskey(d.cons,".zgroup")
check_consolidated_write(i::String) = split(i,'/')[end] in (".zattrs",".zarray",".zgroup") &&
    throw(ArgumentError("Can not modify consolidated metadata, please re-open the dataset with `consolidated=false`"))
function Base.setindex!(d::ConsolidatedStore,v,i::String)
  #Here we check that we don't overwrite consolidated information
  check_consolidated_write(i)
  d.parent.a[i] = v
end
function Base.delete!(d::ConsolidatedStore,i::String)
  check_consolidated_write(i)
  delete!(d.a,i)
end
function subdirs(d::ConsolidatedStore)
  o = Set{String}()
  foreach(collect(keys(d.cons))) do k
    tr = split(k,'/')
    length(tr) > 1 && push!(o,first(tr))
  end
  collect(o)
end
Base.keys(d::ConsolidatedStore) = keys(d.parent)
newsub(d::ConsolidatedStore, n) = newsub(d.parent, n)
function getsub(d::ConsolidatedStore, n)
    newd = Dict{String,Any}()
    for (k,v) in d.cons
        s = split(k,'/')
        if s[1] == n
            newd[join(s[2:end],'/')] = v
        end
    end
    ConsolidatedStore(getsub(d.parent, n), newd)
end
path(d::ConsolidatedStore) = path(d.parent)

function consolidate_metadata(s::AbstractStore,d,prefix)
  for k in (".zattrs",".zarray",".zgroup")
    v = s[k]
    if v !== nothing
      d[string(prefix,k)] = JSON.parse(replace(String(Zarr.maybecopy(v)),": NaN,"=>": \"NaN\","))
    end
  end
  foreach(subdirs(s)) do subname
    ssub = getsub(s,subname)
    consolidate_metadata(ssub,d,string(prefix,subname,"/"))
  end
  d
end
function consolidate_metadata(s::AbstractStore)
  d = consolidate_metadata(s,Dict{String,Any}(),"")
  buf = IOBuffer()
  JSON.print(buf,Dict("metadata"=>d),4)
  s[".zmetadata"] = take!(buf)
  ConsolidatedStore(s)
end

consolidate_metadata(s) = consolidate_metadata(zopen(s))


text:
using HDF5

"""
    h5create(fname)

Create and return an hdf5 file at path `fname`. Files created by this function
will be compatible with `HDF5.start_swmr_write` and other Pope requirements."
"""
h5create(fname) = h5open(fname,"w", "libver_bounds", (HDF5.H5F_LIBVER_LATEST, HDF5.H5F_LIBVER_LATEST))
"""
    d_extend(d::HDF5Dataset, value::Vector, range::UnitRange)

Equivalent to `d[range]=value` on an extendible on dimensional HDF5Dataset `d` except
that the length of `d` is set to `maximum(range)` before writing.
"""
function d_extend(d::HDF5Dataset, value::Vector, range::UnitRange)
  set_dims!(d, (maximum(range),))
	d[range] = value
	d
end
"""
    g_require(parent::Union{HDF5File,HDF5Group}, name)

Return an HDF5Group with name `name` in `parent.` If the group does not exist, create it.
"""
function g_require(parent::Union{HDF5File,HDF5Group}, name)
	exists(parent,name) ? parent[name] : g_create(parent,name)
end

"""
    BufferedHDF5Dataset(ds::HDF5Dataset, v::Vector, lasti)

HDF5 appears to be inefficent for small writes, so this a simple buffer that
allows me to write to HDF5 only once per unit time (typically one second) to
limit the number of small writes. `MassCompatibleBufferedWriters` is based upon
`BufferedHDF5Dataset` and you can call `schedule` on an instand of `MassCompatibleBufferedWriters`
to set the write frequency.
"""
mutable struct BufferedHDF5Dataset{T}
  ds::HDF5Dataset
  v::Vector{T}
  lasti::Int64 # last index in hdf5 dataset
end
function BufferedHDF5Dataset{T}(g::Union{HDF5File,HDF5Group}, name, chunksize) where {T}
  ds = d_create(g, name, T, ((1,),(-1,)), "chunk", chunksize)
  BufferedHDF5Dataset{T}(ds, Vector{T}(), 0)
end
function write_to_hdf5(b::BufferedHDF5Dataset)
  r = b.lasti .+ (1:length(b.v))
  if length(r)>0
    d_extend(b.ds, b.v, r)
    empty!(b.v)
    b.lasti=last(r)
  end
  return
end
Base.write(b::BufferedHDF5Dataset{T},x::T) where {T} = push!(b.v,x)
Base.write(b::BufferedHDF5Dataset{T},x::Vector{T}) where {T} = append!(b.v,x)

abstract type BufferedWriter <: DataSink end
"""
    MassCompatibleBufferedWriters <: DataSink

Construct using `make_buffered_hdf5_writer(h5, channel_number, chunksize=1000, timeout_s=1.0)`.
Contains many `BufferedHDF5Dataset` and organizes writing to them. Supports
`schedule`, `stop`, `wait`, `close`, `flush`, `write`, `write_header`, `hdf5_file`,
`write_header_end`, `write_header_allchannel`. If you have an instance `d`,
then `d(x)` is equivalent to `write(d,x)`.
"""
mutable struct MassCompatibleBufferedWriters <: BufferedWriter
  endchannel        ::Channel{Bool}
  timeout_s         ::Float64
  task              ::Task
  filt_value        ::BufferedHDF5Dataset{Float32}
  filt_phase        ::BufferedHDF5Dataset{Float32}
  timestamp         ::BufferedHDF5Dataset{Float64}
  rowcount          ::BufferedHDF5Dataset{Int64}
  pretrig_mean      ::BufferedHDF5Dataset{Float32}
  pretrig_rms       ::BufferedHDF5Dataset{Float32}
  pulse_average     ::BufferedHDF5Dataset{Float32}
  pulse_rms         ::BufferedHDF5Dataset{Float32}
  rise_time         ::BufferedHDF5Dataset{Float32}
  postpeak_deriv    ::BufferedHDF5Dataset{Float32}
  peak_index        ::BufferedHDF5Dataset{UInt16}
  peak_value        ::BufferedHDF5Dataset{UInt16}
  min_value         ::BufferedHDF5Dataset{UInt16}
end
const mass_fieldnames = collect(fieldnames(MassCompatibleBufferedWriters)[4:end])
"    hdf5file(b::MassCompatibleBufferedWriters)
Return the filename of the hdf5 file associated with `b`."
hdf5file(b::MassCompatibleBufferedWriters) = file(b.filt_value.ds)
function write_to_hdf5(b::MassCompatibleBufferedWriters)
  write_to_hdf5.([b.filt_value, b.filt_phase, b.timestamp, b.rowcount, b.pretrig_mean, b.pretrig_rms, b.pulse_average,
  b.pulse_rms, b.rise_time, b.postpeak_deriv, b.peak_index, b.peak_value, b.min_value])
end
function Base.schedule(b::BufferedWriter)
  b.task=@async begin
    while !isready(b.endchannel)
      write_to_hdf5(b)
      sleep(b.timeout_s)
    end
    write_to_hdf5(b)
  end
end
stop(b::BufferedWriter) = put!(b.endchannel,true)
Base.fetch(b::BufferedWriter) = fetch(b.task)
Base.close(b::BufferedWriter) = (stop(b);fetch(b))
Base.flush(b::BufferedWriter) = flush(hdf5file(b))
function make_buffered_hdf5_writer(h5, channel_number, analzyer::MassCompatibleAnalysisFeb2017, chunksize=1000, timeout_s=1.0)
  g = g_require(h5,"chan$channel_number")
  d=MassCompatibleBufferedWriters(Channel{Bool}(1), timeout_s, Task(nothing),
   [BufferedHDF5Dataset{fieldtype(MassCompatibleDataProductFeb2017,name)}(g,
       string(name), chunksize) for name in mass_fieldnames]...)
  schedule(d)
  d
end

function Base.write(d::MassCompatibleBufferedWriters,x::MassCompatibleDataProductFeb2017)
  write(d.filt_value, x.filt_value)
  write(d.filt_phase, x.filt_phase)
  write(d.timestamp, x.timestamp)
  write(d.rowcount, x.rowcount)
  write(d.pretrig_mean, x.pretrig_mean)
  write(d.pretrig_rms, x.pretrig_rms)
  write(d.pulse_average, x.pulse_average)
  write(d.pulse_rms, x.pulse_rms)
  write(d.rise_time, x.rise_time)
  write(d.postpeak_deriv, x.postpeak_deriv)
  write(d.peak_index, x.peak_index)
  write(d.peak_value, x.peak_value)
  write(d.min_value, x.min_value)
end

function write_header(d::MassCompatibleBufferedWriters,r::LJHReaderFeb2017)
  channelgroup = parent(d.filt_value.ds)
  g_create(channelgroup, "calculated_cuts")
  channelgroup["calculated_cuts"]["pretrig_rms"]=r.analyzer.pretrigger_rms_cuts
  channelgroup["calculated_cuts"]["postpeak_deriv"]=r.analyzer.postpeak_deriv_cuts
  channelattrs = attrs(channelgroup)
  channelattrs["filename"]=r.fname
  channelattrs["pope_preknowledge_file"]=r.analyzer.pk_filename
  channelattrs["channum"]=LJH.channel(r.fname)
  channelattrs["noise_filename"]="analyzed by pope, see `pope_preknowledge_file`"
end
function write_header_end(d::MassCompatibleBufferedWriters,ljh,analyzer::MassCompatibleAnalysisFeb2017)
  # dont add datasets, groups or attributes after SWMR writing is started
  # channelgroup = parent(d.filt_value.ds)
  # channelattrs = attrs(channelgroup)
  # # when write_header_end is called, some values may not be flushed to hdf5 yet
  # # so I add lasti (number pulses written to hdf5) + length(v) (number pulses to be written)
  # channelattrs["npulses"]=d.filt_value.lasti+length(d.filt_value.v)
end
"""
    write_header_allchannel(d::MassCompatibleBufferedWriters, r::LJHReaderFeb2017)

Write "header" information assocaited with all channels to an hdf5 file.
This function is only called once per ljh file, not once per channel.
Also calls `start_swmr_write`. Remember all datasets, groups, and attributes must exist before
calling `start_swmr_write`.
"""
function write_header_allchannel(d::MassCompatibleBufferedWriters, r::LJHReaderFeb2017)
  h5 = hdf5file(d)
  a = attrs(h5)
  # execute only once for the whole HDF5 file
  "nsamples" in names(a) && error("only call_write_header_allchannels once per ljh file set, not once per channel")
  a["nsamples"]=r.analyzer.nsamples
  a["npresamples"]=r.analyzer.npresamples
  a["frametime"]=r.analyzer.frametime
  flush(h5)
  try
    HDF5.start_swmr_write(h5)
  catch
    println("SKIPPING START_SWMR_WRITE, HDF5 VERSION TOO LOW")
  end
end
# function (d::T)(x) where T<:BufferedWriter
#   write(d,x)
# end


text:
import Pkg
Pkg.add("HTTP")
Pkg.add("JSON3")
Pkg.add("PyPlot")

import HTTP
import JSON3

f = xx->sinc(4*xx)*xx^3 # oracle function to fit.

# generate data.
x_range = LinRange(-3, 5, 20)
y = f.(x_range)

# the query positions for plotting.
Nq = 1000
xq_range = LinRange(x_range[1], x_range[end], Nq)

# specify Gaussian process regression tuning parameters.
θ_len = 0.34
σ² = 1e-3
hyperparams = [θ_len; σ²]

# set up and send to microservice.
xq = collect(xq_range)
x = collect(x_range)

input_array = Vector{Vector{Float64}}(undef, 4)
input_array[1] = x
input_array[2] = y
input_array[3] = xq
input_array[4] = hyperparams

# put into the JSON format that the server script is expecting.
q =  """{"input_array":$(input_array)}"""

# z = HTTP.request("POST", "http://localhost:8000/jsonpayload",
#     [("Content-Type", "application/json")], q)

z = HTTP.request("POST", "http://rwresearch.ca/jsonpayload",
    [("Content-Type", "application/json")], q)

#out = convert(Vector{Float64}, z.body)
h = String(z.body)
out_j3 = JSON3.read(h)
yq = convert(Vector{Float64}, out_j3)

# if the server was successful in fitting a curve, i.e., the queried points yq is non-empty,
#   plot.
if !isempty(yq)
    import PyPlot
    PyPlot.close("all")
    fig_num = 1

    PyPlot.figure(fig_num)
    fig_num += 1

    PyPlot.plot(x, y, ".", label = "observed")
    PyPlot.plot(xq, yq, label = "GP mean function")

    title_string = "1D GP regression"
    PyPlot.title(title_string)
    PyPlot.legend()
else
    print("Bad input. The algorithm on the server failed to fit a Gaussian process regression model.")
end

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS