text:
function _random_invertible_matrix(n, B)
  A = identity_matrix(FlintZZ, n)
  if n == 1
    return A
  end
  for k in 1:10
    i = rand(1:(n - 1))
    j = rand((i + 1):n)
    A[i, j] += rand(B)
  end
  @assert det(A) == 1
  C = identity_matrix(FlintZZ, n)
  for k in 1:10
    i = rand(1:(n - 1))
    j = rand((i + 1):n)
    C[i, j] += rand(B)
  end
  @assert det(C) == 1
  return transpose(C) * A
end

const lattices_and_aut_order = [
(([[2]]), 2),
# 2
(([[1, 0], [0, 2]]), 4),
(([[2, -1], [-1, 2]]), 12),
# 3
(([[2, 1, 0],
   [1, 2, 0],
   [0, 0, 26]]), 24),
# 4
(([[1, 0, 0, 0],
   [0, 2, -1, 1],
   [0, -1, 3, -1],
   [0, 1, -1, 3]]), 16),
# 5
(([[2, 1, 1, 1, -1],
   [1, 2, 1, 1, 0],
   [1, 1, 2, 1, -1],
   [1, 1, 1, 2, -1],
   [-1, 0, -1, -1, 2]]), 3840),
(([[1, 0, 0, 0, 0],
   [0, 1, 0, 0, 0],
   [0, 0, 1, 0, 0],
   [0, 0, 0, 2, 1],
   [0, 0, 0, 1, 3]]), 192),
# 6
(([[2, -1, 0, 0, 0, 0],
   [-1, 2, -1, 0, 0, 0],
   [0, -1, 2, -1, 0, -1],
   [0, 0, -1, 2, -1, 0],
   [0, 0, 0, -1, 2, 0],
   [0, 0, -1, 0, 0, 2]]), 103680),
# 6
(([[1, 0, 0, 0, 0, 0],
   [0, 1, 0, 0, 0, 0],
   [0, 0, 2, 1, 0, 1],
   [0, 0, 1, 3, 1, 1],
   [0, 0, 0, 1, 2, 1],
   [0, 0, 1, 1, 1, 3]]), 512),
#(([[8, 4, 4, 4, 4, 4, 4, 2, 4, 4, 4, 2, 4, 2, 2, 2, 4, 2, 2, 2, 0, 0, 0, -3],
#   [4, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 2, 1, 1, 2, 1, 0, 0, -1],
#   [4, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1, 1, 0, 0, -1],
#   [4, 2, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 2, 2, 1, 2, 1, 1, 0, 0, -1],
#   [4, 2, 2, 2, 4, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 1, 0, 0, -1],
#   [4, 2, 2, 2, 2, 4, 2, 2, 2, 2, 2, 1, 2, 2, 1, 1, 2, 1, 2, 1, 0, 0, 0, -1],
#   [4, 2, 2, 2, 2, 2, 4, 2, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, 2, 0, 0, 0, -1],
#   [2, 2, 2, 2, 2, 2, 2, 4, 1, 1, 1, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 0, 0, 1],
#   [4, 2, 2, 2, 2, 2, 2, 1, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, -1],
#   [4, 2, 2, 2, 2, 2, 2, 1, 2, 4, 2, 2, 2, 2, 1, 1, 2, 2, 1, 1, 0, 1, 0, -1],
#   [4, 2, 2, 2, 2, 2, 2, 1, 2, 2, 4, 2, 2, 1, 2, 1, 2, 1, 2, 1, 0, 0, 1, -1],
#   [2, 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 4, 1, 2, 2, 2, 1, 2, 2, 2, 2, 1, 1, 1],
#   [4, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 4, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, -1],
#   [2, 2, 1, 1, 2, 2, 1, 2, 2, 2, 1, 2, 2, 4, 2, 2, 1, 2, 2, 2, 2, 2, 1, 1],
#   [2, 1, 2, 1, 2, 1, 2, 2, 2, 1, 2, 2, 2, 2, 4, 2, 1, 2, 2, 2, 2, 1, 2, 1],
#   [2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 2, 2, 4, 1, 2, 2, 2, 2, 1, 1, 1],
#   [4, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 1, 1, 1, 4, 2, 2, 2, 1, 1, 1, -1],
#   [2, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 4, 2, 2, 2, 2, 1, 1],
#   [2, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 4, 2, 2, 1, 2, 1],
#   [2, 2, 1, 1, 2, 1, 2, 2, 2, 1, 1, 2, 1, 2, 2, 2, 2, 2, 2, 4, 2, 1, 1, 1],
#   [0, 1, 1, 1, 1, 0, 0, 2, 1, 0, 0, 2, 1, 2, 2, 2, 1, 2, 2, 2, 4, 2, 2, 2],
#   [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 4, 2, 2],
#   [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 2, 4, 2],
#   [-3, -1, -1, -1, -1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 2,
#    2, 2, 4]]), 8315553613086720000),
(([[2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [-1, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, -1, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1],
   [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2]]), 711374856192000)]

function test_automorphisms(L, G, ambient_representation)
  if ambient_representation
    @test all(g * gram_matrix(ambient_space(L)) * transpose(g)  == gram_matrix(ambient_space(L)) for g in G)
  else
    @test all(g * gram_matrix(L) * transpose(g) == gram_matrix(L) for g in G)
  end
end

@testset "Zlattices" begin

  # Creation

  G = matrix(QQ, 2, 2, [2, 1, 1, 2])
  B = matrix(ZZ, 1, 2, [1, 0])
  @test (@inferred Zlattice(B ;gram = G)) isa ZLat
  @test (@inferred Zlattice(B)) isa ZLat
  B = matrix(QQ, 1, 2, [1, 0])
  @test (@inferred Zlattice(B ;gram = G)) isa ZLat
  @test (@inferred Zlattice(B)) isa ZLat

  V = quadratic_space(FlintQQ, G)
  B = matrix(ZZ, 1, 2, [1, 0])
  @test (@inferred lattice(V, B)) isa ZLat
  Lr1 = lattice(V, B)
  B = matrix(QQ, 1, 2, [1, 0])
  @test (@inferred lattice(V, B)) isa ZLat

  B = matrix(GF(2), 1, 2, [1, 0])
  @test_throws ErrorException lattice(V, B)

  B = matrix(QQ, 1, 2, [1, 0])
  Lr1 = lattice(V, B)
  B = matrix(ZZ, 2, 2, [2, 0, 0, 2])
  Lr2 = lattice(V, B)
  B = matrix(ZZ, 0, 2, [])
  Lr0 = lattice(V, B)

  @test (@inferred base_ring(Lr0)) isa FlintIntegerRing

  @test !(@inferred issublattice(Lr2, Lr1))
  M = Zlattice(;gram = FlintQQ[2 2; 2 2])
  @test !(@inferred issublattice(Lr0, M))
  @test issublattice(Lr2, Lr0)
  @test issublattice(Lr1, lattice(V, QQ[2 0;]))

  fl, rels = @inferred issublattice_with_relations(Lr1, lattice(V, QQ[2 0;]))
  @test fl
  @test rels == QQ[2;]

  # lattices of rank 0

  B = matrix(QQ, 0, 2, [])
  @test (@inferred lattice(V, B)) isa ZLat

  # Gram matrix

  @test zero_matrix(QQ, 0, 0) == @inferred gram_matrix(Lr0)
  @test zero_matrix(QQ, 0, 0) == @inferred gram_matrix_of_rational_span(Lr0)

  @test matrix(QQ, 1, 1, [2]) == @inferred gram_matrix(Lr1)
  @test matrix(QQ, 1, 1, [2]) == @inferred gram_matrix_of_rational_span(Lr1)

  B = matrix(ZZ, 2, 2, [2, 0, 0, 2])
  L = lattice(V, B)
  @test matrix(QQ, 2, 2, [8, 4, 4, 8]) == @inferred gram_matrix(Lr2)
  @test matrix(QQ, 2, 2, [8, 4, 4, 8]) == @inferred gram_matrix(Lr2)

  # rational span

  for L in [Lr0, Lr1, Lr2]
    V = @inferred rational_span(L)
    @test gram_matrix(V) == gram_matrix(L)
  end

  # orthogonal sum

  for L in [Lr0, Lr1, Lr2]
    for LL in [Lr0, Lr1, Lr2]
      LLL, = @inferred orthogonal_sum(L, LL)
      @test gram_matrix(LLL) == diagonal_matrix(gram_matrix(L), gram_matrix(LL))
    end
  end

  # printing

  s = sprint(show, "text/plain", Lr0)
  @test occursin("lattice", s)

  # isometry testing
  C1 = root_lattice(:A, 2)
  C1m = rescale(C1,-1)
  @test isisometric(C1m, C1m)[1]
  # automorphisms
  C2 = (1//3)*C1

  for (m, o) in lattices_and_aut_order
    n = length(m[1])
    G = matrix(FlintZZ, n, n, reduce(vcat, m))
    L = Zlattice(gram = G)
    Ge = automorphism_group_generators(L, ambient_representation = true)
    test_automorphisms(L, Ge, true)
    Ge = automorphism_group_generators(L, ambient_representation = false)
    test_automorphisms(L, Ge, false)
    @test automorphism_group_order(L) == o

    # L as a non-full lattice
    for C in [C1,C2]
      if rank(L) > 4 # small examples suffice
        continue
      end
      C1L, _, f = orthogonal_sum(C1, L)
      V = ambient_space(C1L)
      imL = lattice(V, f.matrix)
      Ge = automorphism_group_generators(imL, ambient_representation = true)
      test_automorphisms(imL, Ge, true)
      Ge = automorphism_group_generators(imL, ambient_representation = false)
      test_automorphisms(imL, Ge, false)
      @test automorphism_group_order(L) == o
      @test isisometric(imL, imL)[1]
    end
  end

  D = lattice_database()
  l = number_of_lattices(D)
  for i in 1:min(l, 200)
    L = lattice(D, i)
    L = Zlattice(gram = gram_matrix(L)) # to avoid caching
    Ge = automorphism_group_generators(L, ambient_representation = true)
    test_automorphisms(L, Ge, true)
    Ge = automorphism_group_generators(L, ambient_representation = false)
    test_automorphisms(L, Ge, false)
    @test automorphism_group_order(L) == lattice_automorphism_group_order(D, i)
  end

  # isometry

  for (m, o) in lattices_and_aut_order
    n = length(m[1])
    G = matrix(FlintZZ, n, n, reduce(vcat, m))
    L = Zlattice(gram = G)
    X = _random_invertible_matrix(n, -3:3)
    @assert abs(det(X)) == 1
    L2 = Zlattice(gram = X * G * transpose(X))
    b, T = isisometric(L, L2, ambient_representation = false)
    @test b
    @test T * gram_matrix(L2) * transpose(T) == gram_matrix(L)
    L2 = Zlattice(X, gram = G)
    b, T = isisometric(L, L2, ambient_representation = false)
    @test b
    @test T * gram_matrix(L2) * transpose(T) == gram_matrix(L)
    b, T = isisometric(L, L2, ambient_representation = true)
    @test b
    @test T * gram_matrix(ambient_space(L2)) * transpose(T) ==
    gram_matrix(ambient_space(L))
  end

  D = lattice_database()
  l = number_of_lattices(D)
  for i in 1:min(l, 200)
    L = lattice(D, i)
    L = Zlattice(gram = gram_matrix(L)) # to avoid caching
    n = rank(L)
    X = change_base_ring(FlintQQ, _random_invertible_matrix(n, -3:3))
    @assert abs(det(X)) == 1
    L2 = Zlattice(gram = X * gram_matrix(L) * transpose(X))
    b, T = isisometric(L, L2, ambient_representation = false)
    @test b
    @test T * gram_matrix(L2) * transpose(T) == gram_matrix(L)
  end

  #discriminant of a lattice
  L = Zlattice(ZZ[1 0; 0 1], gram = matrix(QQ, 2,2, [2, 1, 1, 2]))
  @test discriminant(L) == 3

  G = matrix(ZZ, 2, 2, [2, 1, 1, 2])
  L = Zlattice(gram=G)
  @test norm(L)==2
  G = (1//4)*matrix(QQ, 2, 2, [2, 1, 1, 2])
  L = Zlattice(gram=G)
  @test norm(L)==1//2
  G = matrix(ZZ, 0, 0, [])
  L = Zlattice(gram=G)
  @test norm(L) == 0
  @test scale(L) == 0

  #orthogonal submodule of a lattice
  V = quadratic_space(QQ, QQ[1 0 0; 0 1 0; 0 0 1])
  L = lattice(V, ZZ[1 -1 0; 0 1 -1])
  S = lattice(V, ZZ[1 -1 0;])
  submod = Hecke.orthogonal_submodule(L, S)
  @test  basis_matrix(submod) == matrix(QQ, 1, 3, [1 1 -2])

  @test isdefinite(L)
  @test ispositive_definite(L)
  @test !isnegative_definite(L)
  @test L+L == L
  @test intersect(L,L) == L
  @test 2*L == L*2
  @test 0*L == L*0
  @test (1//2)L*2 == L
  @test !(L == 2*L)

  gram = QQ[-2 1 0 0 0 0 0 0 0 0; 1 -2 1 1 0 0 0 0 0 0; 0 1 -2 1 0 0 0 0 0 0; 0 1 1 -2 0 0 0 0 0 0; 0 0 0 0 -2 1 0 0 0 0; 0 0 0 0 1 -2 1 0 0 0; 0 0 0 0 0 1 -2 1 0 1; 0 0 0 0 0 0 1 -2 1 0; 0 0 0 0 0 0 0 1 -2 0; 0 0 0 0 0 0 1 0 0 -2]
  BS = QQ[1 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0]
  BH = QQ[3 12 11 11 0 0 0 0 0 0]
  V = quadratic_space(QQ,gram)
  L = lattice(V, BS)
  H = lattice(V, BH)
  R = Hecke.orthogonal_submodule(L,H)
  @test issublattice(L,R)

  # local modification
  L = rescale(Hecke.root_lattice(:A,3),15)
  M = Hecke.maximal_integral_lattice(L)
  for p in prime_divisors(ZZ(discriminant(L)))
    M = Hecke.local_modification(M, L, p)
  end
  @test genus(M) == genus(L)

  # maximal integral lattice
  G = ZZ[3 0 0 0 0 0 0 0 0 0 0 0;
         0 18 0 0 0 0 0 0 0 0 0 0;
         0 0 9 0 0 0 0 0 0 0 0 0;
         0 0 0 9 0 0 0 0 0 0 0 0;
         0 0 0 0 9 0 0 0 0 0 0 0;
         0 0 0 0 0 9 0 0 0 0 0 0;
         0 0 0 0 0 0 9 0 0 0 0 0;
         0 0 0 0 0 0 0 9 0 0 0 0;
         0 0 0 0 0 0 0 0 9 0 0 0;
         0 0 0 0 0 0 0 0 0 9 0 0;
         0 0 0 0 0 0 0 0 0 0 9 0;
         0 0 0 0 0 0 0 0 0 0 0 54]
  L = Zlattice(gram = G)
  LL = @inferred Hecke.maximal_integral_lattice(L)
  @test isone(norm(LL))

  G = QQ[1 0 0 0; 0 2 0 0; 0 0 17 0; 0 0 0 6]
  V = quadratic_space(QQ, G)
  B = QQ[2 0 0 0; 1 1 0 0; 1 0 1 0; 1//2 1//4 1//2 1//4]
  L = lattice(V, B)
  Ld = dual(L)
  @test issublattice(Ld,L)
  discriminant_group(L)

  # Kernel lattice
  L = root_lattice(:A, 2)
  f = matrix(ZZ, 2, 2, [0, 1, 1, 0])
  L1 = lattice(ambient_space(L), ZZ[1 1])
  L2 = lattice(ambient_space(L), ZZ[1 -1])
  L3 = lattice(ambient_space(L), ZZ[1 0])
  M = kernel_lattice(L, f - 1)
  @test basis_matrix(M) == QQ[1 1;]
  M = kernel_lattice(L, f - 1, ambient_representation = false)
  @test basis_matrix(M) == QQ[1 1;]
  M1 = kernel_lattice(L1, f - 1)
  M2 = kernel_lattice(L2, f - 1)
  @test_throws ArgumentError M3 = kernel_lattice(L3, f - 1)
  @test rank(M1) == 1
  @test rank(M2) == 0
  g = matrix(ZZ,1,1,[1])
  M0 = kernel_lattice(L1,g, ambient_representation =false)
  @test rank(M0) == 0

  f = matrix(QQ, 2, 2, [0, 1, 1, 0])
  M = kernel_lattice(L, f - 1)
  @test basis_matrix(M) == QQ[1 1;]
  M = kernel_lattice(L, f - 1, ambient_representation = false)
  @test basis_matrix(M) == QQ[1 1;]

  L = Zlattice(QQ[1 0; 0 2])
  f = matrix(QQ, 2, 2, [0, 1, 0, 0])
  @test_throws ErrorException kernel_lattice(L, f)
  M = kernel_lattice(L, f, ambient_representation = false)
  @test basis_matrix(M) == QQ[0 2;]

  @test_throws ErrorException root_lattice(:F,3)
  @test_throws ErrorException root_lattice(:D,1)

  L = root_lattice(:A, 2)
  @test local_basis_matrix(L, 2) == 1
  @test local_basis_matrix(L, ideal(ZZ,2)) == 1
  @test det(L) == 3
  G = automorphism_group_generators(L)
  N = invariant_lattice(L, G)
  @test ambient_space(N) === ambient_space(L)
  @test rank(N) == 0
  @test basis_matrix(invariant_lattice(L, identity_matrix(QQ, 2))) == basis_matrix(L)

  randlist = rand(2:20,10)
  L = [root_lattice(:D,i) for i in randlist]
  @test any(l -> discriminant(l) == 4, L)
  @test any(iseven, L)
  @test any(l -> norm(l) == 2, L)

  @test root_lattice(:D, 3) != root_lattice(:A, 2)
  # relies on caching
  @test root_lattice(:D, 3) == root_lattice(:A, 3)

  L = root_lattice(:D,4)
  @test norm(L) == 2
  @test automorphism_group_order(L) == 1152

  L = root_lattice(:E,6)
  @test discriminant(L) == 3
  @test iseven(L)
  @test norm(L) == 2
  @test Hecke.kissing_number(L) == 72

  L = root_lattice(:E,7)
  @test discriminant(L) == 2
  @test iseven(L)
  @test norm(L) == 2
  @test Hecke.kissing_number(L) == 126

  L = root_lattice(:E, 8)
  @test discriminant(L) == 1
  @test iseven(L)
  @test norm(L) == 2
  @test norm(L) == 2 # tests caching

  q = quadratic_space(QQ, QQ[2 1; 1 2])
  L = lattice(q, QQ[0 0; 0 0], isbasis = false)
  g = automorphism_group_generators(L)
  @test rank(L) == 0
  @test g == [identity_matrix(QQ, 2)]

  # membership check
  G = QQ[1 0 0 0; 0 2 0 0; 0 0 17 0; 0 0 0 6]
  V = quadratic_space(QQ, G)
  B = QQ[2 0 0 0; 1 1 0 0; 1 0 1 0; 1//2 1//4 1//2 1//4]
  L = lattice(V, B)
  x1 = [27//11, 1, 1//7, 2]
  x2 = [2//1, 14//2, 5//1, 9//3]
  x4 = [2, 1, 0, 1, 2]
  v = [1//2]
  l = Zlattice(matrix(QQ,1,1,[1//2;]))
  @test !(x1 in L)
  @test !(x2 in L)
  @test B[1,:] in L
  @test [B[4,i] for i in 1:ncols(B)] in L
  @test_throws AssertionError x4 in L
  @test v in l
end


text:
mutable struct Simulation{T <: EpidemicModel}
  time::Float64
  iterations::Int64
  population::Population
  risk_functions::RiskFunctions{T}
  risk_parameters::RiskParameters{T}
  disease_states::Vector{DiseaseState}
  transmission_rates::TransmissionRates
  event_rates::EventRates{T}
  events::Events{T}
  transmission_network::TransmissionNetwork

  function Simulation(pop::Population,
                      rf::RiskFunctions{T},
                      rp::RiskParameters{T}) where T <: EpidemicModel
    states = fill(State_S, pop.individuals)
    tr = initialize(TransmissionRates, states, pop, rf, rp)
    rates = initialize(EventRates, tr, states, pop, rf, rp)
    events = Events{T}(pop.individuals)
    net = TransmissionNetwork(states)
    return new{T}(0.0, 0, pop, rf, rp, states, tr, rates, events, net)
  end

  function Simulation(pop::Population,
                      states::Vector{DiseaseState},
                      rf::RiskFunctions{T},
                      rp::RiskParameters{T};
                      time::Float64 = 0.0,
                      skip_checks::Bool=false) where T <: EpidemicModel
    @debug "Initializing $T Simulation with the following starting states:" states
    if !skip_checks
      if length(states) != pop.individuals
        @error "Length of initial disease state vector must match number of individuals"
      elseif !all(in.(states, Ref(_state_progressions[T])))
        @error "All states in initial disease state vector must be valid within specified epidemic model"
      end
    end
    tr = initialize(TransmissionRates, states, pop, rf, rp)
    rates = initialize(EventRates, tr, states, pop, rf, rp)
    events = Events{T}(states)
    net = TransmissionNetwork(states)
    return new{T}(time, 0, pop, rf, rp, copy(states), tr, rates, events, net)
  end
end

function Base.show(io::IO, x::Simulation{T}) where T <: EpidemicModel
  y = "$T epidemic simulation @ time = $(round(x.time, digits=2))\n"
  for s in _state_progressions[T]
    y *= "\n$s = $(sum(x.disease_states .== Ref(s)))"
  end
  return print(io, y)
end

ext: jl


text:
struct ThreadedProjMPOSum
  pm::Vector{ProjMPO}
end

ThreadedProjMPOSum(mpos::Vector{MPO}) = ThreadedProjMPOSum([ProjMPO(M) for M in mpos])

ThreadedProjMPOSum(Ms::MPO...) = ThreadedProjMPOSum([Ms...])

nsite(P::ThreadedProjMPOSum) = nsite(P.pm[1])

Base.length(P::ThreadedProjMPOSum) = length(P.pm[1])

function product(P::ThreadedProjMPOSum, v::ITensor)::ITensor
  Pvs = fill(ITensor(), Threads.nthreads())
  Threads.@threads for M in P.pm
    Pvs[Threads.threadid()] += product(M, v)
  end
  return sum(Pvs)
end

function Base.eltype(P::ThreadedProjMPOSum)
  elT = eltype(P.pm[1])
  for n in 2:length(P.pm)
    elT = promote_type(elT, eltype(P.pm[n]))
  end
  return elT
end

(P::ThreadedProjMPOSum)(v::ITensor) = product(P, v)

Base.size(P::ThreadedProjMPOSum) = size(P.pm[1])

function position!(P::ThreadedProjMPOSum, psi::MPS, pos::Int)
  Threads.@threads for M in P.pm
    position!(M, psi, pos)
  end
  return P
end

function noiseterm(P::ThreadedProjMPOSum, phi::ITensor, dir::String)
  nts = fill(ITensor(), Threads.nthreads())
  Threads.@threads for M in P.pm
    nts[Threads.threadid()] += noiseterm(M, phi, dir)
  end
  return sum(nts)
end


text:
using Revise

using Flux, Zygote
using Flux.Data.MNIST
using Flux: @epochs
using Flux: throttle
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated
using Statistics

using Images: channelview
using Base.Iterators: partition

using CuArrays

include("data.jl")
include("model.jl")
include("train.jl")

function test_BatchNorm()
    using Flux
    using CuArrays
    bn = BatchNorm(32);
    testimg = randn(Float32, 28, 28, 32, 1);
    @time bn(testimg);
    @time gpu(bn)(gpu(testimg));

    Flux.testmode!(bn)
    Flux.testmode!(bn, false)
    sum(bn(testimg))

    p = Flux.param(testimg);

    gs = Tracker.gradient(()->sum(bn(p)), Flux.params(bn))
    gs

    function f()
        sum(bn(p))
    end
    Flux.gradient(f, p)
end

function test_allowscalar()
    CuArrays.allowscalar(false)
    x = rand(Float32, 10, 3) |> gpu;
    y = Flux.onehotbatch(1:3, 1:10) |> gpu;
    accuracy(x, y) = Flux.onecold(x) .== Flux.onecold(y);
    accuracy(x, y)
    CuArrays.allowscalar(true)
end

function test_bn_grad()
    x = randn(Float32, 32, 32, 3, 16) |> gpu
    model = Chain(
        myConv((3,3), 3=>16, stride=(1,1), pad=1),
        BatchNorm(16),
        MeanPool((8,8)),
        x -> reshape(x, :, size(x)[end]),
        Dense(256,10)) |> gpu
    model(x)

    Flux.gradient() do
        sum(model(x))
    end

    Flux.testmode!(model)

    # In test mode, it would fail if I don't patch Flux.jl/src/cuda/cudnn.jl
    # cudnnBNBackward!. The root cause is sum(iter; dims) have dims as explicit
    # kwarg, and thus calling without dims=... will call the wrong sum. Also
    # squeeze is deprecated for dropdims. NOTE: this is in Flux v0.9, the master
    # branch currently removed manual testmode, see discussions:
    #
    # - Train/test mode discussion by MikeInnes https://github.com/FluxML/Flux.jl/issues/643
    # - https://github.com/FluxML/Flux.jl/issues/909
    # - https://github.com/FluxML/Flux.jl/issues/232
    Flux.gradient() do
        sum(model(x))
    end
end

# TODO learning rate?
function test_CIFAR10_ds()
    ds, test_ds = load_CIFAR10_ds(batch_size=128);
    x, y = next_batch!(ds) |> gpu;

    # model = get_CIFAR_CNN_model()
    model = resnet(20) |> gpu
    # model = resnet(32)
    # model = resnet(56)
    # model = resnet(68)

    # model = WRN(28, 10) |> gpu;
    model = WRN(16, 4) |> gpu;

    model(x)

    adv = attack_CIFAR10_PGD_k(7)(model, x, y);
    Flux.testmode!(model)
    adv = attack_CIFAR10_PGD_k(7)(model, x, y);

    size(adv)

    augment = Augment()
    augx = augment(cpu(x)) |> gpu;
    # augment(x)
    model(augx)
    accuracy_with_logits(model(augx), y)

    opt = ADAM(1e-3)

    train!(model, opt, ds, print_steps=50)

    accuracy_with_logits(model(x), y)

end

function test_free_train()
    ds, test_ds = load_MNIST_ds(batch_size=50);
    x, y = next_batch!(ds) |> gpu;

    model = get_Madry_model()
    model(x)
    opt = ADAM(1e-3)

    size(gpu(zeros(size(x)...)) + x)

    free_train!(model, opt, ds, 0.3)
end

function evaluate_AE(ae, ds; cnn=nothing)
    xx,yy = next_batch!(ds) |> gpu

    @info "Sampling BEFORE images .."
    sample_and_view(xx, yy, cnn)
    @info "Sampling AFTER images .."
    # NOTE: I need .data, currently handled in sample_and_view
    rec = ae(xx)
    sample_and_view(rec, yy, cnn)
end

function test_CIFAR_ae()
    ds, test_ds = load_CIFAR10_ds(batch_size=128);
    x, y = next_batch!(ds) |> gpu;
    x, y = next_batch!(test_ds) |> gpu;

    # get a pretrained WRN
    cnn = CIFAR10_pretrained_fn()
    model = CIFAR10_pretrained_fn()

    accuracy_with_logits(cnn(x), y)
    accuracy_with_logits(model(x), y)
    Flux.testmode!(model)
    Flux.testmode!(model, false)
    Flux.testmode!(cnn)
    Flux.testmode!(cnn, false)

    params(cnn[2].μ)

    cnn(x)

    onecold(cpu(cnn(x)))

    model[2].μ
    cnn[2].μ
    cnn[2].σ²

    cnn

    train!(cnn, opt, ds)

    @save "test.bson" MMM=cpu(cnn) opt
    @load "test.bson" MMM

    MMM=gpu(MMM);

    fieldnames(typeof(cnn[5]))
    cnn[5].norm_layers[1].σ²

    # get the AE

    # ae = cifar10_AE()
    ae = dunet()
    encoder, decoder = cifar10_deep_AE()

    ae = Chain(encoder, decoder)

    rx = randn(Float32, 32, 32, 3, 16) |> gpu;
    size(encoder(rx))
    size(decoder(encoder(rx)))
    size(ae(rx))

    opt = ADAM(1e-3)
    # train the AE
    aetrain!(ae, opt, ds)

    evaluate_AE(ae, test_ds, cnn=cnn)

    param_count(ae)
    param_count(cnn)
end

function test_AE()
    ds, test_ds = load_MNIST_ds(batch_size=50);
    x, y = next_batch!(ds) |> gpu;

    cnn = get_Madry_model()
    @show size(cnn(x))

    opt = ADAM(1e-3)

    Flux.testmode!(cnn, false)
    cnn = MNIST_pretrained_fn()

    accuracy_with_logits(cnn(x), y)

    cnn(x)

    size(cnn[1](x))

    # TODO schedule lr simply by setting field of opt, and maintain the
    # state. FIXME But it seems that the state is maintained by a dict with key
    # Param(x). I'm not sure if this aproach can resume the state.
    #
    # opt.eta = 1e-4

    train!(cnn, opt, ds)

    cnn = maybe_train("trained/pretrain-MNIST.bson", get_Madry_model(), train_fn)

    # ae = dense_AE()
    ae = CNN_AE()
    # ae = CNN2_AE()
    size(ae(x))

    # FIXME opt states?
    opt = ADAM(1e-3)
    opt = ADAM(1e-4)
    # FIXME Flux.mse(logits, x)
    # FIXME sigmoid mse?
    aetrain!(ae, opt, ds)

    evaluate_AE(ae, test_ds, cnn=cnn)

    advae_train!(ae, cnn, opt, attack_PGD_k(40), ds, print_steps=20)

    evaluate(cnn, test_ds)
    evaluate(cnn, test_ds, attack_fn=attack_PGD_k(40))

    evaluate(Chain(ae, cnn), test_ds)
    evaluate(Chain(ae, cnn), test_ds, attack_fn=attack_FGSM)
    evaluate(Chain(ae, cnn), test_ds, attack_fn=attack_PGD_k(40))
end

function test_log_image()
    logger = TBLogger("tensorboard_logs/test/image", tb_append, min_level=Logging.Info)

    ds, test_ds = load_CIFAR10_ds(batch_size=128);
    x, y = next_batch!(ds) |> gpu;
    img = x[:,:,:,1];

    # CAUTION if pass in raw img, the error message contains the whole array,
    # which is very slow
    showable("image/png", MyImage(img))
    log_image(logger, "testimage", MyImage(img));

    display(MyImage(img))
end

function sha_arr(x)
    # assuming 3 channel images in a batch
    res = map(collect(1:size(x)[end])) do i
        img = x[:,:,:,i]
        bytes2hex(sha256(repr(img)))
    end
    res
end

function get_wrong_sha(model, ds)
    # evaluate model, and get the wrong sample shas
    res = []
    @showprogress for i in 1:ds.nbatch
        x, y = next_batch!(ds) |> gpu
        logits = model(x)
        wrong_idx = onecold(cpu(logits)) .!= onecold(cpu(y))
        if length(wrong_idx) > 0
            # get sha
            shas = sha_arr(cpu(x)[:,:,:,wrong_idx])
            push!(res, shas...)
        end
    end
    res
end

function load_trained_model(file)
    @load file model
    return gpu(model)
end

# https://github.com/FluxML/Flux.jl/issues/160
function weight_params(m::Chain, ps=Flux.Params())
    map((l)->weight_params(l, ps), m.layers)
    ps
end
weight_params(m::Dense, ps=Flux.Params()) = push!(ps, m.W)
weight_params(m::Conv, ps=Flux.Params()) = push!(ps, m.weight)
weight_params(m::ConvTranspose, ps=Flux.Params()) = push!(ps, m.weight)
weight_params(m, ps=Flux.Params()) = ps

function test_weight_params()
    # get model
    weight_params(model)
end

function test_boundary()
    # test whether all the nail households
    model = CIFAR10_pretrained_fn()
    ds, test_ds = load_CIFAR10_ds(batch_size=128);

    itM = load_trained_model(
        "trained/CIFAR10/itadv-schedule/(0.0004=>40000,4.0e-5=>60000,4.0e-6=>80000).bson")
    dyM = load_trained_model(
        "trained/CIFAR10-dyattack/(0=>2000,1=>3000,2=>4000,3=>5000,4=>6000,7=>8000).bson")
    freeM = load_trained_model("trained/CIFAR10-free/test-2019-11-18T16:57:55.947.bson")

    sha_it = get_wrong_sha(itM, test_ds);
    sha_dy = get_wrong_sha(dyM, test_ds);
    sha_free = get_wrong_sha(freeM, test_ds);

    size(sha_it)
    size(unique(sha_it))
    @show length(sha_it) length(sha_dy) length(sha_free)
    @show length(unique(sha_it)) length(unique(sha_dy)) length(unique(sha_free))
    @show size(intersect(sha_it, sha_dy, sha_free))
    @show size(intersect(sha_it, sha_dy, sha_free))

    ps = params(model);

    model[1].weight in ps

    model

    Flux.mapparams(model) do p
        @show typeof(p)
        p
    end

    model(x)
    typeof(ps)
    typeof(keys(ps.params.dict))
    keytype(ps.params.dict)

    shas[1]
    size(shas)
    size(unique(shas))
    size(unique(shas2))
    # overlap
    size(intersect(shas, shas2))
end

function test()
    ds, test_ds = load_CIFAR10_ds(batch_size=128);
    x, y = next_batch!(ds) |> gpu;
    x, y = next_batch!(test_ds) |> gpu;
    model = load()
    model[2].μ
    # Flux.testmode!(model, false)
    model(x)

    # hash the x
    img = x[:,:,:,1];
    size(img)

    # using SHA
    @time bytes2hex(sha256("test"))

    @time bytes2hex(sha256(repr(img)))

    cat(res..., dims=4)
    size(res)
    res[1]
    size(unique(res))
    res[1]

    ss = map(x) do img
        bytes2hex(sha256(repr(img)))
    end

    for i in 1:10
        i
    end

    accuracy_with_logits(model(x), y)

    sum(onecold(cpu(model(x))) .== onecold(cpu(y)))

    adv = attack_CIFAR10_PGD_k(7)(model, x, y);

    view_preds(model, x, y)
    view_preds(model, adv, y)

    sample_and_view(wrong, cpu(y)[])
    sample_and_view(right)
    displayable("image/png")
end

function view_preds(model, x, y)
    @show accuracy_with_logits(model(x), y)
    right_idx = onecold(cpu(model(x))) .== onecold(cpu(y));
    wrong_idx = onecold(cpu(model(x))) .!= onecold(cpu(y));
    @info "right predictions:"
    sample_and_view(cpu(x)[:,:,:,right_idx], cpu(y)[:,right_idx], model)
    @info "wrong predictions:"
    sample_and_view(cpu(x)[:,:,:,wrong_idx], cpu(y)[:,wrong_idx], model)
    nothing
end

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
__precompile__()
module CurveDiscretization

using Compat
using Compat.Printf

using FastGaussQuadrature
using FMMLIB2D

export DiscreteCurve
export discretize
export multiply_connected
export traparclen
export minmax_panel
export interior_points
export nearest_nb_R

# Curve discretization
struct DiscreteCurve
    panelorder::Int64
    numpanels::Int64
    numpoints::Int64
    edges::Array{Float64, 2}
    n_edges::Array{Float64, 2}
    t_edges::Array{Float64, 2}
    points::Array{Float64, 2}
    normals::Array{Float64, 2}
    weights::Array{Float64}
    dS::Array{Float64}
    curvature::Array{Float64}
    prevpanel::Array{Int64}
    nextpanel::Array{Int64}
    curvenum::Array{Int64}
end

include("legendre.jl")
include("mappings.jl")

"""
    Discretize curve using composite Gauss-Legendre quadrature
    """
function discretize(curve, numpanels, panelorder;
                    equal_arclength=true)
    if isa(curve, Array)
        # TODO: Would be nicer to use multiple dispatch
        curves = curve
        Ncurves = length(curves)
        @assert length(numpanels)==Ncurves
        dcurves = Array{DiscreteCurve}(Ncurves)
        for i=1:Ncurves
            dcurves[i] = discretize(curves[i], numpanels[i], panelorder, equal_arclength=equal_arclength)
        end
        return multiply_connected(dcurves)
    end
    start = 0.0
    stop = 2*pi
    # Setup division in parameter
    if equal_arclength
        tsplit, _, _ = traparclen(curve, numpanels)
        push!(tsplit, 2*pi)
    else
        # Equal parameter length
        tsplit = linspace(start, stop, numpanels+1);
    end
    # Prepare structures
    numpoints = numpanels*panelorder
    glpoints, glweights = gausslegendre(panelorder)
    points = zeros(2, numpoints)
    normals = zeros(2, numpoints)
    derivatives = zeros(2, numpoints)
    derivatives2 = zeros(2, numpoints)
    weights = zeros(numpoints)
    dS = zeros(numpoints)
    curvature = zeros(numpoints)
    t_edges = zeros(2, numpanels)
    edges = zeros(4, numpanels)
    n_edges = zeros(4, numpanels)
    splitcomplex(z) = [real(z), imag(z)]
    L = legendre_matrix(panelorder)
    resolution = 0.0
    # Populate
    for i in 1:numpanels
        # setup quadrature for t in [a,b]
        a, b = tsplit[i:i+1]
        t_edges[:, i] = [a, b]
        edges[1:2, i] = splitcomplex(curve.tau(a))
        edges[3:4, i] = splitcomplex(curve.tau(b))
        n_edges[1:2, i] = splitcomplex(curve.normal(a))
        n_edges[3:4, i] = splitcomplex(curve.normal(b))
        dt = b-a
        for j=1:panelorder
            idx=j + panelorder*(i-1)
            t = a + (1+glpoints[j])*dt/2
            w = glweights[j]*dt/2
            z = curve.tau(t)
            zp = curve.dtau(t)
            zpp = curve.d2tau(t)
            n = curve.normal(t)
            points[:, idx] = splitcomplex(z)
            derivatives[:, idx] = splitcomplex(zp)
            derivatives2[:, idx] = splitcomplex(zpp)
            normals[:, idx] = splitcomplex(n)
            weights[idx] = w
            dS[idx] = w*abs(zp)
            curvature[idx] = imag(conj(zp)*zpp/abs(zp)^3)
        end
        # Get resolution estimate
        idx = (1:panelorder) + panelorder*(i-1)
        coeff = L*(dS[idx]./weights[idx])
        thisresolution = abs(coeff[end]) + abs(coeff[end-1])
        resolution = max(resolution, thisresolution)
    end

    prevpanel = circshift(1:numpanels, 1)
    nextpanel = circshift(1:numpanels, -1)
    curvenum = ones(Int64, numpoints)

    info(@sprintf("  Grid resolution: %.2e\n", resolution))

    DiscreteCurve(
        panelorder,
        numpanels,
        numpoints,
        edges,
        n_edges,
        t_edges,
        points,
        normals,
        weights,
        dS,
        curvature,
        prevpanel,
        nextpanel,
        curvenum
    )
end

"""
Created a multiply connected curve discretization based on several
simply connected ones.
"""
function multiply_connected(grids)
    # Init
    panelorder = grids[1].panelorder
    numpanels = 0
    numpoints = 0
    points = zeros(2, 0)
    normals = zeros(2, 0)
    derivatives = zeros(2, 0)
    derivatives2 = zeros(2, 0)
    weights = zeros(0)
    dS = zeros(0)
    curvature = zeros(0)
    edges = zeros(4, 0)
    n_edges = zeros(4, 0)
    t_edges = zeros(2, 0)
    prevpanel = Int64[]
    nextpanel = Int64[]
    curvenum = Int64[]
    for i=1:length(grids)
        assert(grids[i].panelorder == panelorder)
        # Concatenate
        edges = hcat(edges, grids[i].edges)
        n_edges = hcat(n_edges, grids[i].n_edges)
        t_edges = hcat(t_edges, grids[i].t_edges)
        points = hcat(points, grids[i].points)
        normals = hcat(normals, grids[i].normals)
        weights = vcat(weights, grids[i].weights)
        dS = vcat(dS, grids[i].dS)
        curvature = vcat(curvature, grids[i].curvature)
        prevpanel = vcat(prevpanel, grids[i].prevpanel+numpanels)
        nextpanel = vcat(nextpanel, grids[i].nextpanel+numpanels)
        append!(curvenum, i*ones(Int64, grids[i].numpoints))
        # Update count
        numpanels += grids[i].numpanels
        numpoints += grids[i].numpoints
    end
    # Construct output
    return DiscreteCurve(
        panelorder,
        numpanels,
        numpoints,
        edges,
        n_edges,
        t_edges,
        points,
        normals,
        weights,
        dS,
        curvature,
        prevpanel,
        nextpanel,
        curvenum
    )
end

function glpanels(a, b, numpanels, pedges, order)
    # Discretize interval using composite Gauss-Legendre quadrature
    #
    # Original code by Rikard Ojala
    T, W = gausslegendre(order)
    t = zeros(order*numpanels)
    w = zeros(order*numpanels)
    ptr = 1
    for j = 1:numpanels
        t[ptr:ptr+order-1] = (pedges[j+1]-pedges[j])*T/2+(pedges[j]+pedges[j+1])/2;
        wj = W*(pedges[j+1]-pedges[j])/2
        w[ptr:ptr+order-1] = wj
        ptr = ptr + order
    end
    return t, w
end

function traparclen(curve, N, NL=1000)
    # t,z,W,L = traparclen(f, N, NL=1000)
    #
    # Discretizes the curve given by f(t) t = [0,2*pi] using the trapezoidal rule
    # equidistant in arc-length .
    # N is the number of quadrature points requested and if specified,
    # NL is the number of 16-point Gauss-Legendre panels used to compute the
    # length of the curve. Default is NL = 1000.
    # Non-adaptive. Make sure that the curve is well resolved by N points.
    #
    # Original code by Rikard Ojala
    if N*2 > NL
        NL = 2*N
        info("Using $NL Gauss-Legendre points in traparclen")
    end
    order = 16
    pedges = linspace(0, 2*pi, NL+1)
    T, W = glpanels(0, 2*pi, NL, pedges, order)
    zp = curve.dtau.(T)
    L = W'*abs.(zp)
    L2 = linspace(1/N,1-1/N,N-1)*L
    #Initial guess
    t = linspace(0,2*pi,N+1)
    t = t[2:end-1]
    dt = 1
    iter = 0
    maxiter = 30
    while norm(dt)/norm(t) > 1e-13 && iter < maxiter
        # Set up n-point GL quadrature on each segment.
        T, W = glpanels(0,2*pi,N-1,[0;t],order)
        zpt = curve.dtau.(t)
        zpT = curve.dtau.(T)
        #Compute the cumulative sum of all the segment lengths
        F = cumsum(sum(reshape(W.*abs.(zpT),order,N-1), 1)')
        dt = (F-L2)./abs.(zpt)
        # Sort the parameters just in case.
        t = vec(t - dt)
        sort!(t)
        iter = iter + 1
    end
    if iter == maxiter
        warn("Newton for equal arclength did not converge, bad discretization")
    end
    t = [0;t]
    z = curve.tau.(t)
    W = 2*pi/N*ones(N,1)
    return t, z, W, L
end

function minmax_panel(grid::DiscreteCurve)
    hmin, hmax = Inf, 0.0
    for i=1:grid.numpanels
        idx = (i-1)*grid.panelorder + (1:grid.panelorder)
        h = sum(grid.dS[idx])
        hmin = min(hmin, h)
        hmax = max(hmax, h)
    end
    return hmin, hmax
end

function interior_points(grid::DiscreteCurve, zt::Array{Float64})
    # First pass: Stokes DLP identity works to O(1) very close to bdry
    ndS = grid.normals .* grid.dS'
    # # Hack: Approximate stokes DLP by modstokes DLP with small alpha
    # density = ones(2, grid.numpoints)
    # u = fmm_stresslet_targ(grid.points, zt, density, ndS, 1e-8)
    # marker = u[1,:]
    # Better, use Laplace DLP
    density = ones(grid.numpoints)
    U = rfmm2d(source=grid.points, target=zt, dipstr=density, dipvec=ndS,
                     ifpot=false, ifpottarg=true)
    marker = -1/(2*pi)*U.pottarg
    interior = marker .> 0.5
    # Second pass: find points close to boundary and distrust them
    _, hmax = minmax_panel(grid)
    dx = hmax*3/grid.panelorder # Estimate of largest point gap on bdry
    R = dx
    nnb, _ = nearest_nb_R(grid, zt, R)
    # listed points are within R of a boundary point
    # Third pass: check sign of projection against near normal
    veryclose = zeros(interior)
    N = size(zt, 2)
    for i=1:N
        if nnb[i] != 0
            rvec = zt[:,i] - grid.points[:,nnb[i]]
            n = grid.normals[:,nnb[i]]
            rdotn = rvec[1]*n[1] + rvec[2]*n[2]
            interior[i] = rdotn > 0
            veryclose[i] = abs(rdotn) < R/10
        end
    end
    # Fourth pass: Find projection using parametrization
    coeffs = map_panels(grid)
    for i=1:N
        if veryclose[i] != 0
            # Find corresponding panel
            near_panel = Int64(ceil(nnb[i]/grid.panelorder))
            z = zt[1,i] + 1im*zt[2,i]
            tloc, _ = invert_map(grid, coeffs, near_panel, z)
            interior[i] = imag(tloc) > 0
        end
    end
    return interior
end

function nearest_nb_R(grid::DiscreteCurve, zt::Array{Float64}, R::Float64)
    # Bin sort of boundary points
    delta = 1e-8
    xmax = max(maximum(zt[1,:]), maximum(grid.points[1,:])) + delta
    xmin = min(minimum(zt[1,:]), minimum(grid.points[1,:])) - delta
    ymax = max(maximum(zt[2,:]), maximum(grid.points[2,:])) + delta
    ymin = min(minimum(zt[2,:]), minimum(grid.points[2,:])) - delta
    Lx = xmax-xmin
    Ly = ymax-ymin
    Nx = Int64(floor(Lx/R))
    Ny = Int64(floor(Ly/R))
    hx = Lx/Nx
    hy = Ly/Ny
    Nx += 1
    Ny += 1
    bin_count, bin_lists = binsort_points(grid.points, xmin, ymin, hx, hy, Nx, Ny)
    # For each target points, find neareast neighbor
    R2 = R^2
    N = size(zt, 2)
    nearest_nb = zeros(Int64, N)
    dist = zeros(N)
    for trg_idx=1:N
        z = zt[:,trg_idx]
        ih, jh = bin_idx(z, xmin, ymin, hx, hy) # home bin
        r2min = Inf
        imin = 0
        # Iterate over nb bins
        for inb = ih+(-1:1)
            for jnb = jh+(-1:1)
                if inb>0 && jnb>0 && inb<Nx+1 && jnb<Ny+1
                    # Compare to target points in bin
                    list = bin_lists[inb, jnb]
                    for k=1:length(list)
                        i = list[k]
                        r2 = (grid.points[1,i]-z[1])^2 + (grid.points[2,i]-z[2])^2
                        if r2<r2min && r2<R2
                            imin = i
                            r2min = r2
                        end
                    end
                end
            end
        end
        nearest_nb[trg_idx] = imin
        dist[trg_idx] = sqrt(r2min)
    end # nearest nb list built
    return nearest_nb, dist
end

function binsort_points(zt, xmin, ymin, hx, hy, Nx, Ny)
    N = size(zt, 2)
    # 1. Count points in each bin
    bin_count = zeros(Int64, Nx, Ny)
    for i=1:N
        ix, iy = bin_idx(zt[:,i], xmin, ymin, hx, hy)
        bin_count[ix, iy] += 1
    end
    # 2. Setup bin lists
    fill_ptr = ones(Int64, Nx, Ny)
    bin_lists = Array{Array{Int64}}(Nx, Ny)
    for i=1:Nx
        for j=1:Ny
            bin_lists[i,j] = zeros(Int64, bin_count[i,j])
        end
    end
    # 3. Fill bins
    for i=1:N
        ix, iy = bin_idx(zt[:,i], xmin, ymin, hx, hy)
        ptr = fill_ptr[ix,iy]
        bin_lists[ix,iy][ptr] = i
        fill_ptr[ix,iy] += 1
    end
    return bin_count, bin_lists
end

@inline function bin_idx(z, xmin, ymin, hx, hy)
    ix = Int64(floor( (z[1]-xmin)/hx )) + 1
    iy = Int64(floor( (z[2]-ymin)/hy )) + 1
    return ix, iy
end

end


text:
zlib_test_data = [
    UInt8[
        0x78, 0x5e,
        0x01, 0x03, 0x00, 0xfc, 0xff, 0x66, 0x6f, 0x6f,
        0x02, 0x82, 0x01, 0x45
    ]
]

@testset "Decompression" begin
    indata = zlib_test_data[1]
    decompressor = Decompressor()
    output = zeros(UInt8, 128)

    @test zlib_decompress!(decompressor, output, indata) == 3
    @test String(output[1:3]) == "foo"
    @test zlib_decompress!(decompressor, output, indata, 3) == 3
    @test String(output[1:3]) == "foo"

    @test zlib_decompress!(decompressor, output, indata, 2) == LibDeflateErrors.deflate_insufficient_space
    @test zlib_decompress!(decompressor, output, indata, 4) == LibDeflateErrors.deflate_output_too_short

    @test zlib_decompress!(decompressor, output[1:2], indata, 3) == LibDeflateErrors.deflate_insufficient_space
    @test zlib_decompress!(decompressor, output[1:2], indata) == LibDeflateErrors.deflate_insufficient_space

    cp = copy(indata)

    cp[1] = 0x79
    @test zlib_decompress!(decompressor, output, cp) == LibDeflateErrors.zlib_not_deflate

    cp[1] = 0x98
    @test zlib_decompress!(decompressor, output, cp) == LibDeflateErrors.zlib_wrong_window_size
    cp[1] = 0x78

    cp[2] = 0xff
    @test zlib_decompress!(decompressor, output, cp) == LibDeflateErrors.zlib_needs_compression_dict

    cp[2] = 0x5c
    @test zlib_decompress!(decompressor, output, cp) == LibDeflateErrors.zlib_bad_header_check

    cp[2] = 0x01
    @test zlib_decompress!(decompressor, output, indata, 3) == 3
    cp[2] = 0xda
    @test zlib_decompress!(decompressor, output, indata, 3) == 3
    cp[2] = 0x9c
    @test zlib_decompress!(decompressor, output, indata, 3) == 3

    cp[end] = 0x46
    @test zlib_decompress!(decompressor, output, cp) == LibDeflateErrors.zlib_bad_adler32
end

@testset "Compression" begin
    output = zeros(UInt8, 128)
    compressor = Compressor()

    @test zlib_compress!(compressor, output, "foo") == length(first(zlib_test_data))
    @test output[1:length(first(zlib_test_data))] == first(zlib_test_data)

    @test zlib_compress!(compressor, zeros(Float64, 4), "foo") == length(first(zlib_test_data))
    @test zlib_compress!(compressor, zeros(Int8, 0), "foo") == LibDeflateErrors.zlib_insufficient_space
    @test zlib_compress!(compressor, zeros(Float64, 1), "foo") == LibDeflateErrors.deflate_insufficient_space

    for (compresslevel, header) in [(1, 0x0178), (4, 0x9c78), (LibDeflate.DEFAULT_COMPRESSION_LEVEL, 0x5e78), (12, 0xda78)]
        compressor = Compressor(compresslevel)
        @test zlib_compress!(compressor, output, "foo") > 6
        @test UInt16(output[1]) | (UInt16(output[2]) << 8) == header
    end
end

@testset "Adler32" begin
    @test adler32("") == UInt32(1)
    @test adler32("", 0x98765432) == 0x98765432
    @test adler32("foo") == 0x02820145
end

@testset "Round trip" begin
    compressed = zeros(UInt8, 128)
    decompressed = zeros(UInt8, 128)
    compressor = Compressor()
    decompressor = Decompressor()
    for input in ("Abracadabra", "", collect(codeunits("Power overwhelming")))
        v = Vector{UInt8}(input)
        n_compressed_bytes = zlib_compress!(compressor, compressed, input)
        @test n_compressed_bytes > 6
        @test zlib_decompress!(decompressor, decompressed, compressed[1:n_compressed_bytes], sizeof(input)) isa Int
        @test transcode(ZlibDecompressor, compressed[1:n_compressed_bytes]) == v

        # Can decompress zlib
        zlib_compressed = transcode(ZlibCompressor, v)
        @test zlib_decompress!(decompressor, decompressed, zlib_compressed) == sizeof(input)
        @test decompressed[1:sizeof(input)] == v
    end
end


text:
data_type_map=Dict{String,Int}("Int"=>11,"Decimal"=>9,"String"=>14,"BigInt"=>1,"Boolean"=>3,"DateTime"=>6,"UniqueIdentifier"=>22,
"Int32"=>11,"Int64"=>1,"Double"=>10,"Guid"=>22)
col_type_map=Dict{String,Int}("Key"=>1,"Value"=>4,"TimeStampKey"=>3,"ForeignKey"=>2)
schema_type_map=Dict{String,Int}("DataIngestion"=>1,"TimeSeries"=>2,"Processed"=>3)

"""
Parameters:
indexname::String,indexcolumnnames::Array{String,1}
Required:
indexname
indexcolumnnames
"""
type DestinationTableDefinitionIndex
    IndexName::String
    #ColumnName:____
    IndexColumns::Array{Dict{String,String},1}
    function DestinationTableDefinitionIndex(;indexname::String=nothing,indexcolumnnames::Array{String,1}=nothing)
        cols=Array{Dict{String,String}}(0)
        for col in indexcolumnnames
            push!(cols,Dict("ColumnName"=>col))
        end
        new(indexname,cols)
    end
end

"""
Parameters:
name::String,datatype::String,columntype::String,isrequired::Bool=false
Required:
name
datatype
columntype
"""
type DestinationTableDefinitionColumn
    ColumnDataType::Int
    ColumnName::String
    ColumnType::Int
    IsRequired::Bool
    IsSystemColumn::Bool
    ValidationError::String
    WasRemoved::Bool
    ColumnDataTypePrecision
    ColumnDataTypeScale
    ColumnDataTypeSize
    ForeignKeyTableName
    ForeignKeyColumnName
    Index::Int
    function DestinationTableDefinitionColumn(;name::String=nothing,datatype::String=nothing,columntype::String=nothing,isrequired::Bool=false)
        col=new()
        col.ColumnDataTypePrecision=nothing
        col.ColumnDataTypeScale=nothing
        col.ColumnDataTypeSize=nothing
        col.ForeignKeyTableName=nothing
        col.ForeignKeyColumnName=nothing

        col.ColumnName=name

        if contains(columntype,"ForeignKey")
            col.ColumnType=col_type_map["ForeignKey"]
            col.ForeignKeyTableName=split(columntype,['(',')',','])[2]
            col.ForeignKeyColumnName=split(columntype,['(',')',','])[3]
        else
            col.ColumnType=col_type_map[columntype]
        end

        col.IsRequired=isrequired
        col.IsSystemColumn=false
        col.ValidationError=""
        col.WasRemoved=false

        if contains(datatype,"String")
            col.ColumnDataType=data_type_map["String"]
            col.ColumnDataTypeSize=parse(split(datatype,['(',')',','])[2])
        elseif contains(datatype,"Decimal")
            col.ColumnDataType=data_type_map["Decimal"]
            col.ColumnDataTypePrecision=parse(split(datatype,['(',')',','])[2])
            col.ColumnDataTypeScale=parse(split(datatype,['(',')',','])[3])
        else
            col.ColumnDataType=data_type_map[datatype]
        end
        return col
    end
end

"""
Parameters:
function DestinationTableDefinition(;allowdatachanges::Bool=false,columns::Array{DestinationTableDefinitionColumn,1}=nothing,
        name::String=nothing,tableindexes::Union{Array{DestinationTableDefinitionIndex,1},Void}=nothing,
        schematype::Union{String,Void}=nothing,schematypeid::Union{Int,Void}=nothing,partitionpath::Union{String,Void}=nothing)
Required:
name
columns
schematype or schematypeid
"""
type DestinationTableDefinition
    DestinationTableDefinitionId::String
    AllowDataLossChanges::Bool
    #CreatedBy::String
    #CreatedDate::String
    DestinationTableDefinitionColumns::Array{DestinationTableDefinitionColumn,1}
    DestinationTableDefinitionIndexes::Array{DestinationTableDefinitionIndex,1}
    DestinationTableName::String
    #LastChangedBy::String
    #LastChangedDate::String
    #MappingsCount::Int
    #SchemaError::Bool
    #StorageType::Int
    DataStoreId::Union{String,Void}
    SchemaType::Int
    FilePath::Union{String,Void}
    function DestinationTableDefinition(;allowdatachanges::Bool=false,columns::Array{DestinationTableDefinitionColumn,1}=nothing,
        name::String=nothing,tableindexes::Union{Array{DestinationTableDefinitionIndex,1},Void}=nothing,
        schematype::Union{String,Void}=nothing,schematypeid::Union{Int,Void}=nothing,partitionpath::Union{String,Void}=nothing)
        for ind=1:length(columns)
            columns[ind].Index=ind
        end
        if schematypeid==nothing
            if schematype=="DataIngestion"
                schematypeid=1
            elseif schematype=="TimeSeries"
                schematypeid=2
            elseif schematype=="Processed"
                schematypeid=3
            else
                error("schematype must be \"DataIngestion\", \"TimeSeries\" or \"Processed\"")
            end
        end
        datastoreid=nothing
        indexes=DestinationTableDefinitionIndex[]
        if tableindexes!=nothing
            indexes=tableindexes
        end
       return new("",allowdatachanges,columns,indexes,name,datastoreid,schematypeid,partitionpath)
    end
end

"create_destination_table(;storename::String=val1,tabledefinition::DestinationTableDefinition=val2)"
function create_destination_table(;storename::String=nothing,tabledefinition::DestinationTableDefinition=nothing)
    datastoreid=""
    try
        datastoreid=NeuroJulia.neurocall("80","datastoremanager","GetDataStores",Dict("StoreName"=>storename))["DataStores"][1]["DataStoreId"]
    catch
        error("Data Store name is not valid")
    end

    tabledefinition.DataStoreId=datastoreid
    NeuroJulia.neurocall("80","datapopulation","CreateDestinationTableDefinition",tabledefinition)
end

type GetDestinationTableDefinitionRequest
    TableName
    DataStoreId
end

"get_table_definition(;storename::String=val1,tablename::String=val2)"
function get_table_definition(;storename::String=nothing,tablename::String=nothing)
    datastores=NeuroJulia.neurocall("80","datastoremanager","GetDataStores",Dict("StoreName"=>storename))["DataStores"]
    if length(datastores)<1
        error("Data store doesn't exist")
    end
    datastoreid=datastores[1]["DataStoreId"]
    request=GetDestinationTableDefinitionRequest(tablename,datastoreid)
    table_def=NeuroJulia.neurocall("80","DataPopulation","GetDestinationTableDefinition",request)

    if length(table_def["DestinationTableDefinitions"])==0
        error("Table doesn't exist")
    end

    cols=NeuroData.DestinationTableDefinitionColumn[]
    for col in table_def["DestinationTableDefinitions"][1]["DestinationTableDefinitionColumns"]
        if col["ColumnName"]!="NeuroverseLastModified"
            tmp_col=NeuroData.DestinationTableDefinitionColumn(name=col["ColumnName"],
                    datatype="Int",
                    columntype="Value",
                    isrequired=col["IsRequired"])

            tmp_col.ColumnDataTypeScale=col["ColumnDataTypeScale"]
            tmp_col.ColumnDataType=col["ColumnDataType"]
            tmp_col.ColumnDataTypePrecision=col["ColumnDataTypePrecision"]
            tmp_col.ColumnDataTypeSize=col["ColumnDataTypeSize"]
            tmp_col.ColumnType=col["ColumnType"]
            push!(cols,tmp_col)
        end
    end
    indexes=DestinationTableDefinitionIndex[]
    for ind in table_def["DestinationTableDefinitions"][1]["DestinationTableDefinitionIndexes"]
        push!(indexes,DestinationTableDefinitionIndex(;indexname=ind["IndexName"],indexcolumnnames=[ind["IndexColumns"][i]["ColumnName"] for i = 1:length(ind["IndexColumns"])]))
    end

    schematypeid=table_def["DestinationTableDefinitions"][1]["SchemaType"]

    schematype=""
    if schematypeid==1
        schematype="DataIngestion"
    elseif schematypeid==2
        schematype="TimeSeries"
    elseif schematypeid==3
        schematype="Processed"
    end

    new_table_def=NeuroData.DestinationTableDefinition(allowdatachanges=table_def["DestinationTableDefinitions"][1]["AllowDataLossChanges"],
    columns=cols,name=table_def["DestinationTableDefinitions"][1]["DestinationTableName"],tableindexes=indexes,schematype=schematype)
    new_table_def.DestinationTableDefinitionId=table_def["DestinationTableDefinitions"][1]["DestinationTableDefinitionId"]
    return new_table_def
end

"add_destination_table_indexes(;storename::String=val1,tablename::String=val2,tableindexes::Array{DestinationTableDefinitionIndex,1}=val3)"
function add_destination_table_indexes(;storename::String=nothing,tablename::String=nothing,tableindexes::Array{DestinationTableDefinitionIndex,1}=nothing)
    datastoreid=nothing
    try
        datastoreid=NeuroJulia.neurocall("80","datastoremanager","GetDataStores",Dict("StoreName"=>storename))["DataStores"][1]["DataStoreId"]
    catch
        error("Data Store name is not valid")
    end
    request=GetDestinationTableDefinitionRequest(tablename,datastoreid)
    table_def=NeuroJulia.neurocall("80","DataPopulation","GetDestinationTableDefinition",request)["DestinationTableDefinitions"][1]
    append!(table_def["DestinationTableDefinitionIndexes"],JSON.parse(JSON.json(tableindexes)))
    NeuroJulia.neurocall("80","datapopulation","UpdateDestinationTableDefinition",table_def)
end

"save_table_definition(;tabledefinition::DestinationTableDefinition=val1,filename::String=val2)"
function save_table_definition(;tabledefinition::DestinationTableDefinition=nothing,filename::String=nothing)
    def=JSON.json(tabledefinition)
    file=open(filename,"w+")
    write(file,def)
    close(file)
end

"load_table_definition(;filename::String=val1)"
function load_table_definition(;filename::String=nothing)
    file=open(filename)
    table_def=JSON.parse(readstring(file))
    close(file)
    cols=NeuroData.DestinationTableDefinitionColumn[]
    for col in table_def["DestinationTableDefinitionColumns"]
        tmp_col=NeuroData.DestinationTableDefinitionColumn(name=col["ColumnName"],
                datatype="Int",
                columntype="Value",
                isrequired=col["IsRequired"])

        tmp_col.ColumnDataTypeScale=col["ColumnDataTypeScale"]
        tmp_col.ColumnDataType=col["ColumnDataType"]
        tmp_col.ColumnDataTypePrecision=col["ColumnDataTypePrecision"]
        tmp_col.ColumnDataTypeSize=col["ColumnDataTypeSize"]
        tmp_col.ColumnType=col["ColumnType"]
        push!(cols,tmp_col)
    end
    indexes=DestinationTableDefinitionIndex[]
    for ind in table_def["DestinationTableDefinitionIndexes"]
        push!(indexes,DestinationTableDefinitionIndex(;indexname=ind["IndexName"],indexcolumnnames=[ind["IndexColumns"][i]["ColumnName"] for i = 1:length(ind["IndexColumns"])]))
    end

    datastoreid=table_def["DataStoreId"]
    schematypeid=table_def["SchemaType"]

    new_table_def=NeuroData.DestinationTableDefinition(allowdatachanges=table_def["AllowDataLossChanges"],
    columns=cols,name=table_def["DestinationTableName"],tableindexes=indexes,schematypeid=schematypeid)
    new_table_def.DataStoreId=datastoreid
    return new_table_def
end

type DataPopulationMappingSourceColumn
    DestinationColumnInfo::Dict{String,Any}
    DestinationColumnName::String
    IsMapped::Bool
    SourceColumnName::String
end

type DataPopulationMappingRequest
    DataPopulationMappingSourceColumns
    DestinationTableDefinitionId
    MappingName
    function DataPopulationMappingRequest(tableId,columns,mappingname)
        new(columns,tableId,mappingname)
    end
end

function create_table_mapping(;storename=nothing,tablename=nothing,mappingname=nothing,notmapped=Array{String,1}(),source_dest_name_pairs=Array{Tuple{String,String}}(0))
    if storename==nothing
        error("Supply data store name")
    end
    datastoreid=nothing
    try
        datastoreid=NeuroJulia.neurocall("80","datastoremanager","GetDataStores",Dict("StoreName"=>storename))["DataStores"][1]["DataStoreId"]
    catch
        error("Data Store name is not valid")
    end
    request=GetDestinationTableDefinitionRequest(tablename,datastoreid)
    table_def=NeuroJulia.neurocall("80","DataPopulation","GetDestinationTableDefinition",request)
    columns=NeuroData.DataPopulationMappingSourceColumn[]
    for col in table_def["DestinationTableDefinitions"][1]["DestinationTableDefinitionColumns"]
        if col["ColumnName"]!="NeuroverseLastModified"
            if findfirst(notmapped,col["ColumnName"])==0
                ismapped=true
                destcolumnname=col["ColumnName"]
                sourcecolumnname=destcolumnname
                if length(source_dest_name_pairs)>0
                   tmp_ind=findfirst(map(x->x[2],source_dest_name_pairs),destcolumnname)
                    if tmp_ind>0
                        sourcecolumnname=source_dest_name_pairs[findfirst(map(x->x[2],source_dest_name_pairs),destcolumnname)][1]
                    end
                end

                push!(columns,NeuroData.DataPopulationMappingSourceColumn(
                col,destcolumnname,ismapped,sourcecolumnname))
            end
        end
    end

    data=NeuroData.DataPopulationMappingRequest(table_def["DestinationTableDefinitions"][1]["DestinationTableDefinitionId"],columns,mappingname)
    NeuroJulia.neurocall("80","DataPopulation","CreateDataPopulationMapping",data)
end

"create_processed_table!(datastorename::String,tablename::String,columnnames::Array{String,1},columntypes::Array{String,1};partitionpath::Union{String,Void}=nothing)"
function create_processed_table!(datastorename::String,tablename::String,columnnames::Array{String,1},columntypes::Array{String,1};partitionpath::Union{String,Void}=nothing)
    schematype="Processed"

    columns=NeuroData.DestinationTableDefinitionColumn[]
    for col=1:length(columnnames)
        push!(columns,NeuroData.DestinationTableDefinitionColumn(;name=columnnames[col],datatype=columntypes[col],columntype="Value",isrequired=true))
    end

    table_def=NeuroData.DestinationTableDefinition(;allowdatachanges=false,columns=columns,
            name=tablename,schematype=schematype,partitionpath=partitionpath)

    NeuroData.create_destination_table(storename=datastorename,tabledefinition=table_def)
end
"delete_processed_table!(datastorename::String,tablename::String)"
function delete_processed_table!(datastorename::String,tablename::String)
    table_def=get_table_definition(storename=datastorename,tablename=tablename)
    if table_def.SchemaType!=3
        error("Table schema type is not processed")
    end
    NeuroJulia.neurocall("80","DataPopulation","DeleteDestinationTableDefinition",Dict("DestinationTableDefinitionId"=>table_def.DestinationTableDefinitionId))
end

max_stars_repo_name: SnowdenNeuroverse/NeuroJulia

max_issues_repo_name: SnowdenNeuroverse/NeuroJulia

max_forks_repo_name: SnowdenNeuroverse/NeuroJulia


text:
"""
    ASTCtx

AST code generation context.
"""
mutable struct ASTCtx
    ptr::Int
    help::Symbol
    version::Symbol
end

ASTCtx() = ASTCtx(1, gensym(:help), gensym(:version))

"""
    read_arg(ctx::ASTCtx)

Return the expression that read the next argument from ARGS.
"""
function read_arg(ctx::ASTCtx)
    ex = :(ARGS[$(ctx.ptr)])
    ctx.ptr += 1
    return ex
end

"""
    codegen(ctx, cmd)

Generate target code according to given context `ctx` from a command object `cmd`.
"""
function codegen end

"""
    codegen(cmd)

Generate Julia AST from given command object `cmd`. This will wrap
all the generated AST in a function `command_main`.
"""
function codegen(cmd::AbstractCommand)
    defs = Dict{Symbol,Any}()
    defs[:name] = :command_main
    defs[:args] = [Expr(:kw, :(ARGS::Vector{String}), :ARGS)]

    ctx = ASTCtx()
    defs[:body] = quote
        # let Julia throw InterruptException on SIGINT
        ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)
        $(codegen_scan_glob(ctx, cmd))
        $(codegen(ctx, cmd))
    end
    return combinedef(defs)
end

function codegen(ctx::ASTCtx, cmd::LeafCommand)
    return quote
        $(codegen_help(ctx, cmd))
        $(codegen_body(ctx, cmd))
    end
end

function codegen(ctx::ASTCtx, cmd::NodeCommand)
    return quote
        $(codegen_help(ctx, cmd))
        $(codegen_body(ctx, cmd))
    end
end

function codegen(ctx::ASTCtx, cmd::EntryCommand)
    quote
        $(codegen_help(ctx, cmd.root, xprint_help(cmd)))
        $(codegen_version(ctx, cmd.root, xprint_version(cmd)))
        $(codegen_body(ctx, cmd.root))
    end
end

function codegen_scan_glob(ctx::ASTCtx, cmd::EntryCommand)
    return codegen_scan_glob(ctx, cmd.root)
end

function codegen_scan_glob(ctx::ASTCtx, cmd::LeafCommand)
    # do not prompt help if no required args
    if nrequired_args(cmd.args) == 0
        return quote
            $(ctx.help) = -1
            $(_codegen_scan_glob(ctx, cmd))
        end
    else
        return quote
            $(ctx.help) = isempty(ARGS) ? 1 : -1
            $(_codegen_scan_glob(ctx, cmd))
        end
    end
end

function codegen_scan_glob(ctx::ASTCtx, cmd::NodeCommand)
    return quote
        $(ctx.help) = isempty(ARGS) ? 1 : -1
        $(_codegen_scan_glob(ctx, cmd))
    end
end

function _codegen_scan_glob(ctx::ASTCtx, cmd)
    quote
        $(ctx.version) = -1
        for i in 1:length(ARGS)
            x = ARGS[i]
            if x == "--help" || x == "-h"
                $(ctx.help) = i
                break
            end

            if x == "--version" || x == "-V"
                $(ctx.version) = i
                break
            end
        end
    end
end

function codegen_help(ctx::ASTCtx, cmd::NodeCommand, msg = xprint_help(cmd))
    return quote
        if $(ctx.help) == $(ctx.ptr)
            $msg
            return 0
        end
    end
end

function codegen_help(ctx::ASTCtx, cmd::LeafCommand, msg = xprint_help(cmd))
    return quote
        if $(ctx.help) >= $(ctx.ptr)
            $msg
            return 0
        end
    end
end

function codegen_version(ctx::ASTCtx, cmd::NodeCommand, msg)
    return quote
        if $(ctx.version) == $(ctx.ptr)
            $msg
            return 0
        end
    end
end

function codegen_version(ctx::ASTCtx, cmd::LeafCommand, msg)
    return quote
        if $(ctx.version) >= $(ctx.ptr)
            $msg
            return 0
        end
    end
end

function codegen_body(ctx::ASTCtx, cmd::NodeCommand)
    ex = Expr(:block)
    err_msg = "expect at least one argument for command $(cmd.name)"
    push!(ex.args, Expr(:if, :(length(ARGS) < $(ctx.ptr)), xerror(cmd, err_msg)))

    start = ctx.ptr
    for subcmd in cmd.subcmds
        push!(ex.args, Expr(:if, :($(read_arg(ctx)) == $(subcmd.name)), codegen(ctx, subcmd)))
        # reset ptr
        ctx.ptr = start
    end

    err_msg = :("Error: unknown command $(ARGS[$(ctx.ptr)])")
    push!(ex.args, xerror(cmd, err_msg))
    return ex
end

function codegen_body(ctx::ASTCtx, cmd::LeafCommand)
    parameters = gensym(:parameters)
    n_args = gensym(:n_args)
    nrequires = nrequired_args(cmd.args)
    ret = Expr(:block)
    validate_ex = Expr(:block)

    pushmaybe!(ret, codegen_params(ctx, parameters, cmd))

    if nrequires > 0
        err = xerror(
            cmd,
            :("command $($(cmd.name)) expect at least $($nrequires) arguments, got $($n_args)"),
        )
        push!(validate_ex.args, quote
            if $n_args < $nrequires
                $err
            end
        end)
    end

    # Error: too much arguments
    if isempty(cmd.args) || !last(cmd.args).vararg
        nmost = length(cmd.args)
        err = xerror(cmd, :("command $($(cmd.name)) expect at most $($nmost) arguments, got $($n_args)"))
        push!(validate_ex.args, quote
            if $n_args > $nmost
                $err
            end
        end)
    end

    push!(ret.args, :($n_args = length(ARGS) - $(ctx.ptr - 1)))
    push!(ret.args, validate_ex)
    push!(ret.args, codegen_call(ctx, parameters, n_args, cmd))
    push!(ret.args, :(return 0))
    return ret
end

function codegen_params(ctx::ASTCtx, params::Symbol, cmd::LeafCommand)
    hasparameters(cmd) || return

    regexes, actions = [], []
    arg = gensym(:arg)
    it = gensym(:index)

    for opt in cmd.options
        push!(regexes, regex_flag(opt))
        push!(regexes, regex_option(opt))

        push!(actions, read_forward(params, it, opt))
        push!(actions, read_match(params, it, opt))

        if opt.short
            push!(regexes, regex_short_flag(opt))
            push!(regexes, regex_short_option(opt))

            push!(actions, read_forward(params, it, opt))
            push!(actions, read_match(params, it, opt))
        end
    end

    for flag in cmd.flags
        push!(regexes, regex_flag(flag))
        push!(actions, read_flag(params, it, flag))

        if flag.short
            push!(regexes, regex_short_flag(flag))
            push!(actions, read_flag(params, it, flag))
        end
    end

    return quote
        $params = []
        $it = $(ctx.ptr)
        while !isempty(ARGS) && $(ctx.ptr) <= $it <= length(ARGS)
            $arg = ARGS[$it]
            if startswith($arg, "-") # is a flag/option
                $(xmatch(cmd, regexes, actions, arg))
            else
                $it += 1
            end
        end
    end
end

function codegen_call(ctx::ASTCtx, params::Symbol, n_args::Symbol, cmd::LeafCommand)
    ex_call = Expr(:call, cmd.entry)
    if hasparameters(cmd)
        push!(ex_call.args, Expr(:parameters, :($params...)))
    end

    # since we will check the position of arguments
    # the optional and variational arguments are always in the end
    for (i, arg) in enumerate(cmd.args)
        if arg.require
            push!(ex_call.args, xparse_args(arg, ctx.ptr + i - 1))
        else
            break
        end
    end

    all_required(cmd) && return ex_call

    # handle optional arguments
    ex = Expr(:block)
    if cmd.nrequire >= 0
        push!(ex.args, Expr(:if, :($n_args == $(cmd.nrequire)), ex_call))
    end

    expanded_call_ex = copy(ex_call)
    for i in cmd.nrequire+1:length(cmd.args)
        expanded_call_ex = copy(expanded_call_ex)
        push!(expanded_call_ex.args, xparse_args(cmd.args[i], ctx.ptr + i - 1))
        push!(ex.args, Expr(:if, :($n_args == $i), expanded_call_ex))
    end
    return ex
end

function read_forward(parameters, it, option::Option)
    sym = QuoteNode(cmd_sym(option))

    function action(m)
        arg = xparse_args(option.arg, :($it + 1))
        return quote
            $it < length(ARGS) || error("expect an argument")
            push!($parameters, $sym => $arg)
            deleteat!(ARGS, ($it, $it + 1))
            $it = $it - 1
        end
    end
end

function read_match(parameters, it, option::Option)
    type = option.arg.type
    sym = QuoteNode(cmd_sym(option))

    function action(m)
        arg = xparse(type, :(String($m[1])))
        return quote
            push!($parameters, $sym => $arg)
            deleteat!(ARGS, $it)
            $it = $it - 1
        end
    end
end

function read_flag(parameters, it, flag)
    sym = QuoteNode(cmd_sym(flag))
    function action(m)
        return quote
            push!($parameters, $sym => true)
            deleteat!(ARGS, $it)
            $it = $it - 1
        end
    end
end


text:
# Run the benchmarks for Julia (run.jl) and Python (run.py). Copy the results in the vectors below and plot

using VegaLite
using DataFrames

# Agents.jl benchmark results
jlresults = [
0.0623855,
0.1409112,
0.2434189,
0.3797918]#,
# 0.8206902,
# 0.7178797,
# 0.9326405,
# 1.364939,
# 1.5469985,
# 2.051536]

# Mesa benchmark results
pyresults = [
0.8042,
1.9173,
3.1619,
4.7495
]

size_range = ((100:100:1000) .* 100)[1:4]

dd = DataFrame(runtime=vcat(jlresults, pyresults), lang=vcat(["Agents.jl" for i in 1:4], ["Mesa" for i in 1:4]), nv=vcat(size_range, size_range));

@vlplot(data=dd,
  mark={:line,
    size=4,
    point={filled=false, fill="white", size=60}
  },
  x={:nv,
    type="ordinal",
    title="Grid size",
    axis={titleFontSize=16, labelFontSize=14},
    grid=true,
  },
  y={:runtime,
    type="quantitative",
    title="Run time (seconds)",
    axis={titleFontSize=13,
    labelFontSize=14},
    grid=true
  },
  color={:lang,
    type="nominal",
    title="Framework",
    legend={titleFontSize=15, labelFontSize=14},
    scale={range=["#EA98D2", "#659CCA"]},
  },
  height=400,
  width=500,
)

dd2 = DataFrame(ratio=pyresults./jlresults, nv=size_range);

@vlplot(data=dd2,
  mark={:bar},
  x={:nv,
    type="ordinal",
    title="Grid size",
    axis={titleFontSize=16, labelFontSize=14},
  },
  y={:ratio,
    type="quantitative",
    title="Mesa/Agents run time",
    axis={titleFontSize=13,
    labelFontSize=14}
  },
  height=400,
  width=500,
)


text:
abstract type Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion} <: AbstractFermionfields_5D{NC}  end

"""
Struct for MobiusDomainwallFermion
"""
struct MobiusDomainwallFermion_5D{NC,WilsonFermion} <: Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion}
    w::Array{WilsonFermion,1}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    L5::Int64
    Dirac_operator::String
    NWilson::Int64
    nowing::Bool

    function MobiusDomainwallFermion_5D(L5,NC::T,NX::T,NY::T,NZ::T,NT::T;nowing = false) where T<: Integer
        if nowing
            x = WilsonFermion_4D_nowing(NC,NX,NY,NZ,NT)
        else
            x = WilsonFermion_4D_wing(NC,NX,NY,NZ,NT)
        end
        xtype = typeof(x)
        w = Array{xtype,1}(undef,L5)
        w[1] = x
        for i=2:L5
            w[i] = similar(x)
        end
        #println(w[2][1,1,1,1,1,1])
        NWilson = length(x)
        Dirac_operator = "MobiusDomainwall"
        return new{NC,xtype}(w,NC,NX,NY,NZ,NT,L5,Dirac_operator,NWilson,nowing)
    end

end

function Base.similar(x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion} ) where {NC,WilsonFermion}
    return MobiusDomainwallFermion_5D(x.L5,NC,x.NX,x.NY,x.NZ,x.NT,nowing=x.nowing)
end

function apply_J!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion}) where  {NC,WilsonFermion}
    clear_fermion!(xout)

    L5 = xout.L5
    for i5=1:L5
        j5 = L5-i5+1
        substitute_fermion!(xout.w[i5],x.w[j5])
    end
end

function apply_P!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion}) where  {NC,WilsonFermion}
    L5 = xout.L5
    clear_fermion!(xout)

    for i5=1:L5
        j5 = i5
        #P_- -> P_+ in this definition
        mul_1plusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])

        #P_+ -> P_- in this definition
        j5 = i5+1
        j5 += ifelse(j5 > L5,-L5,0)
        mul_1minusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])
    end
end

function apply_Pdag!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion}) where  {NC,WilsonFermion}
    L5 = xout.L5
    clear_fermion!(xout)

    for i5=1:L5
        j5 = i5
        #P_- -> P_+ in this definition
        mul_1plusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])

        #P_+ -> P_- in this definition
        j5 = i5-1
        j5 += ifelse(j5 < 1,L5,0)
        mul_1minusγ5x_add!(xout.w[i5],x.w[j5],1)
        set_wing_fermion!(xout.w[i5])
    end
end

function apply_1pD!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},L5,U::Array{G,1},A,
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},factor) where  {NC,WilsonFermion,G <: AbstractGaugefields}
    clear_fermion!(xout)

    if factor == 0
        for i5=1:xout.L5
            axpy!(1,x.w[i5],xout.w[i5])
        end
        return
    end

    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end
    ratio = 1

    for i5 in irange
        j5=i5
        D4x!(xout.w[i5],U,x.w[j5],A,4) #Dw*x
        #Dx!(xout.w[i5],U,x.w[j5],A) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #1/(2*A.κ)
        massfactor = -(1/(2*A.κ) + 1)
        set_wing_fermion!(xout.w[i5])
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        add!(ratio,xout.w[i5],factor*ratio*massfactor,x.w[j5]) #D = x + Dw*x
        set_wing_fermion!(xout.w[i5])
    end

    if L5 != xout.L5
        for i5 in irange_out
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end

end

function apply_1pDdag!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},L5,U::Array{G,1},A,
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},factor) where  {NC,WilsonFermion,G <: AbstractGaugefields}
    clear_fermion!(xout)

    if factor == 0
        for i5=1:xout.L5
            axpy!(1,x.w[i5],xout.w[i5])
        end
        return
    end

    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end
    ratio = 1

    for i5 in irange
        j5=i5
        D4dagx!(xout.w[i5],U,x.w[j5],A,4) #Dw*x
        #Dx!(xout.w[i5],U,x.w[j5],A) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #1/(2*A.κ)
        massfactor = -(1/(2*A.κ) + 1)
        set_wing_fermion!(xout.w[i5])
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        add!(ratio,xout.w[i5],factor*ratio*massfactor,x.w[j5]) #D = x + Dw*x
        set_wing_fermion!(xout.w[i5])
    end

    if L5 != xout.L5
        for i5 in irange_out
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end

end

function apply_F!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},L5,m,
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},temp1) where  {NC,WilsonFermion}
    clear_fermion!(xout)

    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end
    ratio = 1

    for i5 in irange
        j5=i5+1
        if 1 <= j5 <= xout.L5
            #-P_- -> - P_+ :gamma_5 of LTK definition
            if xout.L5 != 2
                # ((1 + (b-c)/2 DW )*P_-*x)
                mul_1plusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end
        end

        j5=i5-1
        if 1 <= j5 <= xout.L5
            #-P_+ -> - P_- :gamma_5 of LTK definition
            if xout.L5 != 2
                mul_1minusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end
        end

        if xout.L5 != 1
            if i5==1
                j5 = xout.L5
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1minusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],-m*ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end

            if i5== xout.L5
                j5 = 1
                mul_1plusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],-m*ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end
        end
    end

    #if L5 != xout.L5
    #    for i5 in irange_out
    #        axpy!(1,x.w[i5],xout.w[i5])
    #    end
    #end

end

function apply_Fdag!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},L5,m,
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},temp1) where  {NC,WilsonFermion}
    clear_fermion!(xout)

    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end
    ratio = 1

    for i5 in irange
        j5=i5+1
        if 1 <= j5 <= xout.L5
            #-P_- -> - P_+ :gamma_5 of LTK definition
            if xout.L5 != 2
                # ((1 + (b-c)/2 DW )*P_-*x)
                mul_1minusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end
        end

        j5=i5-1
        if 1 <= j5 <= xout.L5
            #-P_+ -> - P_- :gamma_5 of LTK definition
            if xout.L5 != 2
                mul_1plusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end
        end

        if xout.L5 != 1
            if i5==1
                j5 = xout.L5
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1plusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],-m*ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end

            if i5== xout.L5
                j5 = 1
                mul_1minusγ5x!(temp1.w[i5],x.w[j5])
                add!(ratio,xout.w[i5],-m*ratio,temp1.w[i5])
                set_wing_fermion!(xout.w[i5])
            end
        end
    end

    #if L5 != xout.L5
    #    for i5 in irange_out
    #        axpy!(1,x.w[i5],xout.w[i5])
    #    end
    #end

end

function D5DWx!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},U::Array{G,1},
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},m,A,L5,b,c,
        temp1,temp2) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    coeff_plus = (b+c)/2
    #coeff_plus = 1
    coeff_minus = -(b-c)/2
    #coeff_minus  = 0
    clear_fermion!(xout)
    ratio = 1

    factor = coeff_plus
    #xout = (1 + factor*D)*x
    apply_1pD!(xout,L5,U,A,x,factor)

    #temp2 = F*x
    apply_F!(temp2,L5,m,x,temp1)
    factor = coeff_minus
    #xout = (1 + factor*D)*F*x
    apply_1pD!(temp1,L5,U,A,temp2,factor)
    for i5=1:L5
        axpy!(1,temp1.w[i5],xout.w[i5])
    end

    set_wing_fermion!(xout)

    return
end

function D5DWdagx!(xout::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion} ,U::Array{G,1},
    x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion} ,m,A,L5,b,c,
    temp1,temp2) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)
    ratio = 1
    coeff_plus = (b+c)/2
    coeff_minus = -(b-c)/2
    #coeff_plus = 1
    #coeff_minus  = 0
    #ratio = xout.L5/L5

    factor = coeff_minus
    #xout = Fdag*(1 + factor*Ddag)*x
    apply_1pDdag!(temp2,L5,U,A,x,factor)
    apply_Fdag!(xout,L5,m,temp2,temp1)

    factor = coeff_plus
    #xout = (1 + factor*Ddag)*x
    apply_1pDdag!(temp1,L5,U,A,x,factor)

    for i5=1:L5
        axpy!(1,temp1.w[i5],xout.w[i5])
    end

    set_wing_fermion!(xout)

    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion}) where {NC,WilsonFermion}
    L5 = length(x.w)
    for iL=1:L5
        gauss_distribution_fermion!(x.w[iL])
    end
    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::Abstract_MobiusDomainwallFermion_5D{NC,WilsonFermion},randomfunc,σ) where {NC,WilsonFermion}
    L5 = length(x.w)
    for iL=1:L5
        gauss_distribution_fermion!(x.w[iL],randomfunc,σ)
    end
    return

end


text:
using Plots
using Plots.PlotMeasures
using Colors
using StatsPlots

include("collect_data.jl")

function plot_A()

    fcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)

    fermi = 1000*get_fermi()
    psi4  = 1000*get_psi4()

    groupedbar([fermi psi4], c = [fcolor pcolor], legend=:topright,
                label=["Fermi" "Psi4"], dpi=1200, orientation = :v,
                bar_width=0.8,xticks=[2*i for i = 1:11], framestyle=:box)
    ylabel!("Time (milliseconds)")
    xlabel!("S22 Dimer")

    savefig("A.png")
end

function plot_B()

    fcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)

    fermi = 1000*get_fermi()
    psi4  = 1000*get_psi4()

    nao = [58, 48, 104, 114, 264, 251, 321, 68, 96, 148, 228, 208, 264, 275, 321, 86, 138, 143, 147, 228, 275, 256]
    perm = sortperm(nao)
    nao = nao[perm]
    fermi = fermi[perm]
    psi4 = psi4[perm]
    names = String[]

    for i = 1:22
        N = nao[i]
        S = perm[i]
        push!(names, "$S")
    end

    groupedbar(names, [fermi psi4], c = [fcolor pcolor], legend=:topleft,
                label=["Fermi" "Psi4"], dpi=1200, orientation = :v,
                bar_width=0.8,framestyle=:box,
                right_margin=1cm)

    ylabel!("Time (milliseconds)")
    xlabel!("S22 Dimer")
end

function plot_C()

    fcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)

    fermi = 1000*get_fermi()
    psi4  = 1000*get_psi4()

    nao = [58, 48, 104, 114, 264, 251, 321, 68, 96, 148, 228, 208, 264, 275, 321, 86, 138, 143, 147, 228, 275, 256]
    perm = sortperm(nao)
    nao = nao[perm]
    fermi = fermi[perm]
    psi4 = psi4[perm]
    names = String[]
    snbf = String[]

    for i = 1:22
        N = nao[i]
        S = perm[i]
        push!(names, "S$S")
        push!(snbf, "$N")
    end

    p1 = groupedbar(names, [fermi psi4], c = [fcolor pcolor], legend=:topleft,
                label=["Fermi" "Psi4"], dpi=1200, orientation = :v,
                bar_width=0.8,framestyle=:box,
                right_margin=1cm)

    ylabel!("Time (milliseconds)")
    xlabel!("S22 Dimer")
    p2 = scatter(nao, fermi ./ psi4, label="", yticks=[0.5, 1.0, 1.5], ylims=[0.3,1.8],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:utriangle, markercolor=:teal,
                xlabel="Number of basis functions")
    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]))
    #p1
end

function plot_D()

    fcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)

    fermi = 1000*get_fermi()
    psi4  = 1000*get_psi4()

    nao = [58, 48, 104, 114, 264, 251, 321, 68, 96, 148, 228, 208, 264, 275, 321, 86, 138, 143, 147, 228, 275, 256]
    perm = sortperm(nao)
    nao = nao[perm]

    p1 = groupedbar([fermi psi4], c = [fcolor pcolor], legend=:topright,
                label=["Fermi" "Psi4"], dpi=1200, orientation = :v,
                bar_width=0.8,xticks=([2, 6, 12, 18, 22],["S2","S6", "S12","S18","S22"]), framestyle=:box)

    ylabel!("Time (milliseconds)")
    xlabel!("Dimer")

    p2 = scatter(nao, fermi[perm] ./ psi4[perm], label="", yticks=[0.5, 1.0, 1.5], ylims=[0.0,2.0],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:diamond, markercolor=:navajowhite3,
                xlabel="Number of basis functions")
    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]), leftmargin=0.8cm)
end

function plot_E()

    fcolor = ARGB(186/255,12/255,47/255, 0.8)
    pcolor = ARGB(39/255, 56/255, 150/255, 0.75)

    fermi = 1000 .* get_fermi()
    psi4  = 1000 .* get_psi4()

    nao = [58, 48, 104, 114, 264, 251, 321, 68, 96, 148, 228, 208, 264, 275, 321, 86, 138, 143, 147, 228, 275, 256]
    perm = sortperm(nao)
    nao = nao[perm]

    p1 = groupedbar(([fermi psi4]), c = [fcolor pcolor], legend=:topleft,
                label=["Fermi" "Psi4"], dpi=1200, orientation = :v, yaxis=:log,
                ylim = (10^-0.5,10^4),
                bar_width=0.7,xticks=([2, 6, 12, 18, 22],["S2","S6", "S12","S18","S22"]), framestyle=:box)

    ylabel!("Time (milliseconds)")
    xlabel!("Dimer")

    p2 = scatter(nao, fermi[perm] ./ psi4[perm], label="", yticks=[0.5, 1.0, 1.5], ylims=[0.0,2.0],
                color=:black, ylabel="Relative time\n(Fermi/Psi4)", framestyle=:box, marker=:diamond, markercolor=fcolor,
                xlabel="Number of basis functions")

    hline!(p2, [1], color=:gray, linestyle=:dot, label="")

    plot(p1, p2, layout = grid(2,1, heights=[0.7, 0.3]), leftmargin=0.8cm, dpi=1200)

    savefig("build_fock.png")
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
abstract AbstractCalibration

type CalibrationSpline1D <: AbstractCalibration
	spline::Spline1D
	known_energies::Vector{Float64}
	indicator::Vector{Float64}
	iscalibrated::Bool
end
function CalibrationSpline1D(indicator, known_energies)
	CalibrationSpline1D(calc_spline(indicator, known_energies), known_energies, indicator, true)
end
iscalibrated(c::AbstractCalibration) = c.iscalibrated
donethru(c::AbstractCalibration) = iscalibrated(c)?DONETHRU_MAX:0
apply_calibration(c::CalibrationSpline1D, estimator) = apply_spline(c.spline, estimator)
function PyPlot.plot(cal::CalibrationSpline1D)
	figure()
	plot(cal.indicator, apply_spline(cal.spline, cal.indicator))
	plot(cal.indicator, cal.known_energies,"o")
	xlabel("indicator")
	ylabel("energy")
end

function plotlog(cal::CalibrationSpline1D)
	figure()
	x=log(cal.indicator)
	y=log(cal.known_energies)
	# plot(cal.indicator, apply_spline(cal.spline, cal.indicator))
	plot(x, y,"o")
	p = polyfit(x,y,2)
	plot(x, polyval(p, x))
	xlabel("indicator")
	ylabel("energy")
	title("peakassign_fig_of_merit $(peakassign_fig_of_merit(x,y))")
end

function peakassign(locations_arb, locations_true)
# locations_arb is an array of possible peak locations in arbitrary units, assumed to have an unknown,
# nonlinear (but monotonic and smooth) relationship to locations_true
# locations_arb may contain more peaks than locations_true, including some spurious peaks
# the goal is to find the best possible set of assignments for entries in locations_arb to
# correspond to locations_true by examining all possible assignments and calcualting some figure of merit
@assert(issorted(locations_arb))
@assert(issorted(locations_true))
@assert(length(locations_arb)>=length(locations_true))
combos = combinations(1:length(locations_arb), length(locations_true))
length(combos)==1 && return (locations_arb[first(combos)], 0.0)

assert(!(0 in locations_true)) # don't pass zero and force zero
assert(!(0 in locations_arb))

log_locations_arb = log(locations_arb)
log_locations_true = log(locations_true)

best_combo = first(combos)
best_fom = realmax(Float64)
for c in combos
	cu = log_locations_arb[c]
	fom = peakassign_fig_of_merit(cu, log_locations_true)
	if fom < best_fom
		best_fom = fom
		best_combo = c
	end
end
locations_arb[best_combo], best_fom
end
function peakassign_fig_of_merit(x, y)
	n = length(x) > 2 ? 2 : 1
	p = polyfit(x,y,n)
	std(polyresiduals(p,x,y)./y)
end
# should remove these and depend on Polynomails.jl or CurveFit or something
function polyfit(x, y, n)
  A = [ float(x[i])^p for i = 1:length(x), p = 0:n ]
  A \ y
end
function polyval(p,x::Number)
	y = zero(promote_type(eltype(p), eltype(x)))
	for i=1:length(p)
		y+=p[i]*x^(i-1)
	end
	y
end
function polyval(p,x::AbstractArray)
	y = zeros(promote_type(eltype(p), eltype(x)), length(x))
	for i = 1:length(x)
		y[i] = polyval(p,x[i])
	end
	y
end
function polyresiduals(p,x,y)
	yp = polyval(p,x)
	residuals = y-yp
end

findpeaks(h::Histogram; fwhm=15,n=typemax(Int)) = findpeaks(bin_edges(h), counts(h), fwhm=fwhm,n=n)
function findpeaks(edges::Range, y::Vector; fwhm=15,n=typemax(Int))
	k = UnivariateKDE(edges, convert(Array{Float64,1},y))
	findpeaks_(k,fwhm,n)
end

function findpeaks_(y::UnivariateKDE, fwhm,n)
	# smooth y with a gaussian of FHWM=fwhm, return indicies of n peaks with highest density in smoothed version
	σ = fwhm/sqrt(8*log(2))
	d = Normal(0,σ)
	smoothed = conv(y,d)
	peakinds = findpeaks_(smoothed.density)
	if length(peakinds)>n
		peakinds = peakinds[end-n+1:end]
	end
	peakinds
end

function findpeaks_{T}(y::Vector{T})
	# find all local maxima, return a vector of indicies sorted by value of y at that indey
	peakinds = Array(Int,0)
	for i = 2:length(y)-1
		if y[i-1]<y[i]>y[i+1]
			push!(peakinds, i)
		end
	end
	p = sortperm(y[peakinds])
	peakinds[p]
end

plotpeaks(h::Histogram; fwhm=15, n=typemax(Int)) = plotpeaks(bin_edges(h), counts(h), fwhm=fwhm, n=n)
function plotpeaks(edges::Range, y::Vector; fwhm=15, n=typemax(Int))
	k = UnivariateKDE(edges, convert(Array{Float64,1},y))
	peakinds = findpeaks_(k, fwhm, n)
	peaky = midpoints(edges)[peakinds]
	figure()
	plot(midpoints(edges), y, label="y")
	plot(midpoints(k.x), conv(k,Normal(0,fwhm/2.355)).density, label="smoothed")
	plot(peaky, y[peakinds],"o",label="raw peaks")
	plot(peaky, conv(k,Normal(0,fwhm/2.355)).density[peakinds],"o",label="smoothed peaks")
	legend()
end

function calibrate_nofit(hist::Histogram, known_energies, n_extra=3)
	peak_inds = findpeaks(hist,fwhm=15, n=length(known_energies)+n_extra) # find the largest peaks, at most 3 more than the number of known energies
	# with a 15 unit gaussian applied for smoothing... this needs to be adjustable
	# also we need to check if calibration works, and if not return a failed calibration object with debug info
	peak_x = sort(bin_centers(hist)[peak_inds])
	best_combo, best_fom = peakassign(sort(peak_x), known_energies)
	CalibrationSpline1D(best_combo, known_energies)
end

function calibrate_xy(indicator, known_energies)
	CalibrationSpline1D(indicator, known_energies)
end

function calc_spline(x,y)
	spl = Spline1D(float(x),float(y),k=1,bc="extrapolate")
end
function apply_spline(spl, energy_estimator)
	energy = evaluate(spl, energy_estimator)
end


text:
"""
    Car

    Unicycle model (http://lavalle.pl/planning/)
"""

struct Car{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int
end

function f(model::Car, x, u, w)
    @SVector [u[1] * cos(x[3]),
              u[1] * sin(x[3]),
              u[2]]
end

n, m, d = 3, 2, 3
model = Car{Midpoint, FixedTime}(n, m, d)

function visualize!(vis, model::Car, x;
       Δt = 0.1, traj = false)

	default_background!(vis)

	obj_path = joinpath(pwd(),"models/cybertruck/cybertruck.obj")
	mtl_path = joinpath(pwd(),"models/cybertruck/cybertruck.mtl")

	m = ModifiedMeshFileObject(obj_path,mtl_path,scale=0.05)
	setobject!(vis["car"],m)

   	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    for t = 1:length(x)
        MeshCat.atframe(anim,t) do
			settransform!(vis["car"],
				compose(Translation(x[t][1], x[t][2], 0.0),
					LinearMap(RotZ(x[t][3] + pi) * RotX(pi / 2.0))))
        end
    end

	if traj
		pts_nom = collect(eachcol(hcat(x...)))
		material_nom = LineBasicMaterial(color = colorant"orange", linewidth = 3.0)
		setobject!(vis["traj"], Object(PointCloud(pts_nom), material_nom, "Line"))
	end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, 0.0, 0.0),LinearMap(RotY(-pi/2.5))))

    MeshCat.setanimation!(vis, anim)
end


text:
# this script generates Julia artifacts for debiasing tables from
# Chesley et al. (2010), Farnocchia et al. (2015) and Eggl et al. (2020)
# TODO: adapt script for julia1.6+
# TODO: adapt script for 2010 debiasing tables (requires extra unpacking)

using Pkg.Artifacts
using Pkg.PlatformEngines
using Pkg.PlatformEngines: sha256
using Pkg
# using Downloads: download
Pkg.PlatformEngines.probe_platform_engines!()

# This is the path to the Artifacts.toml we will manipulate
artifact_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")

ftp_jpl = "ftp://ssd.jpl.nasa.gov/pub/ssd/debias/"

names_debias = ["debias", "debias_2014", "debias_2018", "debias_hires2018"]
urls_debias = ftp_jpl .* names_debias .* ".tgz"

names_lsk_spk = ["naif0012", "de430", "TTmTDBde430", "a99942", "ttmtdb_DE430_1995_2030"]
urls_lsk_spk = [
    "https://raw.githubusercontent.com/PerezHz/jpleph/main/naif0012.tar.gz",
    "https://raw.githubusercontent.com/PerezHz/jpleph/main/de430.tar.gz",
    "https://raw.githubusercontent.com/PerezHz/jpleph/main/TTmTDBde430.tar.gz",
    "https://raw.githubusercontent.com/PerezHz/jpleph/main/a99942.tar.gz",
    "https://raw.githubusercontent.com/PerezHz/jpleph/main/ttmtdb_DE430_1995_2030.tar.gz"
]

urls_ = vcat(urls_lsk_spk, urls_debias)
names_ = vcat(names_lsk_spk, names_debias)

for (url, name) in map((x,y)->(x,y), urls_, names_)
    # Query the `Artifacts.toml` file for the hash bound to the name
    # (returns `nothing` if no such binding exists)
    artfct_hash = artifact_hash(name, artifact_toml)

    # If the name was not bound, or the hash it was bound to does not exist, create it!
    if artfct_hash == nothing || !artifact_exists(artfct_hash)
        tb = Tuple[]
        # create_artifact() returns the content-hash of the artifact directory once we're finished creating it
        artfct_hash = create_artifact() do artifact_dir
            # We create the artifact by simply downloading a few files into the new artifact directory
            @show url
            tmp_dir = mktempdir()
            tarball = download(url, joinpath(tmp_dir, basename(url)))
            # @show tarball tmp_dir
            try
                tarball_hash = open(tarball) do file
                    bytes2hex(sha256(file))
                end
                @show tarball_hash
                tb = [(url, tarball_hash)]
                unpack(tarball, artifact_dir)
            finally
                rm(tarball)
            end
        end
        # Now bind that hash within our `Artifacts.toml`. `force = true` means that if it already exists,
        # just overwrite the corresponding entry in `Artifacts.toml` with the new content-hash.
        # Unless the source files change, we do not expect the content hash to
        # change, so this should not cause unnecessary version control churn.
        bind_artifact!(artifact_toml, name, artfct_hash, lazy=true, force=true, download_info=tb)
    end

    # Get the path of the artifact, either newly created or previously generated.
    # this should be something like `~/.julia/artifacts/dbd04e28be047a54fbe9bf67e934be5b5e0d357a`
    artfct_path = artifact_path(artfct_hash)
    @show artfct_path
end