text:
add_verbose_scope(:Fields)
add_assert_scope(:Fields)

add_verbose_scope(:FieldsNonFancy)

mutable struct FieldsTower
  field::AnticNumberField
  generators_of_automorphisms::Vector{NfToNfMor}
  subfields::Vector{NfToNfMor}
  ramified_primes::Vector{fmpz}

  #Solvable embedding problems for the extension
  #They are here to improve the conductor computation
  imgs_autos::Vector{Main.ForeignGAP.MPtr}
  auts_for_conductors::Vector{Main.ForeignGAP.MPtr}
  proj_ext::Main.ForeignGAP.MPtr

  function FieldsTower(K::AnticNumberField, auts::Vector{NfToNfMor}, subfields::Vector{NfToNfMor})
    z = new()
    z.field = K
    z.generators_of_automorphisms = auts
    z.subfields = subfields
    return z
  end

end

include("./brauer.jl")
include("./merge.jl")
include("./abelian_layer.jl")
include("./read_write.jl")

Generic.degree(F::FieldsTower) = degree(F.field)
Hecke.maximal_order(F::FieldsTower) = maximal_order(F.field)

function ramified_primes(F::FieldsTower)
  if !isdefined(F, :ramified_primes)
    f = factor(discriminant(maximal_order(F.field)))
    F.ramified_primes = collect(keys(f.fac))
  end
  return F.ramified_primes
end

###############################################################################
#
#  From automorphisms to permutation groups
#
###############################################################################

function perm_grp(G::Vector{Hecke.NfRel_nsToNfRel_nsMor{nf_elem}})

  n = length(G)
  G1 = closure(G, *)
  Dc = Vector{Tuple{Hecke.NfRel_nsToNfRel_nsMor, Int}}(undef, length(G1))
  for i = 1:length(G1)
    Dc[i] = (G1[i], i)
  end
  D = Dict{Hecke.NfRel_nsToNfRel_nsMor, Int}(Dc)
  permutations = Array{Array{Int, 1},1}(undef, n)
  for s = 1:n
    perm = Array{Int, 1}(undef, length(G1))
    for i = 1:length(G1)
      a = G[s]*G1[i]
      perm[i] = D[a]
    end
    permutations[s] = perm
  end
  return _perm_to_gap_grp(permutations)
end

function perm_grp(G::Array{Hecke.NfToNfMor, 1})
  #TODO: I don't need the closure. I only need the images in the residue ring.
  K = domain(G[1])
  n = length(G)
  G1 = closure(G, degree(K))
  #First, find a good prime
  p = 11
  d = numerator(discriminant(K.pol))
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  pols = gfp_poly[Rx(g.prim_img) for g in G]
  for g in G1
    if g in G
      continue
    end
    push!(pols, Rx(g.prim_img))
  end
  Dcreation = Vector{Tuple{gfp_poly, Int}}(undef, length(pols))
  for i = 1:length(pols)
    Dcreation[i] = (pols[i], i)
  end
  D = Dict{gfp_poly, Int}(Dcreation)
  @assert length(D) == degree(K)
  permutations = Array{Array{Int, 1},1}(undef, n)
  for s = 1:n
    perm = Array{Int, 1}(undef, length(G1))
    for i = 1:length(G1)
      perm[i] = D[Hecke.compose_mod(pols[s], pols[i], fmod)]
    end
    permutations[s] = perm
  end
  return _perm_to_gap_grp(permutations)

end

function _from_autos_to_perm(G::Array{Hecke.NfToNfMor,1})

  K = domain(G[1])
  @assert degree(K) == length(G)
  n = length(G)
  #First, find a good prime
  p = 2
  d = numerator(discriminant(K.pol))
  while mod(d, p) == 0
    p = next_prime(p)
  end
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  fmod = Rx(K.pol)
  pols = Vector{Tuple{gfp_poly, Int}}(undef, n)
  for i = 1:n
    pols[i] = (Rx(G[i].prim_img), i)
  end
  D = Dict{gfp_poly, Int}(pols)
  permutations = Array{Array{Int, 1},1}(undef, n)
  for s = 1:n
    perm = Array{Int, 1}(undef, n)
    for i = 1:n
      perm[i] = D[Hecke.compose_mod(pols[i][1], pols[s][1], fmod)]
    end
    permutations[s] = perm
  end
  return permutations

end

function _perm_to_gap_grp(perm::Array{Array{Int, 1},1})
  g = Main.ForeignGAP.MPtr[]
  for x in perm
    z = _perm_to_gap_perm(x)
    push!(g, z)
  end
  A = GAP.Globals.Group(GAP.julia_to_gap(g))
  return A
end

function _perm_to_gap_perm(x::Array{Int, 1})
  z = GAP.Globals.PermList(GAP.julia_to_gap(x))
  return z
end

function IdGroup(autos::Array{NfToNfMor, 1})
  G = _from_autos_to_perm_grp(autos)
  return GAP.Globals.IdGroup(G)
end

################################################################################
#
#  final computation of the maximal order and automorphisms
#
################################################################################

function _from_relative_to_abs_with_embedding(L::Hecke.NfRel_ns{T}, autL::Array{Hecke.NfRel_nsToNfRel_nsMor{T}, 1}) where T

  S, mS = simple_extension(L)
  K, mK, MK = absolute_field(S, false)

  #First, we compute the maximal order of the absolute field.
  #We start from the maximal orders of the relative extension and of the base field.
  #FALSE: Since the computation of the relative maximal order is slow, I prefer to bring to the absolute field the elements
  # generating the equation order.
  pols = L.pol
  gL = gens(L)
  B = Array{nf_elem, 1}(undef, degree(K))
  B[1] = K(1)
  ind = total_degree(pols[1])
  genjj = mK\(mS\gL[1])
  for i = 2:ind
    B[i] = genjj*B[i-1]
  end
  for jj = 2:length(pols)
    genjj = mK\(mS\gL[jj])
    el = deepcopy(genjj)
    new_deg = total_degree(pols[jj])
    for i = 2:new_deg
      for j = 1:ind
        B[(i-1)* ind + j] = B[j]* el
      end
      mul!(el, el, genjj)
    end
    ind *= new_deg
  end

  #Now, I add the elements of the maximal order
  OB = maximal_order(base_field(S))
  for i = 1:degree(OB)
    el = MK(OB.basis_nf[i])
    for j = 1:ind
      B[(i-1)* ind + j] = B[j] * el
    end
  end
  @vprint :Fields 2 "Product basis computed\n"
  #Now, we compute the maximal order. Then we simplify.
  #We simplify only if the degree of the field is lower than 30

  BasisMat = basis_matrix(B, FakeFmpqMat)
  @vtime :Fields 3 Hecke.hnf_modular_eldiv!(BasisMat.num, BasisMat.den, :lowerleft)
  NewBMat = FakeFmpqMat(BasisMat.num, BasisMat.den)
  @vtime :Fields 3 Ostart = NfAbsOrd(K, NewBMat)
  Ostart.index = divexact(NewBMat.den^degree(K), prod(NewBMat.num[i, i] for i = 1:degree(K)))
  Ostart.gen_index = fmpq(Ostart.index)
  Ostart.disc = divexact(numerator(discriminant(K)), Ostart.index^2)
  ram_primes_rel = numerator(norm(discriminant(L)))
  ram_primes_down = Hecke.ramified_primes(OB)
  for p in ram_primes_down
    if isone(gcd(p, ram_primes_rel))
      push!(Ostart.primesofmaximality, p)
    end
  end
  @vtime :Fields 3 O1 = MaximalOrder(Ostart)
  O1.ismaximal = 1
  Hecke._set_maximal_order_of_nf(K, O1)

  @vtime :Fields 3 Ks, mKs = Hecke.simplify(K)
  #Now, we have to construct the maximal order of this field.
  #I am computing the preimages of mKs by hand, by inverting the matrix.
  @vtime :Fields 3 begin
  arr_prim_img = Array{nf_elem, 1}(undef, degree(Ks))
  arr_prim_img[1] = K(1)
  arr_prim_img[2] = mKs.prim_img
  for i = 3:degree(Ks)
    arr_prim_img[i] = arr_prim_img[i-1]*mKs.prim_img
  end
  M1 = inv(basis_matrix(arr_prim_img, FakeFmpqMat))
  end
  if isdefined(O1, :lllO)
    O2 = NfAbsOrd(Ks, basis_matrix(O1.lllO, copy = false)*M1)
    O2.lllO = O2
  else
    O2 = NfAbsOrd(Ks, basis_matrix(O1, copy = false)*M1)
  end
  O2.ismaximal = 1
  @assert isdefined(O1, :disc)
  O2.disc = O1.disc
  O2.index = root(divexact(numerator(discriminant(Ks)), O2.disc), 2)
  @vtime :Fields 3 OLLL = lll(O2)
  Hecke._set_maximal_order_of_nf(Ks, OLLL)
  #I want also the embedding of the old field into the new one.
  #It is enough to find the image of the primitive element.
  k = base_field(S)
  a = MK(gen(k))
  M = zero_matrix(FlintZZ, 1, degree(Ks))
  Hecke.elem_to_mat_row!(M, 1, denominator(a), a)
  mul!(M, M, M1.num)
  img_a = Hecke.elem_from_mat_row(Ks, M, 1, M1.den*denominator(a))
  embed = NfToNfMor(k, Ks, img_a)
  #@assert iszero(k.pol(img_a))
  @vprint :Fields 3 "MaximalOrder Computed. Now Automorphisms\n"
  #Now, the automorphisms.
  # I need both generators and the whole group.
  autos = Array{NfToNfMor, 1}(undef, length(autL))
  el = mS(mK(mKs.prim_img))
  el1 = mS(mK(gen(K)))
  for i=1:length(autL)
    #@assert iszero(K.pol(mK(mS\(autL[i](el1)))))
    x = mK\(mS\(autL[i](el)))
    Hecke.elem_to_mat_row!(M, 1, denominator(x), x)
    mul!(M, M, M1.num)
    y = Hecke.elem_from_mat_row(Ks, M, 1, M1.den*denominator(x) )
    #@assert Ks.pol(y) == 0
    autos[i] = Hecke.NfToNfMor(Ks, Ks, y)
  end
  #And the generators are done. Now the closure
  @vtime :Fields 3 Hecke._set_automorphisms_nf(Ks, closure(autos, degree(Ks)))
  #Hecke._set_automorphisms_nf(Ks, autsKs)
  @vprint :Fields 2 "Finished\n"
    #@assert codomain(embed) == Ks
  return Ks, autos, embed

end

###############################################################################
#
#  Split Extensions
#
###############################################################################

function _split_extension(G::Array{Hecke.NfToNfMor, 1}, mats::Array{Hecke.GrpAbFinGenMap, 1})

  gtype = map(Int, domain(mats[1]).snf)
  G1 = perm_grp(G)
  gensG1 = GAP.Globals.GeneratorsOfGroup(G1)
  A = GAP.Globals.AbelianGroup(GAP.julia_to_gap(gtype))
  gens = GAP.Globals.GeneratorsOfGroup(A)
  auts = Array{Main.ForeignGAP.MPtr, 1}(undef, length(mats))
  for i = 1:length(mats)
    images = Array{Main.ForeignGAP.MPtr, 1}(undef, length(gtype))
    for j = 1:length(gtype)
      g = GAP.Globals.Identity(A)
      for k = 1:length(gtype)
        if !iszero(mats[i].map[j,k])
          g *= gens[k]^Int(mats[i].map[j,k])
        end
      end
      images[j] = g
    end
    auts[i] = GAP.Globals.GroupHomomorphismByImages(A, A, gens, GAP.julia_to_gap(images))
  end
  AutGrp = GAP.Globals.Group(GAP.julia_to_gap(auts))
  mp = GAP.Globals.GroupHomomorphismByImages(G1, AutGrp, gensG1, GAP.julia_to_gap(auts))
  return GAP.Globals.SplitExtension(G1, mp, A)

end

###############################################################################
#
#  Check group extension
#
###############################################################################

function check_group_extension(TargetGroup::Main.ForeignGAP.MPtr, autos::Array{NfToNfMor, 1}, res_act::Array{GrpAbFinGenMap, 1})

  GS = domain(res_act[1])
  expo = Int(GS.snf[end])
  K = domain(autos[1])
  d = degree(K)
  com, uncom = ppio(expo, d)

  if com == 1
    # I only need to check the split extension, since the second cohomology group is
    # trivial, regardless of the action
    H = _split_extension(autos, res_act)
    return GAP.Globals.IdGroup(H) == TargetGroup
  end

  if uncom == 1
    #Need a cohomological check. Only useful in the prime power case.
    return true
  end

  # I check the split extension related to only uncom
  #Now, I have to check if the split extension is isomorphic to IdH
  Qn, mQn = quo(GS, uncom)
  S1, mS1 = snf(Qn)
  new_res_act = Array{GrpAbFinGenMap, 1}(undef, length(res_act))
  for i = 1:length(res_act)
    Mat = mS1.map*mQn.imap*res_act[i].map*mQn.map*mS1.imap
    Hecke.reduce_mod_snf!(Mat, S1.snf)
    new_res_act[i] = hom(S1, S1, Mat)
  end
  H = _split_extension(autos, new_res_act)
  return GAP.Globals.IdGroup(H) == TargetGroup

end

###############################################################################
#
#  Interface to find Fields
#
###############################################################################

function field_extensions(list::Vector{FieldsTower}, bound::fmpz, IsoE1::Main.ForeignGAP.MPtr, l::Array{Int, 1}, only_real::Bool)

  grp_to_be_checked = Dict{Int, Main.ForeignGAP.MPtr}()
  d = degree(list[1])
  n = lcm(l)
  com, uncom = ppio(n, d)
  fac = factor(n)
  for (p, v) in fac
    grp_to_be_checked[Int(p)] = _construct_grp(IsoE1, Int(p)^v)
  end
  if uncom != 1
    IsoCheck = _construct_grp(IsoE1, uncom)
  else
    IsoCheck = IsoE1
  end
  final_list = FieldsTower[]
  for (j, x) in enumerate(list)
    @vprint :Fields 1 "Field $(j)/$(length(list)): $(x.field.pol)"
    @vprint :FieldsNonFancy 1 "Field $(j)/$(length(list)): $(x.field.pol)\n"
    append!(final_list, field_extensions(x, bound, IsoCheck, l, only_real, grp_to_be_checked, IsoE1))
  end
  return final_list

end

debf = []

function field_extensions(x::FieldsTower, bound::fmpz, IsoE1::Main.ForeignGAP.MPtr, l::Array{Int, 1}, only_real::Bool, grp_to_be_checked::Dict{Int, Main.ForeignGAP.MPtr}, IsoG::Main.ForeignGAP.MPtr)

  list_cfields = _abelian_normal_extensions(x, l, bound, IsoE1, only_real, IsoG)
  if isempty(list_cfields)
    @vprint :Fields 1 "\e[1F$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())Number of new fields found: 0\n\n"
    @vprint :FieldsNonFancy 1 "Number of new fields found: 0\n"
    return Vector{FieldsTower}()
  end
  push!(debf, (list_cfields, x.generators_of_automorphisms, grp_to_be_checked))
  list = from_class_fields_to_fields(list_cfields, x.generators_of_automorphisms, grp_to_be_checked)
  @vprint :Fields 1 "Computing maximal orders"
  @vprint :FieldsNonFancy 1 "Computing maximal orders\n"
  final_list = Vector{FieldsTower}(undef, length(list))
  for j = 1:length(list)
    fld, autos, embed = _from_relative_to_abs_with_embedding(list[j][1], list[j][2])
    previous_fields = Array{NfToNfMor, 1}(undef, length(x.subfields)+1)
    for s = 1:length(x.subfields)
      previous_fields[s] = x.subfields[s]
    end
    previous_fields[end] = embed
    final_list[j] = FieldsTower(fld, autos, previous_fields)
  end
  @vprint :Fields 1 "$(Hecke.set_cursor_col())$(Hecke.clear_to_eol())"
  @vprint :Fields 1 "Number of new fields found: $(length(final_list))\n\n"
  @vprint :FieldsNonFancy 1 "Number of new fields found: $(length(final_list))\n"
  return final_list

end

###############################################################################
#
#  Interface
#
###############################################################################
function fields_direct_product(g1, g2, red, redfirst, absolute_bound; only_real = false)
  b1 = root(absolute_bound, g2[1])
  b2 = root(absolute_bound, g1[1])
  l2 = fields(g2[1], g2[2], b2, only_real = only_real)
  if isempty(l2)
    return FieldsTower[]
  end
  if g1 == g2
    return _merge(l2, l2, absolute_bound, red, redfirst)
  end
  l1 = fields(g1[1], g1[2], b1, only_real = only_real)
  if isempty(l1)
    return FieldsTower[]
  end
  return _merge(l1, l2, absolute_bound, red, redfirst)
end

function fields(a::Int, b::Int, absolute_bound::fmpz; using_direct_product::Bool = true, only_real::Bool = false)

  @assert absolute_bound > 0
  if a == 1
    @assert b == 1
    Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
    K, a = NumberField(x-1, cached = false)
    g = NfToNfMor(K, K, K(1))
    return FieldsTower[FieldsTower(K, NfToNfMor[g], Array{NfToNfMor, 1}())]
  end
  G = GAP.Globals.SmallGroup(a, b)
  @assert GAP.Globals.IsSolvable(G)
  if using_direct_product
    g1, g2, red, redfirst = direct_product_decomposition(G, (a, b))
    if g2 != (1, 1)
      return fields_direct_product(g1, g2, red, redfirst, absolute_bound; only_real = false)
    end
  end
  @vprint :Fields 1 "Doing Group ($a, $b) with bound $absolute_bound\n"
  @vprint :FieldsNonFancy 1 "Doing Group ($a, $b) with bound $absolute_bound\n"
  L = GAP.Globals.DerivedSeries(G)
  lvl = _real_level(L)
  #First step by hand
  l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(G))
  @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
  @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
  #@vtime :Fields 2
  list = abelian_extensionsQQ(l, root(absolute_bound, div(a,prod(l))), (lvl > 1 || only_real))
  if isempty(list)
    return list
  end
  @vprint :Fields 1 "First step completed\n \n"
  @vprint :FieldsNonFancy 1 "First step completed\n \n"
  for i = 2:length(L)-1
    H1 = GAP.Globals.FactorGroup(L[i], L[i+1])
    E1 = GAP.Globals.FactorGroup(L[1], L[i+1])
    l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(H1))
    @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
    @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
    o = divexact(a, GAP.Globals.Size(E1))
    bound = root(absolute_bound, o)
    IsoE1 = GAP.Globals.IdGroup(E1)
    @vprint :Fields 1 "Number of fields at the $i -th step: $(length(list)) \n"
    @vprint :FieldsNonFancy 1 "Number of fields at the $i -th step: $(length(list)) \n"
    lG = snf(DiagonalGroup(l))[1]
    invariants = map(Int, lG.snf)
    onlyreal = (lvl > i || only_real)
    #First, I search for obstruction.
    if iscyclic(lG)
      @vprint :Fields 1 "Computing obstructions\n"
      @vprint :FieldsNonFancy 1 "Computing obstructions\n"
      #@vtime :Fields 1
      list = check_Brauer_obstruction(list, L, i, invariants[1])
      @vprint :Fields 1 "Fields to check: $(length(list))\n\n"
      @vprint :FieldsNonFancy 1 "Fields to check: $(length(list))\n\n"
    end
    if isempty(list)
      return FieldsTower[]
    end
    list = field_extensions(list, bound, IsoE1, invariants, onlyreal)
    @vprint :Fields 1 "Step $i completed\n"
    @vprint :FieldsNonFancy 1 "Step $i completed\n"
    if isempty(list)
      return FieldsTower[]
    end
  end
  return list

end

function fields(list::Vector{FieldsTower}, G, absolute_bound::fmpz; only_real::Bool = false)
  L = GAP.Globals.DerivedSeries(G)
  lvl = _real_level(L)
  #First step by hand
  for i = 2:length(L)-1
    E1 = GAP.Globals.FactorGroup(L[1], L[i+1])
    if GAP.Globals.Size(E1) != degree(list[1].field)
      continue
    end
    H1 = GAP.Globals.FactorGroup(L[i], L[i+1])
    l = GAP.gap_to_julia(Vector{Int64}, GAP.Globals.AbelianInvariants(H1))
    @vprint :Fields 1 "contructing abelian extensions with invariants $l \n"
    @vprint :FieldsNonFancy 1 "contructing abelian extensions with invariants $l \n"
    o = divexact(GAP.Globals.Size(G), GAP.Globals.Size(E1))
    bound = root(absolute_bound, o)
    IsoE1 = GAP.Globals.IdGroup(E1)
    @vprint :Fields 1 "Number of fields at the $i -th step: $(length(list)) \n"
    @vprint :FieldsNonFancy 1 "Number of fields at the $i -th step: $(length(list)) \n"
    lG = snf(DiagonalGroup(l))[1]
    invariants = map(Int, lG.snf)
    onlyreal = (lvl > i || only_real)
    #First, I search for obstruction.
    if iscyclic(lG)
      @vprint :Fields 1 "Computing obstructions\n"
      @vprint :FieldsNonFancy 1 "Computing obstructions\n"
      #@vtime :Fields 1
      list = check_Brauer_obstruction(list, L, i, invariants[1])
      @vprint :Fields 1 "Fields to check: $(length(list))\n\n"
      @vprint :FieldsNonFancy 1 "Fields to check: $(length(list))\n\n"
    end
    if isempty(list)
      return FieldsTower[]
    end
    list = field_extensions(list, bound, IsoE1, invariants, onlyreal)
    @vprint :Fields 1 "Step $i completed\n"
    @vprint :FieldsNonFancy 1 "Step $i completed\n"
    if isempty(list)
      return FieldsTower[]
    end
  end
  return list
end


text:
"""
    get_path_to_example_config_files()::String

Returns the path to the example detector configuration files provided by the package.

See also [`SSD_examples`](@ref).
"""
function get_path_to_example_config_files()::String
    return joinpath(@__DIR__, "../examples/example_config_files")
end

"""
    SSD_examples::Dict{Symbol,String}

Dictionary with the paths to the example detector configuration files provided by the package.

Find the possible keys of the dictionary with `keys(SSD_examples)`.

The example detector configuration files can be loaded via
```
path_to_config_file = SSD_examples[:InvertedCoax]
sim = Simulation(path_to_config_file)
```
"""
const SSD_examples = Dict{Symbol,String}([
    :TwoSpheresCapacitor                     => joinpath(get_path_to_example_config_files(), "two_spheres_capacitor.yaml"),
    :InvertedCoax                            => joinpath(get_path_to_example_config_files(), "public_ivc_config.yaml"),
    :InvertedCoaxInCryostat                  => joinpath(get_path_to_example_config_files(), "ivc_splitted_config/public_ivc_cryostat_config.yaml"),
    :Coax                                    => joinpath(get_path_to_example_config_files(), "public_Coax_config.yaml"),
    :BEGe                                    => joinpath(get_path_to_example_config_files(), "public_SegBEGe_config.yaml"),
    :CGD                                     => joinpath(get_path_to_example_config_files(), "public_CGD_config.yaml"),
    :CGD_CylGrid                             => joinpath(get_path_to_example_config_files(), "public_CGD_config_cyl_grid.yaml"),
    :Spherical                               => joinpath(get_path_to_example_config_files(), "public_spherical_detector_config.yaml"),
    :SigGen                                  => joinpath(get_path_to_example_config_files(), "public_ppc_config_SigGen.config"),
    :InfiniteParallelPlateCapacitor          => joinpath(get_path_to_example_config_files(), "infinite_parallel_plate_capacitor.yaml"),
    :InfiniteCoaxialCapacitor                => joinpath(get_path_to_example_config_files(), "infinite_coaxial_capacitor.yaml"),
    :InfiniteCoaxialCapacitorCartesianCoords => joinpath(get_path_to_example_config_files(), "infinite_coaxial_capacitor_cartesian_coords.yaml"),
    :Hexagon                                 => joinpath(get_path_to_example_config_files(), "minimum_hexagon_config.yaml"),
    :CoaxialTorus                            => joinpath(get_path_to_example_config_files(), "public_coaxial_torus_config.yaml"),
])

ext: jl


text:
"""
    solve(FDDE::FDDEProblem, h, DelayABM())

Use the Adams-Bashforth-Moulton method to solve fractional delayed differential equations.

### References

```tex
@inproceedings{Bhalekar2011APS,
  title={A PREDICTOR-CORRECTOR SCHEME FOR SOLVING NONLINEAR DELAY DIFFERENTIAL EQUATIONS OF FRACTIONAL ORDER},
  author={Sachin Bhalekar and Varsha Daftardar-Gejji},
  year={2011}
}
```
"""
struct DelayABM <: FractionalDiffEqAlgorithm end
#FIXME: There are still some improvments about initial condition
#FIXME: Fix DelayABM method for FDDESystem : https://www.researchgate.net/publication/245538900_A_Predictor-Corrector_Scheme_For_Solving_Nonlinear_Delay_Differential_Equations_Of_Fractional_Order
#FIXME: Also the problem definition f(t, ϕ, y) or f(t, y, ϕ)?
function solve(FDDE::FDDEProblem, h, ::DelayABM)
    @unpack f, ϕ, α, τ, tspan = FDDE
    N::Int = round(Int, tspan/h)
    Ndelay::Int = round(Int, τ/h)
    x1 = zeros(Ndelay+N+1)
    x = zeros(Ndelay+N+1)
    #x1[Ndelay+N+1] = 0

    #x[Ndelay+N+1] = 0

    # History function handling
    #FIXME: When the value of history function ϕ is different with the initial value?
    if typeof(ϕ) <: Number
        x[1:Ndelay] = ϕ*ones(Ndelay)
    elseif typeof(ϕ) <: Function
        x[Ndelay] = ϕ(0)
        x[1:Ndelay-1] .= ϕ(collect(-h*(Ndelay-1):h:(-h)))
    end

    x0 = copy(x[Ndelay])

    x1[Ndelay+1] = x0 + h^α*f(0, x[1], x0)/(gamma(α)*α)

    x[Ndelay+1] = x0 + h^α*(f(0, x[1], x[Ndelay+1]) + α*f(0, x[1], x0))/gamma(α+2)

    @fastmath @inbounds @simd for n=1:N
        M1=(n^(α+1)-(n-α)*(n+1)^α)*f(0, x[1], x0)

        N1=((n+1)^α-n^α)*f(0, x[1], x0)

        @fastmath @inbounds @simd for j=1:n
            M1 = M1+((n-j+2)^(α+1)+(n-j)^(α+1)-2*(n-j+1)^(α+1))*f(0, x[j], x[Ndelay+j])
            N1 = N1+((n-j+1)^α-(n-j)^α)*f(0, x[j], x[Ndelay+j])
        end
        x1[Ndelay+n+1] = x0+h^α*N1/(gamma(α)*α)
        x[Ndelay+n+1] = x0+h^α*(f(0, x[n+1], x[Ndelay+n+1])+M1)/gamma(α+2)
    end

    xresult = zeros(N-Ndelay+1)
    yresult = zeros(N-Ndelay+1)

    xresult[N-Ndelay+1]=0
    yresult[N-Ndelay+1]=0

    @fastmath @inbounds @simd for n=2*Ndelay+1:N+Ndelay+1
       xresult[n-2*Ndelay] = x[n]
       yresult[n-2*Ndelay] = x[n-Ndelay]
    end

    return xresult, yresult
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build TestU01
sources = [
    "http://simul.iro.umontreal.ca/testu01/TestU01.zip" =>
    "bc1d1dd2aea7ed3b3d28eaad2c8ee55913f11ce67aec8fe4f643c1c0d2ed1cac",
    "./src",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd TestU01-1.2.3/
update_configure_scripts
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
  lt_cv_deplibs_check_method=pass_all ./configure --prefix=$prefix --host=$target LDFLAGS="-L/opt/$target/$target/lib/ -lws2_32"
else
  ./configure --prefix=$prefix --host=$target
fi
make tcode.o EXEEXT="" CC=/opt/x86_64-linux-gnu/bin/gcc
make tcode EXEEXT="" CC=/opt/x86_64-linux-gnu/bin/gcc LDFLAGS=""
make -j${nproc} EXEEXT=""
make install EXEEXT=""

# Compile TestU01extractors shim
cd $WORKSPACE/srcdir/TestU01extractors/
make install SHLIB_EXT=$dlext

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line

platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtestu01", :libtestu01),
    LibraryProduct(prefix, "libprobdist", :libprobdist),
    LibraryProduct(prefix, "libtestu01extractors", :libtestu01extractors),
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "TestU01", sources, script, platforms, products, dependencies)


text:
# Tests data_import.jl

using Test
using GeophysicalModelGenerator

# test the import of a csv file with depth given as positive values and units in km

# test loaded structure
# NOTE: currently deactivated these tests -----------------------------
#LoadedData   =   ReadCSV_LatLon("TestData.csv", "positive");
#@test LoadedData.lat.name  ==  "lat"
#@test LoadedData.lat.unit  ==  "deg"
#@test LoadedData.lat.values  ==  [36.424, 36.194, 36.144, 35.486, 36.106, 35.767, 36.24, 36.301, 36.401, 36.321]

#@test LoadedData.lon.name  ==  "lon"
#@test LoadedData.lon.unit  ==  "deg"
#@test LoadedData.lon.values  ==  [68.74, 71.111, 70.757, 70.38, 71.051, 70.55, 68.95, 71.303, 71.567, 70.829]

#@test LoadedData.depth.name  ==  "depth"
#@test LoadedData.depth.unit  ==  "km"
#@test LoadedData.depth.values  ==  [-52.04, -60.51, -60.74, -61.0, -75.5, -81.52, -82.55, -84.36, -90.9, -110.41]

#@test LoadedData.values.varnames == ["local_magnitude", "absError(km)"]
#@test LoadedData.values.vals == [5.5   9.7;5.6  10.0;5.1  10.7;5.5  20.4;5.1   6.5;5.9  10.1;5.2   5.5;5.1   3.2;5.3   3.6;5.3   3.6]
# ---------------------------------------------------------------------

# test loading images (profiles & mapviews)

# Extact & save profile in GeoData format
filename            =   "test.png";             # fake png
Corner_LowerLeft    =   (18.0, 51.0, -590.0)
Corner_UpperRight   =   (9.0, 42.0,    0.0)
data_Image          =   Screenshot_To_GeoData(filename,Corner_LowerLeft, Corner_UpperRight)
@test data_Image.lon[1000]==17.592964824120603
@test data_Image.lat[1000]==50.59296482412061
@test data_Image.depth[1000]==-590km
@test Write_Paraview(data_Image, "Profile_1")==nothing

# Test in CartData
data_Image          =   Screenshot_To_GeoData(filename,Corner_LowerLeft, Corner_UpperRight, Cartesian=true)
@test data_Image.x.val[22] == 18.0km
@test data_Image.y.val[22] == 51.0km
@test data_Image.z.val[22] ≈ -125.15151515151516km

# Test in UTM zone [note that depth should be in m]
data_Image          =   Screenshot_To_GeoData(filename,Corner_LowerLeft, Corner_UpperRight, UTM=true, UTMzone=33, isnorth=true)
@test data_Image.EW.val[22] == 18.0
@test data_Image.NS.val[22] == 51.0
@test data_Image.depth.val[22] ≈ -125.15151515151516m

# Mapview (distorted) in GeoData format
filename            =   "test.png";             # fake png
Corner_LowerLeft    =   (2.0,  40.0, -15.0)
Corner_UpperRight   =   (22.0, 51.0, -15.0)
Corner_LowerRight   =   (20.0, 40.0, -15.0)
Corner_UpperLeft    =   (0.0,  51.0, -15.0)
data_Image          =   Screenshot_To_GeoData(filename,Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight, Corner_UpperLeft=Corner_UpperLeft)
@test data_Image.lon[1000]==2.814070351758794
@test data_Image.lat[1000]==40.00000000000001
@test data_Image.depth[1000]==-15km
@test Write_Paraview(data_Image, "MapView_1") == nothing

# MapView in CartData
data_Image          =   Screenshot_To_CartData(filename,Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight, Corner_UpperLeft=Corner_UpperLeft)
@test data_Image.x.val[22] ≈ 0.42424242424242425km
@test data_Image.y.val[22] ≈ 48.666666666666664km
@test data_Image.z.val[22] ≈ -15km

# MapView in UTMData
data_Image          =   Screenshot_To_UTMData(filename,Corner_LowerLeft, Corner_UpperRight, Corner_LowerRight=Corner_LowerRight, Corner_UpperLeft=Corner_UpperLeft, UTMzone=33, isnorth=true)
@test data_Image.EW.val[22] ≈ 0.42424242424242425
@test data_Image.NS.val[22] ≈ 48.666666666666664
@test data_Image.depth.val[22] ≈ -15.0m


text:
# # Example 5 : Array Types
#
# A `ClimArray` (struct from `ClimateBase.jl`), which includes metadata read from nectdf file, is written back to a netcdf file via `NCTiles.jl`. Then the reverse is done to illustrate consistent workflows.

using ClimateBase, NCDatasets, NCTiles, Dates, Unitful

# File Paths
inputs=NCTiles.NCTILES_TESTCASES
NCTiles.ensure_testcases_installed()

outputs = joinpath(tempdir(),"NCTILES_TESTCASES_OUTPUT/")
if ~ispath(outputs); mkpath(outputs); end

savedir = joinpath(outputs,"ex5")
if ~ispath(savedir); mkpath(savedir); end

# ## 1. Read via NCTiles.jl and write via ClimateBase.jl

# ### Helper function

"""
        NCvar_to_ClimArray(ncvar::NCvar,ncdims::Dict)

Creates a ClimArray struct, from a NCvar struct, which can then be written
to a NetCDF file using `ClimateBase.ncwrite()`.
"""
function NCvar_to_ClimArray(ncvar::NCvar,ncdims::Dict)
    lons=ncdims["lon_c"].values[:]
    lats=ncdims["lat_c"].values[:]
    t=ncdims["tim"].values[:]
    dimensions = (Lon(lons), Lat(lats), ClimateBase.Ti(t))

    name=ncvar.name
    data=ncvar.values[:]
    units=ncvar.units
    return ClimArray(data, dimensions; name = name, attrib = Dict("units" => units))
end

# ### Main call sequence :

fil=joinpath(outputs,"ex1/ex1_NetCDF.nc")
ncvars,ncdims,fileatts = NCTiles.readncfile(fil)

#filout=joinpath(savedir,"ex5_NCTiles.nc")
#write(ncvars["ETAN"],filout,globalattribs=fileatts)

filout=joinpath(savedir,"ex5_ClimateBase.nc")
A=NCvar_to_ClimArray(ncvars["ETAN"],ncdims)
ClimateBase.ncwrite(filout, A)

# ## 2. Read via ClimateBase.jl and write via NCTiles.jl
#
# ### Helper function

"""
        ClimArray_to_NCvar(C::ClimGrid,N::String)

Creates an NCvar struct, from a ClimArray struct, which can then be written
to a NetCDF file using `NCTiles.write()`.

```
C = ClimateBase.ncread(fil, "ETAN")
writefld = ClimArray_to_NCvar(C,"ETAN")
NCTiles.write(writefld,"myfile.nc")
```
"""
function ClimArray_to_NCvar(C::ClimArray,N::String)
        x = C.dims[1][:]
        y = C.dims[2][:]
        timevec = DateTime.(C.dims[3][:])
        timeunit = "days since 1992-01-01"
        timevec = NCDatasets.timeencode(timevec, timeunit)

        dims = [NCvar("lon","degrees_east",size(C.data)[1],x,Dict("long_name" => "longitude"),NCDatasets),
                NCvar("lat","degrees_north",size(C.data)[2],y,Dict("long_name" => "latitude"),NCDatasets),
                NCvar("time",timeunit,Inf,timevec,Dict(("long_name" => "Ti","standard_name" => "time")),NCDatasets)
                ]

        return NCvar(N,C.attrib["units"],dims,C.data,C.attrib,NCDatasets)
end

# ## Main call sequence :
#
# _Note: `ClimateBase.jl` relies on specific dimension names to identify space and time dimensions (as follows). This sometimes lead to warnings when files use different names._
#
# - latitude: lat, latitude, rlat, y, yc
# - longitude: lon, longitude, rlon, x, xc
# - time: time

ETAN = ClimateBase.ncread(fil, "ETAN")
ETAN = ClimArray_to_NCvar(ETAN,"ETAN")
filout=joinpath(savedir,"ex5_ClimArray_to_NCvar.nc")
NCTiles.write(ETAN,filout)

##

using MeshArrays

function ClimArray_to_MeshArray(C::ClimArray)
        u=uparse(C.attrib["units"])
        n=string(C.name)
        ln=C.attrib["long_name"]
        tim=DateTime.(collect(tmp.dims[3][:]))
        m=varmeta(u,fill(0.5,3),tim,n,ln)
        #MeshArray(C.data;meta=m)

        nlon=length(tmp.dims[1][:])
        nlat=length(tmp.dims[2][:])
        XC = MeshArray(tmp.dims[1][:]*ones(1,nlon))
        YC = MeshArray(ones(nlat,1)*tmp.dims[2][:]')
        Γ = (XC=XC,YC=YC)

        MeshArray(C.data;meta=m),Γ
end

function NCvar_to_MeshArray(ncvar::NCvar)
        u=uparse(ncvar.units)
        n=string(ncvar.name)
        ln=ncvar.atts["long_name"]
        tim=ncvar.dims[3].values[:]
        m=varmeta(u,fill(0.5,3),tim,n,ln)
        #MeshArray(ncvar.values[:];meta=m)

        nlon=length(ncvar.dims[1].values[:])
        nlat=length(ncvar.dims[2].values[:])
        XC = MeshArray(ncvar.dims[1].values[:]*ones(1,nlon))
        YC = MeshArray(ones(nlat,1)*ncvar.dims[2].values[:]')
        Γ = (XC=XC,YC=YC)

        MeshArray(ncvar.values[:];meta=m),Γ
end

fil=joinpath(outputs,"ex1/ex1_NetCDF.nc")
tmp = ClimateBase.ncread(fil, "ETAN")
ETAN,Γ = ClimArray_to_MeshArray(tmp)

ncvars,ncdims,fileatts = NCTiles.readncfile(fil)
ETAN,Γ = NCvar_to_MeshArray(ncvars["ETAN"])


text:
export State, GraphModel, ControlAgentState, NetworkAssetState, ModelState, Flow, Flow_Type, MRule, OFS_Action
"""
Self-organising agent
"""

abstract type SOAgent <: AbstractAgent end
abstract type SimAsset <: AbstractAgent end
abstract type State end
abstract type Packet end

abstract type MaintenanceType end
abstract type CorrectiveM <: MaintenanceType end
abstract type PreventiveM <: MaintenanceType end
abstract type PredictiveM <: MaintenanceType end

"""
    "Real" state of the asset
"""
mutable struct SimpleAssetState <: State
    color::Symbol
    condition_trj::Array{Float64,2}
    queue::Channel{Packet}
end

"""
    Agent's state, including observed state of
    monitoring/controlling asset
"""
mutable struct SimpleAgState <: State
    color::Symbol
    condition_trj::Array{Float64,2}
    health_trj::Vector{Float64}
end

#fireflystate
mutable struct FFAgState <: State
    color::Symbol
    #T::Int64 # period to fire
    phase::Float64 #internal clock
    #inbox::Vector{Tuple{Int64,Float64}}
    inbox::Vector{Int64}
    blind::Int64
    maxN::Int64
    maxId::Int64
    condition_trj::Array{Float64,2}
    health_trj::Vector{Float64}
    #pulse_strategy::Symbol # FIFO, NEAR, MIXED
    #ϕ
end

"""
Supported Graphs
"""
@enum GraphModel begin
    CUSTOM=0
    CENTRALISED=1 # Only for control model
    RING=2
    COMPLETE=3
    GRID=4
    STAR=5
    BA_RANDOM=6 # Barrabasi_Albert
    WS_RANDOM=7 # watts_strogatz
    #SM_RANDOM_TOPO=8 # Stochastic Block Model
end

"""
Indicates the type of a flow within a path
"""
@enum Flow_Type begin
    f_SE = -2 # start and end
    f_S = -1 # start
    f_I = 0 # intermediate
    f_E = 1 # end
end

mutable struct ModelState <: State
    tick::Int
    links_load::Dict{Tuple{Int,Int},Int} # key: (src,dst), value: pkts
    active_flows::Vector{Tuple{Int,Int,Flow_Type}} # source,destination,type
end

@enum CTL_Event begin
    NetworkChange = 1 #Network change and reroute traffic
    TriggerDown = 2 #Trigger SNE down
    TriggerUp = 3 #Trigger SNE up
    Prediction = 4 #Update RUL predictions
    InstallFlows = 5 #Install flows already calculated, different from 1, which calculates flow based on network shortest path
end

@enum Ofp_Event begin
    EventOFPPortStatusDown=1
    EventOFPPortStatusUp=2
end
@enum Ofp_Protocol begin
    OFPR_ACTION = 1
    OFPPR_DELETE = 2
    OFPR_ADD_FLOW = 3
    OFPR_NO_MATCH = 4
    OFPPR_JOIN = 5
end
@enum OFS_Action begin
    OFS_Output = 1
    OFS_Drop = 2
end

@enum AG_Protocol begin
    QUERY_PATH = 1
    MATCH_PATH = 2
    NEW_NB = 3
    NE_DOWN = 4
    PREDICTED_NES_DOWN = 5
end

#if has come from in_port and src, going to dst
mutable struct MRule
    in_port::String
    src::String
    dst::String
end

mutable struct DPacket <: Packet
    id::Int64
    src::Int64
    dst::Int64
    size::Float64 # in bytes
    time_sent::Int64
    hop_limit::Int64
end

abstract type CTLMessage end

mutable struct OFMessage <: CTLMessage
    id::Int64
    ticks::Int # ticks when msg was created. Note on forwarding new msgs are created.
    dpid::Int # sender of msg,  aka SimNE.id aka switch.id
    in_port::Int # (Optional) sender's input port, in SimNE to ControlAg comm, this is the port where the packet was received
    reason::Ofp_Protocol
    data::Any
end

mutable struct AGMessage <: CTLMessage
    id::Int
    ticks::Int
    sid::Int # sender id
    rid::Int # receiver id
    reason::AG_Protocol # Type of msg as per enum
    body::Dict{Symbol,Any}
end

mutable struct OFEvent
    msg::OFMessage
end

mutable struct ControlEvent
    #When event is triggered
    time::Int64
    #type of event
    type::CTL_Event
    #affected snes
    snes::Vector{Int64}
end

mutable struct Flow
    dpid::Int64 # datapath id
    match_rule::MRule
    params::Vector{Any}
    action::OFS_Action # out_port, TODO: Check if other actions are needed
end

#State for which we need to save trajectory
mutable struct NetworkAssetState <: State
    ne_id::Int64
    up::Bool
    on_maintenance::Bool
    port_edge_list::Vector{Tuple{Int64,String}}
    in_pkt::Int64
    q_size::Int64
    out_pkt::Int64
    drop_pkt::Int64
    flow_table::Vector{Flow}
    throughput_out::Dict{Int64,Float64} # link/port, value
    # throughput_in::Float64
    condition_ts::Array{Float64,1} # sensor data related to the condition of the asset
    rul::Float64 # "real" rul
    rul_e::Float64 # estimated rul
    maintenance_due::Int64 #
end

mutable struct ControlAgentState <: State
    a_id::Int64
    up::Bool
    active_paths::Dict{Tuple{Int64,Int64},Array{Int64}}
    in_ag_msg::Float64
    out_ag_msg::Float64
    in_of_msg::Float64
    out_of_msg::Float64
    q_queries::Float64
    path_scores::Array{Tuple{Int64,Int64,Float64}}
end

"""
 It represents data structure for maintenance related info
"""
mutable struct MaintenanceInfo
    #Current maintenance policy
    policy::Type{<:MaintenanceType}
    # Expected useful life. Standard manufacturer time-to-failure
    eul::Int64
    # When last job started
    job_start::Int64
    # duration of current/next job
    duration::Int64
    # Cost of this type of maintenance
    cost::Float64
    #To trigger preventive maintenance ahead of breakdown
    threshold::Int64
    #How often rul predictions are run
    predictive_freq::Int64
    # how many steps ahead the prediction is going to be for
    prediction_window::Int64
    deterioration_parameter::Float64
    #Duration of the corrective maintenance
    reference_duration::Int64
    #Cost of the corrective maintenance
    reference_cost::Float64
end

"""
    Control Agent
"""
mutable struct Agent <: SOAgent
    id::Int64
    pos::Int64
    color::Symbol
    size::Float16
    pending::Vector{Tuple{Int64,OFMessage,Bool}} # <-time initially processed, ofmsg, reprocess in next tick?) reprocess only if a match is done
    # key(src,dst): value(tick_found,confidence,score,path)
    paths::Dict{Tuple{Int64,Int64},Array{Tuple{Int64,Float64,Float64,Array{Int64}}}} # pre-calculated paths for operation
    state_trj::Vector{ControlAgentState}
    msgs_links::Array{Vector{AGMessage},2}
    msgs_in::Vector{AGMessage}
    queue::Channel{OFMessage}
    previous_queries::Dict{Tuple{Int64,Int64},Tuple{Int64,Array{Int64}}} # (src,dst):(tick last queried,[nb ag queried])
    matched_queries::Dict{Tuple{Int64,Int64,Array{Int64}},Int64}
    ctl_paths::Vector{Array{Int64}}
    #Maintenance info for all assets controlled by this agent
    maintenance::MaintenanceInfo
    rul_predictions::Matrix{Float64}
    events::Dict{Int64,Vector{ControlEvent}}
    ntw_graph::MetaGraph
    base_ntw_graph::MetaGraph
    ctl_graph::MetaGraph
    last_cache_cleared::Int64
    params::Dict{Symbol,Any}
end

function Agent(id::Int64,nid::Int64,maintenance::MaintenanceInfo,rul_predictions,params::Dict{Symbol,Any})::Agent
    s0 = ControlAgentState(id,true,Dict(),0,0,0,0,0,[])
    Agent(id,nid,:lightblue,0.1,Vector{Tuple{Int64,OFMessage,Bool}}(),Dict{Tuple{Int64,Int64},Array{Tuple{Int64,Float64,Array{Int64}}}}(),[s0],Array{Vector{AGMessage}}(undef,1,1),[],Channel{OFMessage}(500),Dict(),Dict(),[],maintenance,rul_predictions,Dict{Int64,Vector{ControlEvent}}(),MetaGraph(),MetaGraph(),MetaGraph(),0,params)
end

"""
    Simulated Physical Network Element
"""
mutable struct SimNE <: SimAsset
    id::Int64
    pos::Int64
    size::Float16
    queue::Channel{OFMessage} #
    pending::Vector{OFMessage}
    requested_ctl::Dict{Tuple{Int64,Int64},Int64} # flows requested to controller: key{src,dst}:value{tick}
    state_trj::Vector{NetworkAssetState}
    one_way_time_pkt::Dict{Int64,Array{Int64}}
    controller_id::Int64
    #Individual Maitenance info for this asset
    maintenance::MaintenanceInfo
    params::Dict{Symbol,Any}
end

function SimNE(id,nid,params,max_q,maintenance)
    SimNE(id,nid,0.3,Channel{OFMessage}(max_q),Vector{OFMessage}(),Dict{Tuple{Int64,Int64},Int64}(),[NetworkAssetState(id)],Dict(),-1,maintenance,params) #initialise SimNE with a placeholder in the controller
end

function ModelState(tick::Int)
    ModelState(tick,Dict(),[])
end

function SimpleAgState(condition_trj::Array{Float64,2}, health_trj::Vector{Float64})
    SimpleAgState(:red,condition_trj,health_trj)
end

function OFMessage(id::Int64,ticks::Int,dpid::Int,reason::Ofp_Protocol)
    return OFMessage(id,ticks,dpid,-1,reason,nothing)
end

"""
    Message with default reason: Forward
"""
function OFMessage(id::Int64,ticks::Int,dpid::Int,in_port::Int,data::DPacket)
    return OFMessage(id,ticks,dpid,in_port,OFPR_ACTION,data)
end

"""
    Message with no input port
"""
function OFMessage(id::Int64,ticks::Int,dpid::Int,ofp::Ofp_Protocol,data::Int)
    return OFMessage(id,ticks,dpid,-1,ofp,data)
end

"""
    Message with pair in data
"""
function OFMessage(id::Int64,ticks::Int,dpid::Int,ofp::Ofp_Protocol,data::Tuple{Int64,Int64})
    return OFMessage(id,ticks,dpid,-1,ofp,data)
end

function NetworkAssetState(ne_id::Int)
    NetworkAssetState(ne_id,true,false,Vector{Tuple{Int64,String}}(),0,0,0,0,Vector{Flow}(),Dict(),Array{Float64,1}(),0.0,0.0,0)
end


text:
module Cytof5
#= Note:
Julia uses Gamma(shape, scale) and InverseGamma(shape, scale).
Note that InverseGamma(shape, scale) IS my preferred parameterization.
It has a mean of scale / (shape - 1) for shape > 1.
=#

export MCMC

include("Util/Util.jl")
include("MCMC/MCMC.jl")
include("Model/Model.jl")
include("VB/VB.jl")

# TODO: get rid of RCall, use Plots with GR backend.
# include("Plot/Plot.jl")

# For R plots
include("R/Plots.jl")

end # module


text:
"""
    function mincost_flow(g::AbstractGraph,
    		node_demand::AbstractVector,
    		edge_capacity::AbstractMatrix,
    		edge_cost::AbstractMatrix,
    		optimizer;
    		edge_demand::Union{Nothing,AbstractMatrix} = nothing,
    		source_nodes = (),
    		sink_nodes = ()
    		)

Find a flow over a directed graph `g` satisfying the `node_demand` at each
node and `edge_capacity` constraints for each edge while minimizing `dot(edge_cost, flow)`.

Returns a sparse flow matrix, flow[i,j] corresponds to the flow on the (i,j) arc.

# Arguments

- `g` is a directed `LightGraphs.AbstractGraph`.
-`node_demand` is a vector of nodal demand values, which should be negative for sources,
positive for sink nodes, and zero for all other nodes.
- `edge_capacity::AbstractMatrix` sets an upper bound on the flow of each arc.
- `edge_cost::AbstractMatrix` the cost per unit of flow on each arc.
- `optimizer` is an optimizer constructor, like `Clp.Optimizer` or `() -> Clp.Optimizer()` passed at the construction of the `JuMP.Model`.

# Keyword arguments

- `edge_demand::Union{Nothing,AbstractMatrix}`: require a minimum flow for edges, or nothing.
- `source_nodes` Collection of sources at which the nodal netflow is allowed to be greater than nodal demand, defaults to an empty tuple.
- `sink_nodes` Collection of sinks at which the nodal netflow is allowed to be less than nodal demand, defaults to an empty tuple.

`source_nodes` & `sink_nodes` are only needed when nodal flow are not explictly set in node_demand

### Usage Example:

```julia
julia> import LightGraphs
julia> const LG = LightGraphs
julia> using LightGraphsFlows: mincost_flow
julia> import Clp # use your favorite LP solver here
julia> using SparseArrays: spzeros
julia> g = LG.DiGraph(6) # Create a flow-graph
julia> LG.add_edge!(g, 5, 1)
julia> LG.add_edge!(g, 5, 2)
julia> LG.add_edge!(g, 3, 6)
julia> LG.add_edge!(g, 4, 6)
julia> LG.add_edge!(g, 1, 3)
julia> LG.add_edge!(g, 1, 4)
julia> LG.add_edge!(g, 2, 3)
julia> LG.add_edge!(g, 2, 4)
julia> cost = zeros(6,6)
julia> cost[1,3] = 10
julia> cost[1,4] = 5
julia> cost[2,3] = 2
julia> cost[2,4] = 2
julia> demand = spzeros(6)
julia> demand[5] = -2
julia> demand[6] = 2
julia> capacity = ones(6,6)
julia> flow = mincost_flow(g, demand, capacity, cost, Clp.Optimizer)
```
"""
function mincost_flow end

@traitfn function mincost_flow(g::AG::lg.IsDirected,
		node_demand::AbstractVector,
		edge_capacity::AbstractMatrix,
		edge_cost::AbstractMatrix,
		optimizer;
		edge_demand::Union{Nothing,AbstractMatrix} = nothing,
		source_nodes = (), # Source nodes at which to allow a netflow greater than nodal demand
		sink_nodes = ()	   # Sink nodes at which to allow a netflow less than nodal demand
		) where {AG <: lg.AbstractGraph}

	m = JuMP.Model(optimizer)
	vtxs = vertices(g)

	source_nodes = [v for v in vtxs if v in source_nodes || node_demand[v] < 0]
	sink_nodes = [v for v in vtxs if v in sink_nodes || node_demand[v] > 0]

	@variable(m, 0 <= f[i=vtxs,j=vtxs; (i,j) in lg.edges(g)] <= edge_capacity[i, j])
	@objective(m, Min, sum(f[src(e),dst(e)] * edge_cost[src(e), dst(e)] for e in lg.edges(g)))

	for v in lg.vertices(g)
	    if v in source_nodes
            @constraint(m,
                sum(f[v, vout] for vout in outneighbors(g, v)) - sum(f[vin, v] for vin in lg.inneighbors(g, v)) >= -node_demand[v]
            )
	    elseif v in sink_nodes
            @constraint(m,
                sum(f[vin, v] for vin in lg.inneighbors(g, v)) - sum(f[v, vout] for vout in outneighbors(g, v)) >= node_demand[v]
            )
	    else
            @constraint(m,
                sum(f[vin, v] for vin in lg.inneighbors(g, v)) == sum(f[v, vout] for vout in outneighbors(g, v))
            )
        end
	end

    if edge_demand isa AbstractMatrix
        for e in lg.edges(g)
            (i,j) = Tuple(e)
            JuMP.set_lower_bound(f[i,j], edge_demand[i,j])
        end
    end
    optimize!(m)
    ts = termination_status(m)
    result_flow = spzeros(nv(g), nv(g))
    if ts != MOI.OPTIMAL
        @warn "Problem does not have an optimal solution, status: $(ts)"
        return result_flow
    end
    for e in lg.edges(g)
        (i,j) = Tuple(e)
        result_flow[i,j] = JuMP.value(f[i,j])
    end
    return result_flow
end


text:
import Pkg
Pkg.instantiate()
using Documenter, FlatBuffers

makedocs(
    modules = [FlatBuffers],
    format = :html,
    sitename = "FlatBuffers.jl",
    pages = ["Home" => "index.md"]
)

deploydocs(
    repo = "github.com/JuliaData/FlatBuffers.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module Poisson

using IterativeSolvers: cg
using AlgebraicMultigrid

using STLCutters
using Gridap
using GridapEmbedded
using Test

function main(filename;n=20,δ=0.2,output=nothing)

  # Manufactured solution
  u(x) = x[1] + x[2] - x[3]
  f(x) = - Δ(u)(x)
  ud(x) = u(x)

  geo = STLGeometry( filename )

  n = 20
  δ = 0.2
  pmin,pmax = get_bounding_box(geo)
  diagonal = pmax-pmin
  pmin = pmin - diagonal*δ
  pmax = pmax + diagonal*δ
  partition = (n,n,n)

  bgmodel = CartesianDiscreteModel(pmin,pmax,partition)

  # Cut the background model

  cutgeo,facet_to_inoutcut = cut(bgmodel,geo)

  strategy = AggregateAllCutCells()
  aggregates = aggregate(strategy,cutgeo,facet_to_inoutcut)

  # Setup integration meshes
  Ω_bg = Triangulation(bgmodel)
  Ω = Triangulation(cutgeo,PHYSICAL)
  Γd = EmbeddedBoundary(cutgeo)

  # Setup normal vectors
  n_Γd = get_normal_vector(Γd)

  #writevtk(Ω,"trian_O")
  #writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
  #writevtk(Triangulation(bgmodel),"bgtrian")

  # Setup Lebesgue measures
  order = 1
  degree = 2*order
  dΩ = Measure(Ω,degree)
  dΓd = Measure(Γd,degree)

  # Setup FESpace
  Ω_act = Triangulation(cutgeo,ACTIVE)
  Vstd = FESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)

  V = AgFEMSpace(Vstd,aggregates)
  U = TrialFESpace(V)

  # Weak form
  γd = 10.0
  h = (pmax - pmin)[1] / partition[1]

  a(u,v) =
    ∫( ∇(v)⋅∇(u) ) * dΩ  +
    ∫( (γd/h)*v*u  - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u ) * dΓd

  l(v) =
    ∫( v*f ) * dΩ +
    ∫( (γd/h)*v*ud - (n_Γd⋅∇(v))*ud ) * dΓd

  # FE problem
  op = AffineFEOperator(a,l,U,V)

  A = get_matrix(op)
  b = get_vector(op)

  p = aspreconditioner(smoothed_aggregation(A))

  x = cg(A,b,verbose=true,Pl=p,reltol=1e-12)

  uh = FEFunction(U,x)

  e = u - uh

  # Postprocess
  l2(u) = sqrt(sum( ∫( u*u )*dΩ ))
  h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ ))

  el2 = l2(e)
  eh1 = h1(e)
  ul2 = l2(uh)
  uh1 = h1(uh)

  if output !== nothing
    colors = color_aggregates(aggregates,bgmodel)
    writevtk(Ω_bg,output*"_trian",celldata=["aggregate"=>aggregates,"color"=>colors],cellfields=["uh"=>uh])
    writevtk(Ω,output*"_results",cellfields=["uh"=>uh])
  end

  @test el2/ul2 < 1.e-9
  @test eh1/uh1 < 1.e-9
end

end # module


text:
function main()
    N = parse(Int, readline())
    AB = [parse.(Int, split(readline())) for _ in 1:N]
    AB_s = sort(AB, by = x -> x[2])

    ans = "Yes"
    tmp = 0
    for ab in AB_s
        tmp += ab[1]
        if tmp > ab[2]
            ans = "No"
            break
        end
    end

    println(ans)
end

main()


text:
# load requried packages
using IncrementalInference

## parameters
lm_prior_noise = 0.01
meas_noise = 0.25
odom_noise = 0.1
n_samples = 100

# initialize mean landmark locations
l0 = 0.0
l1 = 10.0
l2 = 40.0

# "Ground-truth" robot poses
x0 = 0.0
x1 = 10.0
x2 = 20.0
x3 = 40.0

## Initialize empty factor graph
fg = initfg()

# Place strong prior on locations of three "doors"
addVariable!(fg, :l0, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l0], Prior(Normal(l0, lm_prior_noise)))

addVariable!(fg, :l1, ContinuousScalar, N=n_samples)
addFactor!(fg, [:l1], Prior(Normal(l1, lm_prior_noise)))

# Add first pose
addVariable!(fg, :x0, ContinuousScalar, N=n_samples)

# Make first "door" measurement
# addFactor!(fg, [:x0; :l0], LinearRelative(Normal(0, meas_noise)))
addFactor!(fg, [:x0; :l0; :l1], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/2.0; 1.0/2.0])

# Add second pose
addVariable!(fg, :x1, ContinuousScalar, N=n_samples)

# Gaussian transition model
addFactor!(fg, [:x0; :x1], LinearRelative(Normal(x1-x0, odom_noise)))

# Make second "door" measurement
# addFactor!(fg, [:x1; :l1], LinearRelative(Normal(0, meas_noise)) )
addFactor!(fg, [:x1; :l0; :l1], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/2.0; 1.0/2.0])

## Add one more pose/odometry to invoke issue #236

# Add third pose
addVariable!(fg, :x2, ContinuousScalar, N=n_samples)
addFactor!(fg, [:x1; :x2], LinearRelative(Normal(x2-x1, odom_noise)))

# Add fourth pose
# addVariable!(fg, :x3, ContinuousScalar, N=n_samples)

# Add odometry transition and new landmark sighting
# addFactor!(fg, [:x2, :x3], LinearRelative(Normal(2, odom_noise)))
# addFactor!(fg, [:x3; :l0; :l1; :l2], LinearRelative(Normal(0, meas_noise)), multihypo=[1.0; 1.0/3.0; 1.0/3.0; 1.0/3.0])

## Do some debugging
initAll!(fg)

##

drawGraph(fg, show=true)

tree = buildTreeReset!(fg, drawpdf=true, show=true)

## Solve graph
tree = solveTree!(fg)

# tree = buildTreeReset!(fg, drawpdf=true, show=true)

## Plotting functions below

using RoMEPlotting

pl = plotKDE(fg, [:x0;:x1])

pl = plotKDE(fg, [:x0;:x1;:x2])
pl |> PNG("/tmp/test.png")

pl = plotKDE(fg, [:l0; :l1])

spyCliqMat(tree, :l0)

spyCliqMat(tree, :x2)

## specialized debugging

stuff = treeProductUp(fg, tree, :l0, :x0)
plotKDE(manikde!(stuff[1], (:Euclid,)) )

# plotTreeProductUp(fg, tree, :x1)

## Do one clique inference only

tree = buildTreeReset!(fg, drawpdf=true, show=true)
urt = doCliqInferenceUp!(fg, tree, :l0, false, iters=1, drawpdf=true)
upmsgs = urt.keepupmsgs

plotKDE([upmsgs[:x0]; upmsgs[:l1]; upmsgs[:x1]], c=["red";"green";"blue"])

## swap iteration order
#TODO guess order from bellow
# getCliqueData(tree.cliques[2]).itervarIDs = [5;7;3;1]
getCliqueData(tree.cliques[2]).itervarIDs = [:x0, :x1, :l0, :l1]

solveTree!(fg, tree)

## manually build the iteration scheme for second clique

# iter order:  x0, x1, l0, l1

stuff = treeProductUp(fg, tree, :l0, :x0)
X0 = manikde!(stuff[1], (:Euclid,))
plotKDE([X0; getBelief(fg, :x0)], c=["red";"green"])
setValKDE!(fg, :x0, X0)

stuff = treeProductUp(fg, tree, :l0, :x1)
X1 = manikde!(stuff[1], (:Euclid,))
plotKDE([X1; getBelief(fg, :x1)], c=["red";"green"])
setValKDE!(fg, :x1, X1)

stuff = treeProductUp(fg, tree, :l0, :l0)
L0 = manikde!(stuff[1], (:Euclid,))
plotKDE([L0; getBelief(fg, :l0)], c=["red";"green"])
setValKDE!(fg, :l0, L0)

stuff = treeProductUp(fg, tree, :l0, :l1)
L1 = manikde!(stuff[1], (:Euclid,))
plotKDE([L1; getBelief(fg, :l1)], c=["red";"green"])
setValKDE!(fg, :l1, L1)

## Reconstruct individual steps for broader clique factor selection

# Cliq 2:
# L0,X0,L1,X1
# x ,  ,  ,
# x ,x ,x ,
# x ,  ,x ,x
#   ,x ,  ,x
#   ,  ,x ,

# choose iteration order (priors last): :x0, :x1, :l0, :l1

# for new initialization format
# cliq 2: init :l0, :l1 directly from priors, them proceed with regular order
# cliq 1: initialize :x2 from incoming message singleton and proceed with regular order

# get factors for :x0 in clique2:
# :x0l0l1f1, :x0x1f1
ptsX0, = predictbelief(fg, :x0, [:x0l0l1f1; :x0x1f1])
X0 = manikde!(ptsX0, (:Euclid,))
plotKDE([X0; getBelief(fg, :x0)], c=["red";"green"])
setValKDE!(fg, :x0, X0)

# get factors for :x1 in clique2:
# :x1l0l1f1, :x0x1f1
ptsX1, = predictbelief(fg, :x1, [:x1l0l1f1, :x0x1f1])
X1 = manikde!(ptsX1, (:Euclid,))
plotKDE([X1; getBelief(fg, :x1)], c=["red";"green"])
setValKDE!(fg, :x1, X1)

# get factors for :l0
# :x0l0l1f1, :x1l0l1f1, :l0f1
ptsL0, = predictbelief(fg, :l0, [:x0l0l1f1, :x1l0l1f1, :l0f1])
L0 = manikde!(ptsL0, (:Euclid,))
plotKDE([L0; getBelief(fg, :l0)], c=["red";"green"])
setValKDE!(fg, :l0, L0)

# get factors for :l1
# :x0l0l1f1, :x1l0l1f1, :l1f1
ptsL1, = predictbelief(fg, :l1, [:x0l0l1f1, :x1l0l1f1, :l1f1])
L1 = manikde!(ptsL1, (:Euclid,))
plotKDE([L1; getBelief(fg, :l1)], c=["red";"green"])
setValKDE!(fg, :l1, L1)

## double check the upmessages are stored properly

##

plotLocalProduct(fg, :x0)
plotLocalProduct(fg, :x1)
plotLocalProduct(fg, :l0)
plotLocalProduct(fg, :l1)

##

initAll!(fg)

tree = buildTreeReset!(fg, drawpdf=true, show=true)
cliqorder = getCliqOrderUpSolve(tree)

spyCliqMat(cliqorder[end])

## Development zone

# treel = deepcopy(tree)
# fgl = deepcopy(fg)
# cliql = deepcopy(cliq)

##  develop better factor selection method

varlist = [:l0; :x0; :l1; :x1]
getFactorsAmongVariablesOnly(fg, varlist)

#


text:
using Meganet, BenchmarkTools

const history = Pkg.dir("Meganet")*"//benchmarks//micro//bm_resNN.jld"

TYPE = Float32
K = getConvGEMMKernel(TYPE,[96,96],[3,3,64,64])
nex = 20
# Bin = randn(TYPE,18,4)
# Bout = randn(TYPE,18,3)
nLayer = getBatchNormLayer(TYPE,[96*96,64],isTrainable=true)
lay     = getSingleLayer(TYPE,K,nLayer)

L = getResNN(TYPE,lay,4)

theta = initTheta(L)

Y     = randn(TYPE,nFeatIn(L),nex)

function benchmarkApply()
    funcName = "apply"
    Zd,Z,tmp  = apply(L,theta,Y,true);
    @code_warntype apply(L,theta,Y,true)

    trial = @benchmark apply(L,theta,Y,true)
    Meganet.updatehistory!(history, trial, funcName)
    hist = JLD.load(history, funcName)
    judge(hist)
end

###
function benchmarkJYTmv()
    funcName = "JYTmv"
    _,_,tmp = apply(L,theta,Y,true)
    Zout = randn(TYPE,nFeatOut(L),nex)

    #Warmup
    Z1 =  JYTmv(L,copy(Zout),(TYPE)[],theta,Y,tmp);
    trial = @benchmark JYTmv(L,copy(Zout),(TYPE)[],theta,Y,tmp)
    Meganet.updatehistory!(history, trial, funcName)
    hist = JLD.load(history, funcName)
    judge(hist)
end

function benchmarkJTmv()
    funcName = "JTmv"
    _,_,tmp = apply(L,theta,Y,true)
    Zout = randn(TYPE,nFeatOut(L),nex)

    #Warmup
    Z1 =  JTmv(L,copy(Zout),(TYPE)[],theta,Y,tmp);
    trial = @benchmark JTmv(L,copy(Zout),(TYPE)[],theta,Y,tmp)
    Meganet.updatehistory!(history, trial, funcName)
    hist = JLD.load(history, funcName)
    judge(hist)
end

function benchmarkJYmv()
    funcName = "JYmv"
    _,_,tmp = apply(L,theta,Y,true)
    Zout = randn(TYPE,nFeatOut(L),nex)

    #Warmup
    Z1 =  JYmv(L,copy(Zout),theta,Y,tmp);
    trial = @benchmark JYmv(L,copy(Zout),theta,Y,tmp)
    Meganet.updatehistory!(history, trial, funcName)
    hist = JLD.load(history, funcName)
    judge(hist)
end

function benchmarkJmv()
    funcName = "Jmv"
    _,_,tmp = apply(L,theta,Y,true)
    Zout = randn(TYPE,nFeatOut(L),nex)
    dtheta = randn(TYPE, size(theta))
    #Warmup
    Z1 =  Jmv(L,copy(dtheta),copy(Zout),theta,Y,tmp);
    trial = @benchmark Jmv(L,copy(dtheta),copy(Zout),theta,Y,tmp)
    Meganet.updatehistory!(history, trial, funcName)
    hist = JLD.load(history, funcName)
    judge(hist)
end

function benchmarkJthetamv()
    funcName = "Jthetamv"
    _,_,tmp = apply(L,theta,Y,true)
    dtheta = randn(TYPE, size(theta))

    #Warmup
    Z1 =  Jthetamv(L,dtheta,theta,Y,tmp);
    trial = @benchmark Jthetamv(L,dtheta,theta,Y,tmp)
    Meganet.updatehistory!(history, trial, funcName)
    hist = JLD.load(history, funcName)
    judge(hist)
end


text:
using CUDArt: CudaArray, HostArray
using CUBLAS
using CUBLAS: statuscheck, cublasStatus_t, cublasHandle_t, cublasOperation_t, BlasChar, cublasop, cublashandle
using CUBLAS: gemm!, gemv!, syr!, gemm_batched!, axpy!

const libcublas = CUBLAS.libcublas

for (fname, elty) in
        ((:cublasDgemmStridedBatched, :Float64),
         (:cublasSgemmStridedBatched, :Float32),
         (:cublasZgemmStridedBatched, :Complex128),
         (:cublasCgemmStridedBatched, :Complex64))
    @eval begin
        #=
        cublasStatus_t  cublasDgemmStridedBatched(
        cublasHandle_t handle,
        cublasOperation_t transa, cublasOperation_t transb,
        int m, int n, int k,
        const double *alpha,
        const double *A, int lda, long long int strideA,
        const double *B, int ldb, long long int strideB,
        const double *beta,
        double *C, int ldc, long long int strideC, int batchCount);
        =#
        function gemm_strided_batched!(transA::BlasChar,
                                       transB::BlasChar,
                                       alpha::($elty),
                                       A::CudaArray{$elty, 3},
                                       B::CudaArray{$elty, 3},
                                       beta::($elty),
                                       C::CudaArray{$elty, 3})

            m = size(A, transA == 'N' ? 1 : 2)
            k = size(A, transA == 'N' ? 2 : 1)
            n = size(B, transB == 'N' ? 2 : 1)
            if m != size(C, 1) || n != size(C, 2) || k != size(B, transB == 'N' ? 1 : 2)
                throw(DimensionMismatch(""))
            end

            strideA = size(A, 1) * size(A, 2)
            strideB = size(B, 1) * size(B, 2)
            strideC = size(C, 1) * size(C, 2)
            batchCount = size(A, 3)

            cutransA = cublasop(transA)
            cutransB = cublasop(transB)
            lda = max(1,stride(A,2))
            ldb = max(1,stride(B,2))
            ldc = max(1,stride(C,2))

            statuscheck(ccall(($(string(fname)), libcublas),
                cublasStatus_t,
                (cublasHandle_t, cublasOperation_t, cublasOperation_t,
                Cint, Cint, Cint, # m, n, k
                Ptr{$elty},       # alpha
                Ptr{$elty}, Cint, Clonglong, # A, lda, strideA
                Ptr{$elty}, Cint, Clonglong, # B, ldb, strideB
                Ptr{$elty},       # beta
                Ptr{$elty}, Cint, Clonglong, # C, ldc, strideC
                Cint),  # batchCount
                cublashandle[1], cutransA, cutransB,
                m, n, k,
                [alpha],
                A, lda, strideA,
                B, ldb, strideB,
                [beta],
                C, ldc, strideC,
                batchCount))
            C
        end
    end
end

"""
```
function takelast!{F <: AbstractFloat, T <: Integer}(
    dest::CudaArray{F}, src::CudaArray{F},
    idxes::Union{Tuple{Vararg{T}}, Array{T}})
```
get subarray by last dimension.
"""
function takelast!{F <: AbstractFloat, T <: Integer}(
    dest::CudaArray{F}, src::CudaArray{F},
    idxes::Union{Tuple{Vararg{T}}, Array{T}})

    siz = size(src)[1:end-1]
    strides_ = prod(siz) * sizeof(eltype(src))

    for (i, k) in enumerate(idxes)
        dest_ptr = dest.ptr + (i-1) * strides_
        src_ptr = src.ptr + (k-1) * strides_
        src_ = CudaArray(src_ptr, siz, src.dev)
        dest_ = CudaArray(dest_ptr, siz, dest.dev)
        copy!(dest_, src_)
    end
    dest
end

function viewlast{F <: AbstractFloat}(src::CudaArray{F}, idx::Integer)
    siz = size(src)[1:end-1]
    strides_ = prod(siz) * sizeof(eltype(src))

    src_ptr = src.ptr + (idx-1) * strides_
    dest = CudaArray(src_ptr, siz, src.dev)
    dest
end

function splitarray(a::CudaArray, shape1::Tuple, shape2::Tuple) :: NTuple{2, CudaArray}
    l1 = prod(shape1)
    l2 = prod(shape2)
    @assert length(a) == l1 + l2
    p1_ptr = a.ptr
    p2_ptr = p1_ptr + l1 * sizeof(eltype(a))
    p1 = CudaArray(p1_ptr, shape1, a.dev)
    p2 = CudaArray(p2_ptr, shape2, a.dev)
    return (p1, p2)
end


text:
#= Fields for OceanTurb.jl

A `Field` is an abstraction of a variable or function defined on a staggered grid.
Our hope is that it simplifies differential operations on a staggered grid.

OceanTurb.jl solves one-dimensional PDEs on a staggered grid.
The geometry of a grid with `N=3` is

```
      ▲ z
      |

         i=4           *
                j=4   ===       ▲
         i=3           *        | Δf (i=3)
                j=3   ---       ▼
         i=2           *    ▲
                j=2   ---   | Δc (j=2)
         i=1           *    ▼
                j=1   ===
         i=0           *
```

where the i's index cells and the j's index faces.
The variable Δc gives the separation between
cell centers, and Δf gives the separation between faces.
Ghost cells at i=0 and i=N+1 bound the domain.

There are two types of fields:

  1. Fields defined at cell centers with dimension `N+2`: `Field{Cell}`
  2. Fields defined at cell interfaces with dimension `N+1`: `Field{Face}`
=#
import Base: +, *, -, ^, setindex!, getindex, eachindex, lastindex, similar,
             eltype, length, @propagate_inbounds

import Statistics: mean

default_arraytype(T) = Array{T, 1}

struct CellField{A, G, T} <: AbstractField{A, G, T}
    data :: OffsetArray{T, 1, A}
    grid :: G
end

struct FaceField{A, G, T} <: AbstractField{A, G, T}
    data :: A
    grid :: G
end

arraytype(::AbstractField{A}) where A = A
eltype(::AbstractField{A}) where A = eltype(A)

#
# (legacy) Field Location and Constructors
#

abstract type FieldLocation end
struct Cell <: FieldLocation end
struct Face <: FieldLocation end

Field(::Type{Face}, grid) = FaceField(grid)
Field(::Type{Cell}, grid) = CellField(grid)

"""
    FaceField(grid)

Return a `FaceField` on `grid` with its data initialized to 0.
"""
function FaceField(A::DataType, grid)
    data = convert(A, fill(0, face_size(grid)))
    FaceField{typeof(data), typeof(grid), eltype(data)}(data, grid)
end

"""
    CellField(grid)

Return a `CellField` on `grid` with its data initialized to 0.
"""
function CellField(A::DataType, grid)
    data = convert(A, fill(0, cell_size(grid)))
    offset_data = OffsetArray(data, 0:grid.N+1)
    CellField{typeof(data), typeof(grid), eltype(data)}(offset_data, grid)
end

CellField(grid) = CellField(default_arraytype(eltype(grid)), grid)
FaceField(grid) = FaceField(default_arraytype(eltype(grid)), grid)

"""
    CellField(data, grid)

Return a `CellField` with its `data` located on the `grid`.
if `data` is an array, it must be broadcastable to `c.data`, where
`c` is a `CellField`.
"""
function CellField(data, grid)
    c = CellField(grid)
    set!(c, data)
    return c
end

"""
    FaceField(data, grid)

Return a `FaceField` with its `data` located on the `grid`.
if `data` is an array, it must be broadcastable to `f.data`, where
`f` is a `FaceField`.
"""
function FaceField(data, grid)
    f = FaceField(grid)
    set!(f, data)
    return f
end

#
# Basic 'Field' functionality
#

data(c::FaceField) = c.data
data(c::CellField) = view(c.data, 1:c.grid.N)

nodes(c::CellField) = c.grid.zc
nodes(f::FaceField) = f.grid.zf

length(c::CellField) = c.grid.N
length(f::FaceField) = f.grid.N + 1

height(c::AbstractField) = height(c.grid)

# All indices
eachindex(c::CellField) = 1:c.grid.N
eachindex(f::FaceField) = 1:f.grid.N + 1

lastindex(c::CellField) = c.grid.N
lastindex(f::FaceField) = f.grid.N + 1

# Interior indices, omitting boundary-adjacent values
interiorindices(c::CellField) = 2:c.grid.N - 1
interiorindices(f::FaceField) = 2:f.grid.N
boundaryindices(c::CellField) = (1, c.grid.N)

# Sugary sweet: access indices of c.data by indexing into c.
@propagate_inbounds getindex(c::AbstractField, inds...) = getindex(c.data, inds...)
@propagate_inbounds setindex!(c::AbstractField, d, inds...) = setindex!(c.data, d, inds...)
@propagate_inbounds setindex!(c::AbstractField, d::AbstractField, inds...) = setindex!(c.data, d.data, inds...)

#
# Ways to specify a field's data
#

set!(c::AbstractField, data::Number) = fill!(c.data, data)
set!(c::AbstractField{Ac, G}, d::AbstractField{Ad, G}) where {Ac, Ad, G} = c.data .= convert(Ac, d.data)
set!(c::FaceField, fcn::Function) = c.data .= fcn.(nodes(c))

function set!(c::CellField, func::Function)
    data = func.(nodes(c))
    set!(c, data)
    # Set ghost points to get approximation to first derivative at boundary
    data_bottom = func(c.grid.zf[1])
    data_top = func(c.grid.zf[end])

    # Set ghost values so that
    # ∂z(c, 1) = (c[1] - c[0]) / Δc(c, 1) = (c[1] - c_bottom) / 0.5*Δc(c, 1)
    #
    # and
    # ∂z(c, N+1) = (c[N+1] - c[N]) / Δc(c, N+1) = (c_top - c[N]) / 0.5*Δc(c, N)

    N = c.grid.N
    @inbounds begin
        c[0] = c[1] - 2 * (c[1] - data_bottom)
        c[N+1] = c[N] + 2 * (data_top - c[N])
    end

    return nothing
end

set!(f::FaceField, data::AbstractArray) = f.data .= data

function set!(c::CellField, data::AbstractArray)
    for i in eachindex(data)
        @inbounds c[i] = data[i]
    end
    # Default boundary conditions if data is not an OffsetArray
    typeof(data) <: OffsetArray || set_default_bcs!(c)
    return nothing
end

function set_default_bcs!(c)
    @inbounds begin
        c[0] = c[1]
        c[c.grid.N+1] = c[c.grid.N]
    end
    return nothing
end

function integral(fn::Function, c::CellField)
    total = zero(eltype(c))
    for i in eachindex(c)
        @inbounds total += fn(c[i]) * Δf(c, i)
    end
    return total
end

integral(c::CellField) = integral(x->x, c)

"""
    mean([f], c::CellField)

Compute the mean of the field `c` over its domain,
applying the function `f` to each element.
`f` is the identity function f(x) = x by default.
"""
mean(fn::Function, c::CellField) = integral(fn, c) / height(c)
mean(c::CellField) = mean(x->x, c)

function integrate_range(c::CellField, i₁::Int, i₂::Int)
    total = 0
    for i = i₁:i₂
        @inbounds total += c[i] * Δf(c.grid, i)
    end
    return total
end

function integral(c::CellField, z₋, z₊=0)

    @assert z₊ > c.grid.zf[1] "Integration region lies outside the domain."
    @assert z₊ > z₋ "Invalid integration range: upper limit greater than lower limit."

    # Find region bounded by the face ≤ z₊ and the face ≤ z₁
    i₁ = searchsortedfirst(c.grid.zf, z₋) - 1
    i₂ = searchsortedfirst(c.grid.zf, z₊) - 1

    if i₂ ≠ i₁
        # Calculate interior integral, recalling that the
        # top interior cell has index i₂-2.
        total = integrate_range(c, i₁+1, i₂-1)

        # Add contribution to integral from fractional bottom part,
        # if that region is a part of the grid.
        if i₁ > 0
            total += c[i₁] * (c.grid.zf[i₁+1] - z₋)
        end

        # Add contribution to integral from fractional top part
        total += c[i₂] * (z₊ - c.grid.zf[i₂])
    else
        total = c[i₁] * (z₊ - z₋)
    end

    return total
end

function set!(c1::CellField{A1, G2}, c2::CellField{A2, G2}) where {A1, G1, A2, G2}
    @assert height(c1) == height(c2) "Physical domains differ between the two fields."
    for i in eachindex(c1)
        @inbounds c1[i] = integral(c2, c1.grid.zf[i], c1.grid.zf[i+1]) / Δf(c1, i)
    end
    return nothing
end

similar(c::CellField) = CellField(c.grid)
similar(f::FaceField) = FaceField(f.grid)

# Define +, -, and * on fields as element-wise calculations on their data. This
# is only true for fields of the same type. So far, we haven't found use for
# these sweets because we tend to write element-wise kernels for operations.
for op in (:+, :-, :*)
    @eval begin
        # +, -, * a Field by a Number on the left
        function $op(num::Number, f::AbstractField)
            ff = similar(f)
            @. ff.data = $op(num, f.data)
            ff
        end

        # +, -, * a Field by a Number on the right.
        $op(f::AbstractField, num::Number) = $op(num, f)

        # Binary two-field operations
        function $op(f1::F, f2::F) where {F <: AbstractField}
            f3 = similar(f1)
            @. f3.data = $op(f1.data, f2.data)
            f3
        end
    end
end

function ^(c::AbstractField, b::Number)
    d = similar(c)
    set!(d, c.data.^b)
    return d
end

#
# Differential operators and such for fields
#

@propagate_inbounds Δc(c::AbstractField, i_face) = Δc(c.grid, i_face)
@propagate_inbounds Δf(c::AbstractField, i_cell) = Δf(c.grid, i_cell)

"Upwards-biased difference."
∂z⁺(i, grid, f::Function, args...) = (f(i+1, grid, args...) - f(i, grid, args...)) / Δc(grid, i+1)

"""
    ∂z(a, i)

Return the discrete derivative of `a` at grid point `i`.

The derivative of a `CellField` is computed at face points,
and the derviative of a `FaceField` is computed at cell points.
"""
∂z(a, i) = throw("∂z is not defined for arbitrary fields.")

"Return ∂c/∂z at face index i."
@propagate_inbounds ∂z(c::CellField, i) = (c.data[i] - c.data[i-1]) / Δc(c, i)

"Return ∂c/∂z at face index i."
@propagate_inbounds ∂z(c::FaceField, i) = (c.data[i+1] - c.data[i]) / Δc(c, i)
@propagate_inbounds ∂²z(c::AbstractField, i) = (∂z(c, i+1) - ∂z(c, i)) / Δf(c, i)

"Calculate `f = ∂c/∂z` in the grid interior."
function ∂z!(f::FaceField, c::CellField)
    for i = eachindex(f)
        @inbounds f.data[i] = ∂z(c, i)
    end
    return nothing
end

"Calculate `c = ∂f/∂z` in the grid interior."
function ∂z!(c::CellField, f::FaceField)
    for i = eachindex(c)
        @inbounds c.data[i] = ∂z(f, i)
    end
    return nothing
end

"Return the `FaceField` ∂c/∂z, where `c` is a `CellField`."
function ∂z(c::CellField)
    f = FaceField(c.grid)
    ∂z!(f, c)
    return f
end

"Return the `CellField` ∂f/∂z, where `f` is a `FaceField`."
∂z(a) = throw("∂z is not defined for arbitrary fields.")

function ∂z(f::FaceField)
    c = CellField(f.grid)
    ∂z!(c, f)
    return c
end

#
# Advection and diffusion operators
#

@propagate_inbounds K∂z(K, ϕ, i) = K * ∂z(ϕ, i)

"Return the diffusive flux divergence at cell i."
@propagate_inbounds ∇K∇ϕ(Kᵢ₊₁, Kᵢ, ϕ, i) = (K∂z(Kᵢ₊₁, ϕ, i+1) - K∂z(Kᵢ, ϕ, i)) / Δf(ϕ, i)

"Return the upwind advective flux divergence at cell i for M<0."
@propagate_inbounds ∂zM(Mᵢ₊₁, Mᵢ, ϕ, i) = (Mᵢ₊₁ * ϕ[i+1] - Mᵢ * ϕ[i]) / Δc(ϕ, i+1)

"Return the total flux (advective + diffusive) across face i."
@propagate_inbounds flux(M, K, ϕ, i) = M * onface(ϕ, i) - K * ∂z(ϕ, i)

#
# Convenience functions
#

top(a) = a
top(a::AbstractArray) = @inbounds a[end]
top(a::CellField) = @inbounds a[a.grid.N]
top(a::FaceField) = @inbounds a[a.grid.N+1]

bottom(a) = throw("bottom(a) Not implemented for typeof(a) = $(typeof(a)).")
bottom(a::Number) = a
bottom(a::Union{AbstractField, AbstractArray}) = @inbounds a[1]

"""
    onface(c, i)

Return the interpolation of `c` onto face point `i`.
"""
@propagate_inbounds onface(c::CellField, i) = (c.data[i] + c.data[i-1])/2
@propagate_inbounds onface(f::FaceField, i) = f[i]

"""
    oncell(f, i)

Return the interpolation of `f` onto cell point `i`.
"""
@propagate_inbounds oncell(f::FaceField, i) = (f.data[i+1] + f.data[i])/2
@propagate_inbounds oncell(c::CellField, i) = c[i]

@inline oncell(i, grid, f::Function, args...) = ( f(i+1, grid, args...) + f(i, grid, args...) ) / 2
@inline onface(i, grid, f::Function, args...) = ( f(i, grid, args...) + f(i-1, grid, args...) ) / 2

"""
    absolute_error(c, d, p=2)

Compute the absolute error between `c` and `d` with norm `p`, defined as

``\\mathrm{abs \\, error} = \\left ( L^{-1} \\int_{-L}^0 (c-d)^p \\, \\mathrm{d} z \\right )^(1/p)``.
"""
function absolute_error(c::CellField, d::CellField, p=2)
    if length(c) != length(d)
        𝒹 = similar(c)
        set!(𝒹, d)
    else
        𝒹 = d
    end

    total = zero(eltype(c))
    for i in eachindex(c)
        @inbounds total += (c[i] - 𝒹[i])^p * Δf(c, i)
    end

    return  ( total / height(c) )^(1/p)
end

"""
    relative_error(c, d, p=2)

Compute the relative error between `c` and `d` with norm `p`, defined as

```math
\\beq
\\mathrm{rel \\, error} = \\frac{\\left ( int_{-L}^0 (c-d)^p \\, \\mathrm{d} z \\right )^(1/p)}
                             {\\left ( int_{-L}^0 d^p \\, \\mathrm{d} z \\right )^(1/p)}
\\eeq
```
"""
relative_error(c::CellField, d::CellField, p=2) = absolute_error(c, d, p) / mean(x -> x^p, d)^(1/p)


text: hist(randn(100); ecol="red", fill="wheat")


text:
# HSA dispatch packet offsets
_packet_names = fieldnames(HSA.KernelDispatchPacket)
_packet_offsets = fieldoffset.(HSA.KernelDispatchPacket, 1:length(_packet_names))

@generated function _intr(::Val{fname}, out_arg, inp_args...) where {fname,}
    Context() do ctx
        inp_exprs = [:( inp_args[$i] ) for i in 1:length(inp_args)]
        inp_types = [inp_args...]
        out_type = convert(LLVMType, out_arg.parameters[1]; ctx)

        # create function
        bool_types = map(x->x===Bool, inp_types)
        T_bool = LLVM.Int1Type(ctx)
        param_types = LLVMType[convert.(LLVMType, inp_types; ctx=ctx)...]
        llvm_f, _ = create_function(out_type, param_types)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            # call the intrinsic
            intr_types = map(x->bool_types[x[1]] ? T_bool : x[2], enumerate(param_types))
            intr_out_type = out_arg === Type{Bool} ? T_bool : out_type
            intr_ftype = LLVM.FunctionType(intr_out_type, intr_types)
            intr = LLVM.Function(mod, string(fname), intr_ftype)
            for idx in 1:length(param_types)
                if bool_types[idx]
                    attrs = parameter_attributes(intr, idx)
                    push!(attrs, EnumAttribute("zeroext", 0; ctx))
                end
            end
            params = map(x->bool_types[x[1]] ? trunc!(builder, x[2], T_bool) : x[2], enumerate(parameters(llvm_f)))
            value = call!(builder, intr, [params...])
            if out_arg === Type{Bool}
                ret!(builder, zext!(builder, value, convert(LLVMType, Bool; ctx)))
            else
                ret!(builder, value)
            end
        end

        call_function(llvm_f, out_arg.parameters[1], Tuple{inp_args...}, inp_exprs...)
    end
end

struct GCNIntrinsic
    jlname::Symbol
    rocname::Symbol
    isbroken::Bool # please don't laugh...
    isinverted::Bool
    # FIXME: Input/output types
    inp_args::Tuple
    out_arg::Type
    roclib::Symbol
    suffix::Symbol
    tobool::Bool
end

GCNIntrinsic(jlname, rocname=jlname; isbroken=false, isinverted=false,
             inp_args=(), out_arg=(), roclib=:ocml, suffix=fntypes[first(inp_args)], tobool=false) =
    GCNIntrinsic(jlname, rocname, isbroken, isinverted, inp_args, out_arg, roclib, suffix, tobool)

import ExprTools

macro device_override(ex)
    esc(quote
        $GPUCompiler.@override($method_table, $ex)
    end)
end

macro device_function(ex)
    ex = macroexpand(__module__, ex)
    def = ExprTools.splitdef(ex)

    # generate a function that errors
    def[:body] = quote
        error("This function is not intended for use on the CPU")
    end

    esc(quote
        $(ExprTools.combinedef(def))
        @device_override $ex
    end)
end


text:
#################
# LinearIndices #
#################
toindices(x::Int) = LinearIndices((x,))
ncomponents(x::LinearIndices) = length(x)

####################
# SymmetricIndices #
####################
struct SymmetricIndices{order, dim} <: AbstractArray{Int, order}
    sym::Symmetry{NTuple{order, dim}}
end
SymmetricIndices(dims::Int...) = SymmetricIndices(Symmetry(dims))
SymmetricIndices(dims::Tuple{Vararg{Int}}) = SymmetricIndices(Symmetry(dims))

toindices(sym::Symmetry) = SymmetricIndices(sym)
ncomponents(x::SymmetricIndices) = ncomponents(x.sym)

Base.size(x::SymmetricIndices) = Dims(x.sym)
function Base.getindex(s::SymmetricIndices{order}, I::Vararg{Int, order}) where {order}
    @boundscheck checkbounds(s, I...)
    sorted = sort!(collect(I), rev = true) # `reverse` is for column-major order
    @inbounds LinearIndices(size(s))[CartesianIndex{order}(sorted...)]
end

struct TensorIndices{order, S <: Tuple{Vararg{Union{Int, Symmetry}}}} <: AbstractArray{Int, order}
    dims::NTuple{order, Int}
    size::S
end

TensorIndices(S::Tuple{Vararg{Union{Int, Symmetry}}}) = TensorIndices(flatten_tuple(map(Dims, S)), S)
TensorIndices(S::Vararg{Union{Int, Symmetry}}) = TensorIndices(S)

Base.size(x::TensorIndices) = x.dims

function Base.getindex(t::TensorIndices{order}, I::Vararg{Int, order}) where {order}
    @boundscheck checkbounds(t, I...)
    S = t.size
    st = 1
    inds = Vector{Int}(undef, length(S))
    @inbounds begin
        for (i, s) in enumerate(S)
            x = toindices(s)
            n = ndims(x)
            inds[i] = x[I[st:(st+=n)-1]...]
        end
        LinearIndices(prod.(Dims.(S)))[inds...]
    end
end

function _independent_indices(inds::TensorIndices)
    dict = Dict{Int, Int}()
    map(inds) do i
        get!(dict, i, length(dict) + 1)
    end
end

function _indices(inds::TensorIndices)
    unique(inds)
end

function _duplicates(inds::TensorIndices)
    dups = Dict{Int, Int}()
    for i in inds
        if !haskey(dups, i)
            dups[i] = 1
        else
            dups[i] += 1
        end
    end
    map(x -> x.second, sort(collect(dups), by = x->x[1]))
end

for func in (:independent_indices, :indices, :duplicates)
    @eval @generated function $func(::Space{S}) where {S}
        arr = $(Symbol(:_, func))(TensorIndices(S))
        quote
            SArray{Tuple{$(size(arr)...)}, Int}($arr)
        end
    end
end


text:
using BisectPy:
    bisect_left,
    bisect_right,
    index,
    find_lt,
    find_le,
    find_gt,
    find_ge,
    insort_left,
    insort_right
using Test

@testset "BisectPy.jl" begin
    @test bisect_left([1, 2, 3, 4, 5], 3.5) == 4

    @test bisect_right([1, 2, 3, 4, 5], 3.5) == 4

    @testset "Equal tests" begin
        @test bisect_left([1, 2, 3, 4, 5], 2) == 2
        @test bisect_left([1, 2, 3, 3, 3, 5], 3) == 3
        @test bisect_right([1, 2, 3, 4, 5], 2) == 3
        @test bisect_right([1, 2, 3, 3, 3, 5], 3) == 6
    end

    @testset "Boundary tests" begin
        @test bisect_left([1, 2, 3, 4, 5], 0) == 1
        @test bisect_left([1, 2, 3, 4, 5], 5) == 5
        @test bisect_left([1, 2, 3, 4, 5], 6) == 6
        @test bisect_right([1, 2, 3.0, 4, 5], 5) == 6
        @test bisect_right([1, 2, 3.0, 4, 5], 6) == 6
        @test bisect_right([1, 2, 3, 4, 5], 1) == 2
        @test bisect_right([1, 2, 3, 4, 5], 0) == 1
    end

    @testset "`Python example 1`" begin
        function grade(score, breakpoints = [60, 70, 80, 90], grades = "FDCBA")
            i = bisect_right(breakpoints, score)
            return grades[i]
        end
        @test [grade(score) for score in [33, 99, 77, 70, 89, 90, 100]] == ['F', 'A', 'C', 'C', 'B', 'A', 'A']
    end

    @testset "`Python example 2`" begin
        data = [("red", 5), ("blue", 1), ("yellow", 8), ("black", 0)]
        sort!(data)
        keys = [r[2] for r in sort(data)]
        @test data[bisect_left(keys, 0)] == ("black", 0)
        @test data[bisect_left(keys, 1)] == ("blue", 1)
        @test data[bisect_left(keys, 5)] == ("red", 5)
        @test data[bisect_left(keys, 8)] == ("yellow", 8)
    end

    @test insort_left([1, 2, 3, 4, 5], 0) == [0, 1, 2, 3, 4, 5]
    @test insort_left([1, 2, 3, 4, 5], 6) == [1, 2, 3, 4, 5, 6]
    @test insort_right([1, 2, 3, 4, 5], 5) == [1, 2, 3, 4, 5, 5]
    @test insort_right([1, 2, 3, 4, 5], 6) == [1, 2, 3, 4, 5, 6]

    @test index([1, 2, 3, 4, 5], 3.5) === nothing
    @test index([1, 2, 2, 3, 4, 5], 2) == 2
    @test index([1, 2, 3, 4, 5], 5) == 5
    @test index([1, 2, 3, 4, 5], 6) === nothing

    @test find_lt([1, 2, 3, 4, 5], 2) == 1
    @test find_lt([1, 2, 2, 3, 4, 4, 5], 3) == 2
    @test find_lt([1, 2, 2, 3, 3.5, 4, 4, 5], 3.5) == 3
    @test find_lt([1, 2, 3, 4, 5], 0) === nothing
    @test find_lt([1, 2, 3, 4, 5], 6) == 5

    @test find_le([1, 2, 3, 4, 5], 2) == 2
    @test find_le([1, 2, 2, 3, 4, 4, 5], 3) == 3
    @test find_le([1, 2, 2, 3, 4, 5, 5], 3.5) == 3
    @test find_le([1, 2, 3, 4, 5], 0) === nothing
    @test find_le([1, 2, 3, 4, 5], 6) == 5

    @test find_gt([1, 2, 3, 4, 5], 2) == 3
    @test find_gt([1, 2, 2, 3, 4, 4, 5], 3) == 4
    @test find_gt([1, 2, 2, 3, 4, 5, 5], 3.5) == 4
    @test find_gt([1, 2, 3, 4, 5], 0) == 1
    @test find_gt([1, 2, 3, 4, 5, 5], 6) === nothing

    @test find_ge([1, 2, 3, 4, 5], 2) == 2
    @test find_ge([1, 2, 2, 3, 4, 4, 5], 3) == 3
    @test find_ge([1, 2, 2, 3, 4, 5, 5], 3.5) == 4
    @test find_ge([1, 2, 3, 4, 5], 0) == 1
    @test find_ge([1, 2, 3, 4, 5, 5], 6) === nothing
end


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#john-addition-test.jl

include("../unum.jl")
using Unums

#do addition of random numbers in a set of environments

const iters = 100000

#set up arrays of values so that generating arrays is not part of the speed calc

function exprand()
  exp(randn() * 100) * (rand() > 0.5 ? 1: - 1)
end

function add_test(T::Type)

  a = [convert(T, exprand()) for i=1:iters]
  b = [convert(T, exprand()) for i=1:iters]

  println("testing $T")
  tic()
    #i don't think that llvm aggressively optimizes this out since we don't use the result.
    for (idx = 1:iters)
      a[idx] + b[idx]
    end
  toc()
end

function random_with_ulp(T::Type)
  x = convert(T, rand())
  (rand() > 0.5) ? x : unum_unsafe(x, x.flags | Unums.UNUM_UBIT_MASK)
end

function add_test_ulp(T::Type)
  a = [random_with_ulp(T) for i=1:iters]
  b = [random_with_ulp(T) for i=1:iters]

  println("testing with ulps $T")
  tic()
    #i don't think that llvm aggressively optimizes this out since we don't use the result.
    for (idx = 1:iters)
      a[idx] + b[idx]
    end
  toc()
end

Unums.__unum_release_environment()
#map the add_test function onto an array of types we want to run it on!
map(add_test, [Float64, BigFloat, Unum{0,0}, Unum{1,1}, Unum{3,3}, Unum{4,6}, Unum{4,7}, Unum{4,8}])
map(add_test_ulp, [Unum{0,0}, Unum{1,1}, Unum{3,3}, Unum{4,6}, Unum{4,7}, Unum{4,8}])

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using StateSpaceModels, Plots

y = collect(1:0.25:20) + 0.5*randn(77)
y[10:20] .= NaN
model = linear_trend(y)
ss = statespace(model)

plot([y y], legend = :topleft, label = "observations", lw = 2, layout = (2, 1), xticks = 0:10:77, color = "black", grid = false)
plot!([ss.filter.a[:, 1] ss.smoother.alpha[:, 1]], label = ["filtered state" "smoothed state"],
        color = ["indianred" "green"], lw = 2)


text:
export func_transform, euler_matrix,euler_matrix_update, nonlinear_transform, make_input_vector, make_hermitian
"""
    make_input_vector(x::Vector{T}) -> Vector{T}

Generates input vector for `QuNLDE` iterations.
"""
function make_input_vector(x::Vector{T}) where T
    if !(norm(x) ≈ 1)
        throw(ArgumentError("Input vector is not normalized"))
    end
    len = length(x)
    siz = nextpow(2,len + 1)
    z = zeros(T, siz)
    z[1] = one(T)
    z[2:len+1] = x
    normalize!(z)
    reg = [1,0] ⊗ z ⊗ z
    return reg
end

"""
    make_hermitian(A::Matrix) -> Matrix

Returns hermitian matrix containing A.
"""
function make_hermitian(A::Matrix)
    if ishermitian(A)
        return A
    end

    return [zero(A) im*A'; -im*A zero(A)]
end

function nonlinear_transform(H::Matrix, x::Vector, k::Int, ϵ::Real = 1e-4)
    r, N = taylorsolve(im*H,x,k,ϵ)
    n = log2i(length(x))
    nb = Int((n-1)/2)
    r = relax!(r) |> focus!(n)|> select!(1) |> focus!(1:nb...,) |> select!(0)
    return r,sqrt(2)*N/ϵ
end

"""
    func_transform(A::Matrix, x::Vector, k::Int) -> ArrayReg, <: Complex

Function transform sub-routine. Returns state register and inverse probability of finding  it.
"""
function func_transform(A::Matrix, x::Vector, k::Int,ϵ::Real = 1e-4)
    reg = make_input_vector(x)
    H = make_hermitian(A)
    r, N = nonlinear_transform(H,reg,k,ϵ)
    return r, N
end

"""
    euler_matrix(A::Matrix{CPType},b::Vector{CPType},h::Real) -> Matrix

Generates matrix for forward Euler iteration.
"""
function euler_matrix(A::Matrix{CPType},b::Vector{CPType},h::Real) where CPType
    n = length(b)
    A = CPType(h)*A
    A[1,1] = 1
    @inbounds for i in 1:n
        A[4*i+ 1, i+1] += 1
    end
    return A
end

"""
    euler_matrix(A::Matrix{CPType},b::Vector{CPType},h::Real) -> Matrix

Updates euler matrix for forward Euler iteration.
"""
function euler_matrix_update(A::Matrix{CPType}, b::Vector{CPType}, nrm::Real) where CPType
    n = length(b)
    A = CPType(nrm^2)*A
    A[1,1] = 1
    @inbounds for i in 1:n
        A[4*i+ 1, 1] = A[4*i+ 1, 1]/(nrm^2)
        @. A[4*i+ 1, 2:n+1] = A[4*i+ 1, 2:n+1]/nrm
        for j in 1:n
            A[4*i + 1, 4*j + 1] = A[4*i + 1, 4*j + 1]/nrm
        end
    end
    return A
end

function DiffEqBase.solve(prob::QuODEProblem,alg::QuNLDE; dt = (prob.tspan[2]-prob.tspan[1])/10, kwargs...)
    A = prob.A
    b = prob.b
    k = alg.k
    ϵ = alg.ϵ
    tspan = prob.tspan
    len = round(Int,(tspan[2] - tspan[1])/dt) + 1
    siz = length(b)
    res = Array{eltype(b),2}(undef,len,siz)
    res[1,:] = b
    A = euler_matrix(A,b,dt)
    H = make_hermitian(A)
    reg = make_input_vector(b)
    r, N = nonlinear_transform(H,reg,k,ϵ)
    ntem = 1
    @views for step in 2:len - 1
        tem = vec(state(r))*N*sqrt(2)
        res[step,:] = tem[2:siz+1]
        ntem = norm(res[step,:])
        C = euler_matrix_update(A,b,ntem)
        H = make_hermitian(C)
        reg = res[step,:]/ntem
        reg = make_input_vector(reg)
        r, N = nonlinear_transform(H,reg,k,ϵ)
    end
    tem = vec(state(r))*N*sqrt(2)
    @views res[len,:] = tem[2:siz+1]
    return res
end


text:
function test_oilmm(rng, kernels, H::Orthogonal, x_train, x_test, y_train, y_test)
    fs = independent_mogp(map(GP, kernels))

    ilmm = ILMM(fs, collect(H))
    oilmm = ILMM(fs, H)

    ilmmx = ilmm(x_train, 0.1)
    oilmmx = oilmm(x_train, 0.1)

    @test isapprox(mean(ilmmx), mean(oilmmx))
    @test isapprox(var(ilmmx), var(oilmmx))
    @test isapprox(cov(ilmmx), cov(oilmmx))
    @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train))
    @test _is_approx(marginals(ilmmx), marginals(oilmmx))
    @test length(rand(rng, oilmmx)) == size(H, 1) * length(x_train.x)

    p_ilmmx = posterior(ilmmx, y_train)
    p_oilmmx = posterior(oilmmx, y_train)

    pi = p_ilmmx(x_test, 0.1)
    po = p_oilmmx(x_test, 0.1)

    @test isapprox(mean(pi), mean(po))
    @test isapprox(var(pi), var(po))
    @test isapprox(logpdf(pi, y_test), logpdf(po, y_test))
    @test _is_approx(marginals(pi), marginals(po))
    @test length(rand(rng, po)) == size(H, 1) * length(x_test.x)

    test_sampling_consistency(rng, oilmm, x_train)

    @test gradient(logpdf, oilmmx, y_train) isa Tuple
    @test gradient(logpdf, po, y_test) isa Tuple

    @testset "primary_public_interface" begin
        test_finitegp_primary_and_secondary_public_interface(rng, oilmmx)
        test_finitegp_primary_and_secondary_public_interface(rng, po)
    end
end

@testset "oilmm" begin
    rng = Random.seed!(04161999)
    x_train, x_test, y_train, y_test = generate_toy_data(rng)

    @testset "Full Rank, Dense H" begin
        U, S, _ = svd(rand(rng, 3, 3))
        H = Orthogonal(U, Diagonal(S))
        kernels = [SEKernel(), Matern32Kernel(), Matern32Kernel()]
        test_oilmm(rng, kernels, H, x_train, x_test, y_train, y_test)
    end

    @testset "M Latent Processes" begin
        U, S, _ = svd(rand(rng, 3, 2))
        H = Orthogonal(U, Diagonal(S))
        kernels = [SEKernel(), Matern32Kernel()]
        test_oilmm(rng, kernels, H, x_train, x_test, y_train, y_test)
    end

    @testset "1 Latent Processes" begin
        U, S, _ = svd(rand(rng, 3, 1))
        H = Orthogonal(U, Diagonal(S))
        kernels = [SEKernel()]
        test_oilmm(rng, kernels, H, x_train, x_test, y_train, y_test)
    end
end
@info "Ran oilmm tests."