text:
export short_vectors, shortest_vectors

################################################################################
#
#  Auto- and isomorphism computation of lattices
#
################################################################################

# This is a port of the program ISOM and AUTO by Bernd Souvignier
# which implemented an algorithm published in
# W. PLESKEN, B. SOUVIGNIER, Computing Isometries of Lattices,
# Journal of Symbolic Computation, Volume 24, Issues 3-4, September 1997,
# Pages 327-334, ISSN 0747-7171, 10.1006/jsco.1996.0130.
# (https://www.sciencedirect.com/science/article/pii/S0747717196901303)

mutable struct SCPComb
	rank::Int
	trans::fmpz_mat
	coef::fmpz_mat
  F::Vector{fmpz_mat}

  SCPComb() = new()
end

mutable struct VectorList{S, T}
  vectors::Vector{S}
  lengths::Vector{Vector{T}}
  lookup::Dict{S, Int}
  issorted::Bool
  use_dict::Bool

  function VectorList{S, T}() where {S, T}
    return new{S, T}()
  end
end

function VectorList(vectors::Vector{S}, lengths::Vector{Vector{T}},
                    use_dict::Bool = true) where {S, T}

  V = VectorList{S, T}()
  if use_dict
    V.lookup = Dict{S, Int}(vectors[i] => i for i in 1:length(vectors))
    V.lengths = lengths
    V.vectors = vectors
    V.use_dict = true
  else
    p = sortperm(vectors)
    permute!(vectors, p)       # apply the permutation to V
    permute!(lengths, p)       # apply the permutation to lengths
    V.use_dict = false
    V.vectors = vectors
    V.lengths = lengths
  end

  return V
end

length(V::VectorList) = length(V.vectors)

Base.issorted(V::VectorList) = V.issorted

getindex(V::VectorList, i::Int) = i > 0 ? V.vectors[i] : -V.vectors[-i]

function isnormalized(w::fmpz_mat)
  for k in 1:ncols(w)
    if !iszero(w[1, k])
      return w[1, k] > 0
    end
  end
end

function isnormalized(w::Vector{Int})
  for k in 1:length(w)
    if !iszero(w[k])
      return w[k] > 0
    end
  end
end

function neg!(w::Vector{Int})
  w .*= -1
end

function neg!(w::fmpz_mat)
  ccall((:fmpz_mat_neg, libflint), Cvoid, (Ref{fmpz_mat}, Ref{fmpz_mat}), w, w)
  return w
end

function find_point(w, V::VectorList)
  fl, k = has_point(w, V)
  @assert fl
  return k
end
#  positive = isnormalized(w)
#
#  if positive
#    k = V.lookup[w]
#    #k = findfirst(isequal(w), V.vectors)
#    #@assert k !== nothing
#    #@assert V[k] == w
#    #return k
#    return k
#  end
#
#  neg!(w)
#  #k = findfirst(isequal(w), V.vectors)
#  #@assert k !== nothing
#  k = V.lookup[w]
#  neg!(w)
#  @assert V[-k] == w
#  return -k
#end

function has_point(w, V::VectorList)
  positive = isnormalized(w)

  if positive
    k = get(V.lookup, w, 0)
    if k == 0
      return false, 0
    else
      @assert V[k] == w
      return true, k
    end
  end

  neg!(w)
  #k = findfirst(isequal(w), V.vectors)
  #@assert k !== nothing
  k = get(V.lookup, w, 0)
  neg!(w)
  if k == 0
    return false, 0
  else
    @assert V[-k] == w
    return true, -k
  end
end

function _find_point(w::fmpz_mat, V::VectorList{fmpz_mat, T}) where T
  positive = false
  for k in 1:length(w)
    if !iszero(w[1, k])
      positive = w[1, k] > 0
      break
    end
  end
  if positive
    if issorted(V)
      k = searchsortedfirst(V.vectors, w)
    else
      k = findfirst(isequal(w), V.vectors)
    end
    @assert k !== nothing
    return k
  else
    mw = -w
    if issorted(V)
      k = searchsortedfirst(V.vectos, mw)
    else
      k = findfirst(isequal(mw), V.vectors)
    end
    @assert k !== nothing
    return -k
  end
end

mutable struct ZLatAutoCtx{S, T, V}
  G::Vector{T}
  Gtr::Vector{T}
  dim::Int
  max::S
  V::VectorList{V, S}
  v::Vector{T}
  per::Vector{Int}
  fp::Matrix{Int}
  fp_diagonal::Vector{Int}
  std_basis::Vector{Int}
  scpcomb::SCPComb

  orders::Vector{Int}
  ng::Vector{Int}
  nsg::Vector{Int}
  g::Vector{Vector{T}}
  prime::S

  issymmetric::BitArray{1}
  operate_tmp::V

  function ZLatAutoCtx(G::Vector{fmpz_mat})
    z = new{fmpz, fmpz_mat, fmpz_mat}()
    z.G = G
    z.Gtr = fmpz_mat[transpose(g) for g in G]
    z.dim = nrows(G[1])
    z.issymmetric = falses(length(G))
    z.operate_tmp = zero_matrix(FlintZZ, 1, ncols(G[1]))

    for i in 1:length(z.G)
      z.issymmetric[i] = issymmetric(z.G[i])
    end

    return z
  end

  function ZLatAutoCtx{S, T, V}() where {S, T, V}
    return new{S, T, V}()
  end
end

function Base.show(io::IO, C::ZLatAutoCtx)
  print(io, "Automorphism context for ", C.G)
end

dim(C::ZLatAutoCtx) = C.dim

function LinearAlgebra.issymmetric(M::MatElem)
  for i in 1:nrows(M)
    for j in i:ncols(M)
      if M[i, j] != M[j, i]
        return false
      end
    end
  end
  return true
end

function init(C::ZLatAutoCtx, auto::Bool = true, bound::fmpz = fmpz(-1), use_dict::Bool = true)
  # Compute the necessary short vectors

  r = length(C.G)

  n = nrows(C.G[1])

  if bound == -1
    bound = maximum(diagonal(C.G[1]))
    C.max = bound
  end

  @assert bound > 0

  @vprint :Lattice 1 "Computing short vectors of length $bound\n"

  @vtime :Lattice 1 V = _short_vectors_gram_integral(C.G[1], bound)

  vectors = Vector{fmpz_mat}(undef, length(V))

  lengths = Vector{Vector{fmpz}}(undef, length(V))

  tmp = zero_matrix(FlintZZ, 1, n)

  for i in 1:length(V)
    # First canonicalize them
    cand = V[i]
    v = cand[1]
    k = 1
    while iszero(v[k])
      k += 1
    end
    if v[k] < 0
      v .*= -1
    end

    vfmpz = matrix(FlintZZ, 1, n, v)

    w = Vector{fmpz}(undef, r)
    w[1] = cand[2]
    for k in 2:r
      w[2] = _norm(vfmpz, C.G[k], tmp)
    end

    lengths[i] = w
    vectors[i] = vfmpz
  end

  V = VectorList(vectors, lengths, use_dict)

  for i in 1:length(C.G)
    C.issymmetric[i] = issymmetric(C.G[i])
  end

  @assert C.issymmetric[1]

  C.V = V

  # Compute the fingerprint
  @vprint :Lattice 1 "Computing fingerprint: "
  if auto
    @vtime :Lattice 1 fingerprint(C)
    @vprint :Lattice 1 "$(C.fp_diagonal)\n"
  end

  if auto
    # Find the standard basis vectors
    C.std_basis = Vector{Int}(undef, dim(C))
    z = zero_matrix(FlintZZ, 1, dim(C))
    for i in 1:dim(C)
      z[1, C.per[i]] = 1
      k = find_point(z, C.V)
      C.std_basis[i] = k
      z[1, C.per[i]] = 0
    end
  end

  #

  C.v = Vector{fmpz_mat}(undef, length(C.G))

  for i in 1:length(C.G)
    A = zero_matrix(FlintZZ, length(C.V), dim(C))
    for j in 1:length(C.V)
      for k in 1:dim(C)
        A[j, k] = _dot_product_with_row(C.V[j], C.G[i], k)
      end
    end
    C.v[i] = A
  end

  if false
    for i in 1:length(C.G)
      for j in 1:length(C.V)
        for k in 1:length(C.V)
          @assert _dot_product(C.V[j], C.v[i], k) == (C.V[j] * C.G[i] * transpose(C.V[k]))[1, 1]
        end
      end
    end
  end

  C.g = Vector{Vector{fmpz_mat}}(undef, dim(C))
  for i in 1:dim(C)
    C.g[i] = fmpz_mat[]
  end
  C.ng = zeros(Int, dim(C))
  C.nsg = zeros(Int, dim(C))
  C.orders = Vector{Int}(undef, dim(C))

  # -Id is always an automorphism
  C.g[1] = fmpz_mat[-identity_matrix(FlintZZ, dim(C))]
  C.ng[1] = 1

  # Calculate orbit lengths

  nH = 0

  for i in 1:dim(C)
    nH += C.ng[i]
  end

  H = Vector{fmpz_mat}(undef, nH)

  if auto
    for i in 1:dim(C)
      if C.ng[i] > 0
        nH = 0
        for j in i:dim(C)
          for k in 1:C.ng[j]
            nH += 1
            H[nH] = C.g[j][k]
          end
        end
        #@assert _orbitlen_naive(C.std_basis[i], C.fp_diagonal[i], H, nH, C.V) == _orbitlen(C.std_basis[i], C.fp_diagonal[i], H, nH, C.V)
        C.orders[i] = _orbitlen(C.std_basis[i], C.fp_diagonal[i], H, nH, C.V, C)
      else
        C.orders[i] = 1
      end
    end
  end

  return C
end

# The following functions tries to initialize a ZLatAutoCtx with entries in Int.
# The return value is flag, Csmall
function try_init_small(C::ZLatAutoCtx, auto::Bool = true, bound::fmpz = fmpz(-1), use_dict::Bool = true)
  # Compute the necessary short vectors
  @vprint :Lattice 1 "Computing short vectors of length $max\n"
  automorphism_mode = bound == fmpz(-1)

  Csmall = ZLatAutoCtx{Int, Matrix{Int}, Vector{Int}}()

  if bound == -1
    bound = maximum(diagonal(C.G[1]))
    Csmall.max = Int(bound)
  else
    Csmall.max = -1
  end
  @assert bound > 0

  @vtime :Lattice 1 V = _short_vectors_gram_integral(C.G[1], bound)

  vectors = Vector{Vector{Int}}(undef, length(V))

  lengths = Vector{Vector{Int}}(undef, length(V))

  r = length(C.G)

  Gsmall = Matrix{Int}[Matrix{Int}(g) for g in C.G]

  Gsmall_nbits = 0
  for k in 1:r
    Gsmall_nbits = max(Gsmall_nbits, maximum(nbits, Gsmall[k]))
  end

  n = nrows(C.G[1])

  Gsmall_nbits = Gsmall_nbits + 1 # for the sign

  vectors_nbits = 0

  nrows_nbits = nbits(size(Gsmall[1], 1))

  bitbound = Int == Int64 ? 64 : 32

  abs_maxbits_vectors = Int == Int64 ? 30 : 15

  tmp = Vector{Int}(undef, n)

  for i in 1:length(V)
    # First canonicalize them
    cand = V[i]
    v = cand[1]
    vectors_nbits = max(vectors_nbits, maximum(nbits, v) + 1)
    k = 1

    while iszero(v[k])
      k += 1
    end

    if v[k] < 0
      v .*= -1
    end

    if vectors_nbits > abs_maxbits_vectors
      return false, Csmall
    end

    if Gsmall_nbits + vectors_nbits + nrows_nbits + 1 > bitbound
      return false, Csmall
    end

    _v = Vector{Int}(undef, n)

    for i in 1:n
      _v[i] = v[i]
    end

    w = Vector{Int}(undef, r)
    w[1] = cand[2]
    for k in 2:r
      w[k] = _norm(_v, Gsmall[k], tmp)
    end

    lengths[i] = w
    vectors[i] = _v
  end

  V = VectorList(vectors, lengths, use_dict)

  Csmall.V = V

  Csmall.prime = next_prime(2^(vectors_nbits + 1) + 1)

  Csmall.G = Matrix{Int}[Matrix{Int}(g) for g in C.G]
  Csmall.Gtr = Matrix{Int}[transpose(g) for g in Gsmall]
  Csmall.dim = n
  Csmall.issymmetric = C.issymmetric
  Csmall.operate_tmp = zeros(Int, n)

  @assert C.issymmetric[1]

  # Compute the fingerprint
  if automorphism_mode
    @vprint :Lattice 1 "Computing fingerprint: "
    @vtime :Lattice 1 fingerprint(Csmall)
    @vprint :Lattice 1 "$(Csmall.fp_diagonal)\n"
  end

  if automorphism_mode
    # Find the standard basis vectors
    Csmall.std_basis = Vector{Int}(undef, dim(Csmall))
    z = zeros(Int, dim(Csmall))
    for i in 1:dim(Csmall)
      z[Csmall.per[i]] = 1
      k = find_point(z, V)
      Csmall.std_basis[i] = k
      z[Csmall.per[i]] = 0
    end
  end

  #

  Csmall.v = Vector{Matrix{Int}}(undef, length(C.G))

  # Here needs to be another overflow check
  @inbounds for i in 1:length(Csmall.G)
    A = zeros(Int, length(Csmall.V.vectors), dim(C))
    for j in 1:length(Csmall.V.vectors)
      for k in 1:dim(Csmall)
        A[j, k] = _dot_product_with_row(Csmall.V.vectors[j], Csmall.G[i], k)
      end
    end
    Csmall.v[i] = A
  end

  if false
    for i in 1:length(Csmall.G)
      for j in 1:length(Csmall.V.vectors)
        for k in 1:length(Csmall.V.vectors)
          @assert  _dot_product_with_row(Csmall.V.vectors[j], Csmall.v[i], k) == dot(reshape(Csmall.V[j], (1, dim(C))) * Csmall.G[i], Csmall.V[k])
        end
      end
    end
  end

  Csmall.g = Vector{Vector{Matrix{Int}}}(undef, dim(C))
  for i in 1:dim(Csmall)
    Csmall.g[i] = Matrix{Int}[]
  end
  Csmall.ng = zeros(Int, dim(Csmall))
  Csmall.nsg = zeros(Int, dim(Csmall))
  Csmall.orders = Vector{Int}(undef, dim(Csmall))

  # -Id is always an automorphism
  mid = zeros(Int, dim(Csmall), dim(Csmall))
  for i in 1:dim(Csmall)
    mid[i, i] = -1
  end
  Csmall.g[1] = Matrix{Int}[mid]
  Csmall.ng[1] = 1

  # Calculate orbit lengths

  nH = 0

  for i in 1:dim(Csmall)
    nH += Csmall.ng[i]
  end

  H = Vector{Matrix{Int}}(undef, nH)

  if automorphism_mode
    for i in 1:dim(Csmall)
      if Csmall.ng[i] > 0
        nH = 0
        for j in i:dim(Csmall)
          for k in 1:Csmall.ng[j]
            nH += 1
            H[nH] = Csmall.g[j][k]
          end
        end
        #@assert _orbitlen_naive(Csmall.std_basis[i], Csmall.fp_diagonal[i], H, nH, Csmall.V) == _orbitlen(Csmall.std_basis[i], Csmall.fp_diagonal[i], H, nH, Csmall.V)
        Csmall.orders[i] = _orbitlen(Csmall.std_basis[i], Csmall.fp_diagonal[i], H, nH, Csmall.V, Csmall)
      else
        Csmall.orders[i] = 1
      end
    end
  end

  return true, Csmall
end

#function compute_scpcomb(C::ZAutLatCtx, depth::Int)
#  comb = Vector{SCPComp}(undef, dim(C))
#  for i in 1:dim(C)
#  end
#end

#	if (flags.DEPTH > 0)
#	{
#		if ((comb = (scpcomb*)malloc(dim * sizeof(scpcomb))) == 0)
#			exit (1);
#		for (i = 0; i < dim; ++i)
#		{
#/* compute the list of scalar product combinations and the corresponding
#   vector sums */
#			scpvecs(&comb[i].list, &sumveclist, i, fp.e, flags.DEPTH, V, F);
#/* compute a basis for the lattice that is generated by the vector sums and
#   a transformation matrix that expresses the basis in terms of the
#   vector sums */
#			base(&comb[i], &sumvecbase, sumveclist, F.A[0], dim);
#			if (flags.PRINT == 1)
#/* if the -P option is given, print the rank of the lattice generated by the
#   vector sums on level i on AUTO.tmp */
#			{
#				outfile = fopen("AUTO.tmp", "a");
#				fprintf(outfile, "comb[%d].rank = %d\n", i, comb[i].rank);
#				fclose(outfile);
#			}
#/* compute the coefficients of the vector sums in terms of the basis */
#			coef(&comb[i], sumvecbase, sumveclist, F.A[0], dim);
#			for (j = 0; j <= comb[i].list.n; ++j)
#				free(sumveclist[j]);
#			free(sumveclist);
#/* compute the scalar products of the base-vectors */
#			scpforms(&comb[i], sumvecbase, F);
#			for (j = 0; j < comb[i].rank; ++j)
#				free(sumvecbase[j]);
#			free(sumvecbase);
#		}
#	}

@doc Markdown.doc"""
    short_vectors(L, ub) -> Vector{Tuple{Vector{Int}, fmpq}}

Returns all tuples `(v, n)` such that `v G v^t = n <= ub`, where `G` is the
Gram matrix of `L`.

Note that the vectors are computed up to sign (so only one of `v` and `-v`
appears).
"""
short_vectors(L::ZLat, ub)

function short_vectors(L::ZLat, ub)
  _G = gram_matrix(L)
  return _short_vectors_gram(_G, ub)
end

function short_vectors(L::ZLat, lb, ub)
  _G = gram_matrix(L)
  return _short_vectors_gram(_G, lb, ub)
end

function shortest_vectors(L::ZLat, ::Type{Vector{Int}})
  _G = gram_matrix(L)
  min, V = _shortest_vectors_gram(_G)
  L.minimum = min
  return V
end

function shortest_vectors(L::ZLat)
  _G = gram_matrix(L)
  min, V = _shortest_vectors_gram(_G)
  W = Vector{fmpz_mat}(undef, length(V))
  n = rank(L)
  for i in 1:length(V)
    W[i] = matrix(FlintZZ, 1, n, V[i])
  end
  L.minimum = min
  return W
end

function minimum(L::ZLat)
  if !isdefined(L, :minimum)
    shortest_vectors(L)
  end

  return L.minimum
end

function kissing_number(L::ZLat)
  return 2 * length(shortest_vectors(L))
end

function compute_short_vectors(C::ZLatAutoCtx{Int, Matrix{Int}, Vector{Int}}, max = fmpz(-1))
  #V = enumerate_using_gram(G, R(max))

  if max == -1
    max = maximum(C.G[1][i, i] for i in 1:dim(C))
  end

  @vprint :Lattice 1 "Computing short vectors of actual length $max\n"
  V = _short_vectors_gram_integral(C.G[1], max)
  return V
end

function compute_short_vectors(C::ZLatAutoCtx, max::fmpz = fmpz(-1))
  #V = enumerate_using_gram(G, R(max))

  if max == -1
    max = maximum(C.G[1][i, i] for i in 1:dim(C))
  end
  @vprint :Lattice 1 "Computing short vectors of actual length $max\n"
  V = _short_vectors_gram_integral(C.G[1], max)
  n = ncols(C.G[1])
  C.V = Vector{fmpz_mat}(undef, length(V))
  C.V_length = Vector{Vector{fmpz}}(undef, length(V))
  for i in 1:length(V)
    z = Vector{fmpz}(undef, length(C.G))
    z[1] = V[i][2]
    m = matrix(FlintZZ, 1, n, V[i][1])
    mt = transpose(m)
    for k in 2:length(C.G)
      z[k] = (m * C.G[k] * mt)[1, 1]
    end
    C.V[i] = m
    C.V_length[i] = z
  end
  #@show length(C.V)
  C.max = max
  return C
end

function _get_vectors_of_length(G::Union{fmpz_mat, FakeFmpqMat}, max::fmpz)
  C = enum_ctx_from_gram(G)
  enum_ctx_start(C, max)
  res = Tuple{fmpz_mat, fmpz}[]
  while enum_ctx_next(C)
    push!(res, (deepcopy(C.x), (C.x * G * transpose(C.x))[1, 1]))
    push!(res, (-deepcopy(C.x), (C.x * G * transpose(C.x))[1, 1]))
  end
  return res
end

function _get_vectors_of_length(G::ZLat, max::fmpz)
  return _get_vectors_of_length(FakeFmpqMat(gram_matrix(G)), max)
end

function possible(C::ZLatAutoCtx, per, I, J)
  V = C.V.vectors
  W = C.V.lengths
  F = C.G
  Ftr = C.Gtr
  n = length(W)
  f = length(F)
  _issymmetric = C.issymmetric
  return possible(V, W, F, Ftr, _issymmetric, n, f, per, I, J)
end

function possible(V, W, F, Ftr, _issymmetric, n, f, per, I, J)
  count = 0

  tmp1 = zero(eltype(V[1]))
  tmp2 = zero(eltype(V[1]))

  for j in 1:n
    Wj = W[j]
    Vj = V[j]
    good_scalar = true
    good_length = true
    @inbounds for k in 1:f
      if Wj[k] != F[k][J, J]
        good_length = false
        break
      end
    end

    if !good_length
      continue
    end

    @inbounds for k in 1:f
      for i in 1:I
        if !(_dot_product_with_column!(tmp1, Vj, F[k], per[i], tmp2) == F[k][J, per[i]]) ||
              (!_issymmetric[k] && _dot_product_with_row!(tmp1, Vj, F[k], per[i], tmp2) != F[k][per[i], J])
          good_scalar = false
          break
        end
      end

      if !good_scalar
        break
      end
    end

    if good_length && good_scalar
      count = count + 1
    end

    if !good_length
      continue
    end

    # length is correct

    good_scalar = true

    @inbounds for k in 1:f
      for i in 1:I
        if !(_dot_product_with_column!(tmp1, Vj, F[k], per[i], tmp2) == -F[k][J, per[i]]) ||
              (!_issymmetric[k] && _dot_product_with_row!(tmp1, Vj, F[k], per[i], tmp2) != -F[k][per[i], J])
          good_scalar = false
          break
        end

      end

      if !good_scalar
        break
      end
    end

    if good_scalar
      count = count + 1
    end
  end
  return count
end

#function _dot_product(V::Vector, M, i)
#  z = zero(base_ring(V))
#  for j in 1:length(V)
#    z = z + V[j] * M[i, j]
#  end
#  return z
#end
#
#function _dot_product(V::fmpz_mat, M, i)
#  z = zero(base_ring(V))
#  for j in 1:length(V)
#    z = z + V[1, j] * M[i, j]
#  end
#  return z
#end

# a permutation per for the
#	order of the basis-vectors is chosen
#	such that in every step the number of
#	possible continuations is minimal,
#	for j from per[i] to per[dim-1] the
#	value f[i][j] in the fingerprint f is
#	the number of vectors, which have the
#	same scalar product with the
#	basis-vectors per[0]...per[i-1] as the
#	basis-vector j and the same length as
#	this vector with respect to all
#	invariant forms

function fingerprint(C::ZLatAutoCtx)
  V = C.V
  n = dim(C)
  k = length(C.G)
  per = Vector{Int}(undef, n)
  for i in 1:n
    per[i] = i
  end

  fp = zeros(Int, n, n)

  # fp[1, i] = # vectors v such that v has same length as b_i for all forms
  @inbounds for i in 1:n
    for j in 1:length(V)
      good = true
      cvl = @inbounds V.lengths[j]
      for l in 1:k
        if cvl[l] != C.G[l][i, i]
          good = false
          break
        end
      end

      if good
        fp[1, i] += 2 # the negative vector also has the correct length
      end
    end
  end

  @inbounds for i in 1:(n - 1)
    # Find the minimal non-zero entry in the i-th row
    mini = i
    for j in (i+1):n
      if fp[i, per[j]] < fp[i, per[mini]]
        mini = j
      end
    end

    per[mini], per[i] = per[i], per[mini]

    # Set entries below the minimal entry to zero
    for j in (i + 1):n
      fp[j, per[i]] = 0
    end

    # Now compute row i + 1

    for j in (i + 1):n
      fp[i + 1, per[j]] = possible(C, per, i, per[j])
    end
  end

  # Extract the diagonal

  res = Vector{Int}(undef, n)

  @inbounds for i in 1:n
    res[i] = fp[i, per[i]]
  end

  C.fp = fp

  C.fp_diagonal = res

  C.per = per

  return per, fp, res
end

# computes min(#O, orblen), where O is the orbit of pt under the first nG matrices in G
function _orbitlen(point::Int, orblen::Int, G::Vector{T}, nG, V, C) where {T}
  n = length(V)
  orb = ones(Int, orblen)
  orb[1] = point
  flag = zeros(Bool, 2*n + 1)
  flag[point + n + 1] = true
  # if flag[i + n+1] = 1, -n <= i <= n, then the point i is already in the orbit
  #flag = zero_Flv(2*n + 1)+n+1;
  #orb  = zero_Flv(orblen);
  #orb[1] = pt;
  #flag[pt] = 1;
  len = 1
  cnd = 1
  #@show G, nG
  while cnd <= len && len < orblen
    i = 1
    while i <= nG && len < orblen
      imag = _operate(orb[cnd], G[i], V, C.operate_tmp)
      if !flag[imag + n + 1]
        # the image is a new point in the orbit
        len += 1
        orb[len] = imag
        flag[imag + n + 1] = true
      end
      i += 1
    end
    cnd += 1
  end
  return len
end

function _operate(point, A::Matrix{Int}, V)
  return _operate(point, A, V, zeros(Int, size(A, 2)), sorted)
end

function _operate(point, A::fmpz_mat, V)
  return _operate(point, A, V, zero_matrix(FlintZZ, 1, ncols(A)))
end

Base.replace!(::typeof(-), m::fmpz_mat) = -m

function _operate(point, A, V, tmp)
# 	V.v is a sorted list of length V.n of vectors
#				of dimension V.dim, the number of V.v[nr]*A in
#				the list is returned, where a negative number
#				indicates the negative of a vector
  tmp = _vec_times_matrix!(tmp, V[point], A)
  #w = V[abs(point)] * A
  #if point < 0
  #  if tmp isa fmpz_mat
  #    for i in 1:ncols(tmp)
  #      tmp[1, i] = -tmp[1, i]
  #    end
  #  else
  #    tmp .*= -1 # tmp = -tmp
  #  end
  #end
  k = find_point(tmp, V)
  #@assert V[k] == tmp
  return k
end

function _find_point(w::Vector{Int}, V)
  positive = false
  for k in 1:length(w)
    if !iszero(w[k])
      positive = w[k] > 0
      break
    end
  end
  if positive
    if sorted
      k = searchsortedfirst(V, w)
    else
      k = findfirst(isequal(w), V)
    end
    @assert k !== nothing
    return k
  else
    w .*= -1 # w = -w
    if sorted
      k = searchsortedfirst(V, w)
    else
      k = findfirst(isequal(w), V)
    end
    @assert k !== nothing
    w .*= -1 # w = -w
    return -k
  end
end

function _find_point(w::fmpz_mat, V)
  positive = false
  for k in 1:length(w)
    if !iszero(w[1, k])
      positive = w[1, k] > 0
      break
    end
  end
  if positive
    if sorted
      k = searchsortedfirst(V, w)
    else
      k = findfirst(isequal(w), V)
    end
    @assert k !== nothing
    return k
  else
    mw = -w
    if sorted
      k = searchsortedfirst(V, mw)
    else
      k = findfirst(isequal(mw), V)
    end
    @assert k !== nothing
    return -k
  end
end

function _orbitlen_naive(point::Int, orblen::Int, G::Vector{fmpz_mat}, nG::Int, V)
  working_list = Int[point]
  orbit = Int[point]
  while !isempty(working_list)
    current_point = pop!(working_list)
    for i in 1:nG
      if current_point < 0
        new_point_coord = -V[abs(current_point)] * G[i]
      else
        new_point_coord = V[current_point] * G[i]
      end
      new_point = find_point(new_point_coord, V)
      if !(new_point in orbit)
        push!(orbit, new_point)
        push!(working_list, new_point)
      end
    end
  end
  return min(orblen, length(orbit))
end

function auto(C::ZLatAutoCtx{S, T, U}) where {S, T, U}
  dim = Hecke.dim(C)

  candidates = Vector{Vector{Int}}(undef, dim) # candidate list for the image of the i-th basis vector

  for i in 1:dim
    candidates[i] = zeros(Int, C.fp_diagonal[i])
  end

  x = Vector{Int}(undef, dim)
  bad = Vector{Int}(undef, 2 * length(C.V))

  sta = 1
  for step in sta:dim
    @vprint :Lattice 1 "Entering step $step\n"
    nH = 0
    for i in step:dim
      nH += C.ng[i]
    end
    #@show nH
    H = Vector{T}(undef, nH)
    nH = 0
    @inbounds for i in step:dim
      for j in 1:C.ng[i]
        nH += 1
        #@show C.g[i]
        H[nH] = C.g[i][j]
      end
    end
    @inbounds for i in 1:2*length(C.V)
      bad[i] = 0
    end
    nbad = 0
    @inbounds for i in 1:(step - 1)
      x[i] = C.std_basis[i]
    end
    #@show C.fp_diagonal[step]
    #@show candidates
    if C.fp_diagonal[step] > 1
      nC = cand(candidates[step], step, x, C, 0)#comb)
    else # there is only one candidates
      candidates[step] = Int[C.std_basis[step]]
      nC = 1
    end
    #@show nC
    #@show candidates
    orb = orbit(C.std_basis[step], 1, H, nH, C.V, C)
    C.orders[step] = length(orb)
    # delete the orbit of the step-th basis vector from the candidates
    #nC = delete(candidates[step], nC, orb, C.orders[step])
    setdiff!(candidates[step], orb)
    nC = length(candidates[step])
    #@show step, nC
    while nC > 0 && ((im = candidates[step][1]) != 0)
      @vprint :Lattice 1 "Step $(step), number of candidates left $(nC)\n"
      #@show im
      found = 0
      # try C.V[im] as the image of the step-th basis vector
      x[step] = im
      if step < dim
        #@show candidates
        if cand(candidates[step + 1], step + 1, x, C, 0) == C.fp_diagonal[step + 1]
          #@show candidates
          #@show "right before aut"
          #@show step + 1, x
          found = aut(step + 1, x, candidates, C, 0)#comb)
        else
          found = 0
        end
      else
        found = 1
      end

      #@show found

      if found == 0
        # x[1],...,x[step] cannot be continued
        oc = orbit(im, 1, H, nH, C.V, C)
        # delete the orbit of im from the candidates for x[step]
        #
        # This could go very bad ...
        candidates[step] = setdiff!(candidates[step], oc)
        nC = length(candidates[step])
        #nC = delete(candidates[step], nC, oc, noc)
        #empty!(oc)
        nbad += 1
        bad[nbad] = im
      else
        #@show x, step
        # a new generator has been found
        C.ng[step] += 1
        # append the new generator to C.>g[step]
        #@show "================================"
        ##@show C.g, step
        Gstep = resize!(C.g[step], C.ng[step])
        ##@show C.g, step
        matgen(x, dim, C.per, C.V)
        Gstep[C.ng[step]] = matgen(x, dim, C.per, C.V)
        C.g[step] = Gstep
        nH += 1
        H = Vector{T}(undef, nH)
        nH = 0
        @inbounds for i in step:dim
          for j in 1:C.ng[i]
            nH += 1
            H[nH] = C.g[i][j]
          end
        end
        # compute the new orbit of std_basis[step]
        orb = orbit(C.std_basis[step], 1, H, nH, C.V, C)
        C.orders[step] = length(orb)
        # delete the orbit from the candidates for x[step]
        setdiff!(candidates[step], orb)
        nC = length(candidates[step])
        #nC = delete(candidates[step], nC, orb, C.orders[step])
        # compute the new orbit of the vectors, which could be continued to an automorphism
        oc = orbit(bad, nbad, H, nH, C.V, C)
        # delete the orbit from the candidates
        setdiff!(candidates[step], oc)
        nC = length(candidates[step])
        #nC = delete(candidates[step], nC, oc, noc)
        #empty!(oc)
      end
    end
    if step == sta
      # test, whether on step flags.STAB some generators may be omitted
      tries = C.nsg[step]
      while tries <= C.ng[step]
      #for tries in C.nsg[step]:C.ng[step]
        nH = 0
        for j in 1:(tries-1)
          nH += 1
          H[nH] = C.g[step][j]
        end
        for j in (tries + 1):(C.ng[step]-1)
          nH += 1
          H[nH] = C.g[step][j]
        end
        for i in (step + 1):dim
          for j in 1:C.ng[i]
            nH += 1
            H[nH] = C.g[i][j]
            if false #_orbitlen(C.std_basis[step], C.orders[step], H, nH, C.V) == C.orders[step]
              # /* the generator g[step][tries] can be omitted */
              C.ng[step] -= 1
              for i in tries:(C.ng[step] - 1)
                C.g[step][i] = C.g[step][i + 1]
              end
              tries -= 1
            end
          end
        end
        tries += 1
      end
    end
    if step < dim && C.orders[step] > 1
     # /* calculate stabilizer elements fixing the basis-vectors
     #    C.std_basis[1]...fp.e[step] */
      stab(step, C)
    end
  end
  return _get_generators(C)
end

function _get_generators(C::ZLatAutoCtx{S, T, U}) where {S, T, U}
  # Extract generators

  gens = T[]

  orde = prod(fmpz.(C.orders))

  for i in 1:dim(C)
    for j in (C.nsg[i] + 1):C.ng[i]
      push!(gens, C.g[i][j])
    end
  end

  return gens, orde
end

function aut(step, x, candidates, C, comb)
  dim = Hecke.dim(C)
  found = 0
  x[step + 1:length(x)] .= 0
  while candidates[step][1] != 0 && found == 0
    if step < dim
      x[step] = candidates[step][1]
      #/* check, whether x[1]...x[step] is a partial automorphism and compute the candidates for x[step + 1]
			if (cand(candidates[step + 1], step + 1, x, C, comb) == C.fp_diagonal[step + 1])
        found = aut(step + 1, x, candidates, C, comb)
      end
      if found == 1
        return found
      end
      orb = Int[x[step]]
      # delete the chosen vector from the list of candidates
      #delete(candidates[step], C.fp_diagonal[step], orb, 1)
      k = findfirst(isequal(x[step]), candidates[step])
      #setdiff!(candidates[step], orb)
      # This should be made faster to not always go to the end
      # Actually I should copy the delete function
      #@show candidates[step], k
      for i in (k + 1):length(candidates[step])
        candidates[step][i - 1] = candidates[step][i]
      end
      candidates[step][end] = 0
      #@show candidates[step]
    else
      x[dim] = candidates[dim][1]
      found = 1
    end
  end
  return found
end

function cand(candidates, I, x, C::ZLatAutoCtx{S, T, U}, comb) where {S, T, U}
  #@show candidates, I, x, C, comb
  DEP = 0 # this is bs
  dim = Hecke.dim(C)
  len = length(C.G) * DEP
  vec = Vector{S}(undef, dim)
  vec2 = Vector{S}(undef, dim)
  scpvec = Vector{Int}(undef, len)
  if I >= 2 && DEP > 0
    com = comb[I - 1]
    rank = com.rank
    n = com.list.n
    # xvec is the list of vector sums which are computed with respect to the partial basis in x
    xvec = Vector{Vector{fmpz}}(undef, n + 1)
    for i in 1:(n + 1)
      xvec[i] = Vector{fmpz}(undef, dim)
      for j in 1:dim
        xvec[i][j] = zero(fmpz)
      end
    end
#/* xbase should be a basis for the lattice generated by the vectors in xvec,
#   it is obtained via the transformation matrix comb[I-1].trans */
    #xbase = zero_matrix(FlintZZ, rank, dim)
    #Fxbase = zero_matrix(FlintZZ, rank, dim)
  end
  # candidates is the list for the candidates
  #@show C.fp_diagonal[I], length(candidates)
  for i in 1:C.fp_diagonal[I]
    candidates[i] = 0
  end

  nr = 0
  fail = 0
  @inbounds for j in 1:length(C.V)
    if fail != 0
      break
    end
    Vvj = C.V[j]
    okp = 0
    okm = 0
    for k in 1:len
      scpvek[k] = 0
    end
    #@show C.V[j]
    for i in 1:length(C.G)
      _issym = C.issymmetric[i]
      CAiI = C.G[i][C.per[I]]
      Cvi = C.v[i]
      #@show Cvi

    # vec is the vector of scalar products of V.v[j] with the first I base vectors
    #   x[1]...x[I]

      for k in 1:(I - 1)
        #@show x[k]
        xk = x[k]
        if xk > 0
          #vec[k] = _dot_product(Vvj, C.G[i], C.V[xk])
          vec[k] = _dot_product_with_row(Vvj, C.v[i], xk)
          #@assert _tutut == vec[k]
          if !_issym
            #vec2[k] = _dot_product(C.V[xk], C.G[i], Vvj)
            vec2[k] = _dot_product_with_row(C.V[xk], C.v[i], j)
            #@assert vec2[k] == _tutut2
          end
        else
          #vec[k] = -_dot_product(Vvj, C.G[i], C.V[-xk])
          vec[k] = -_dot_product_with_row(Vvj, C.v[i], -xk)
          #@assert _tutut == vec[k]

          if !_issym
            #vec2[k] = -_dot_product(C.V[-xk], C.G[i], Vvj)
            vec2[k] = -_dot_product_with_row(C.V[-xk], C.v[i], j)
            #@assert vec2[k] == _tutut2
          end
        end
      end

      good = true
      for k in 1:(I - 1)
        if vec[k] != C.G[i][C.per[I], C.per[k]] || (!_issym && vec2[k] != C.G[i][C.per[k], C.per[I]])
          good = false
          break
        end
      end

      #@show "pos", Vvj, good

      if good && C.V.lengths[j][i] == C.G[i][C.per[I], C.per[I]]
        # C.V[j] is a candidate for x[I] with respec to the form C.G[i]
        okp += 1
      end

      good = true
      for k in 1:(I - 1)
        if vec[k] != -C.G[i][C.per[I], C.per[k]] || (!_issym && vec2[k] != -C.G[i][C.per[k], C.per[I]])
          good = false
          break
        end
      end

      #@show "neg", Vvj, good

      if good && C.V.lengths[j][i] == C.G[i][C.per[I], C.per[I]]
        # C.V[j] is a candidate for x[I] with respec to the form C.G[i]
        #@show "here"
        okm += 1
      end

      if okp < i && okm < i
        break
      end

      #if I >= 2 && DEP > 0
      #  for k in I-1:-1:1
      #    if k <= I - 1 - DEP
      #      continue
      #    end
      #    scpvec[(i - 1) * DEP + I - k] = vec[k]
      #  end
      #end
    end

    #if I >= 2 && DEP > 0
    #  # check whether the scalar product combination scpvec is contained in the list comb[I - 1].list
    #  if all(iszero, scpvec)
    #    num = 0
    #  else
    #    num = find_vector(scpvec, com.list)
    #    sign = num > 0 ? 1 : -1
    #    num = sign * num
    #  end

    #  if num > n
    #    # scpvec is not found, hence x[1],...,x[I - 1] is not a partial automorphism
    #    fail = 1
    #  elseif num > 0
    #    # scpvec is found and the vector is added to the corresponding vector sum
    #    xnum = xvec[num]
    #    for k in 1:dim
    #      xnum[k] += sign * Vvj[k]
    #    end
    #  end
    #end

    if okp == length(C.G)
      # V.v[j] is a candidate for x[I]
      if nr < C.fp_diagonal[I]
        nr += 1
        candidates[nr] = j
      else
        # there are too many candidates
        fail = 1
      end
    end

    #@show nr

    #@show okm == length(C.G)

    if okm == length(C.G)
      # -V.v[j] is a candidate for x[I]
      if nr < C.fp_diagonal[I]
        nr += 1
        candidates[nr] = -j
      else
        # there are too many candidates
        fail = 1
      end
    end

    #@show nr
  end

  #@show fail

  if fail == 1
    nr = 0
  end

  if nr == C.fp_diagonal[I] && I >= 2 && DEP > 0
    # compute the basis of the lattice generated by the vectors in xvec via the transformation matrix comb[I - 1].trans
    for i in 1:rank
      comtri = com.trans[i]
      for j in 1:dim
        xbij = FlintZZ(0)
        for k in 1:(n+1)
          xbij += comtri[k] * xvec[k][j]
        end
        xbase[i, j] = xbij
      end
    end
  end

  if nr == C.fp_diagonal[I] && I >= 2 && DEP > 0
    for i in 1:length(C.G)
      if !(nr > 0)
        break
      end
      for j in 1:rank
        for k in 1:dim
          Fxbase[j, k] = _dot_product(xbase[j], C.A[i], k)
        end
      end
      for j in 1:rank
        if !(nr > 0)
          break
        end
        for k in 1:j
          if !(nr > 0)
            break
          end
          if _dot_product(xbase[j], Fxbase[k]) != com.F[i][j][k]
            # scalar product is wrong
            nr = 0
          end
        end
      end
    end
  end

  if nr == C.fp_diagonal[I] && I >= 2 && DEP > 0
    for i in 1:(n + 1)
      if !(nr > 0)
        break
      end
    end
    comcoi = com.coeff[i]
    for j in 1:dim
      vj = zero(fmpz)
      for k in 1:rank
        vj += comcoi[k] * xbase[k][j]
      end
      if vj != xvec[i][j]
        # entry wrong
        nr = 0
        break
      end
    end
  end
  return nr
end

function orbit(pt, npt, G, nG, V, C::ZLatAutoCtx{S, T, U}) where {S, T, U}
  # Assumes that V is sorted
  orb = Vector{Int}(undef, npt)
  n = length(V)
  flag = zeros(Bool, 2*n + 1)
  #/* if flag[i + length(V)] is true, then the point i is already in the orbit */
  for i in 1:npt
    orb[i] = pt[i]
    flag[pt[i] + n + 1] = true
  end
  norb = npt
  cnd = 1
  while cnd <= norb
    for i in 1:nG
      im = _operate(orb[cnd], G[i], V, C.operate_tmp)
      if !flag[im + n + 1]
        # this is a new point
        norb += 1
        push!(orb, im)
        flag[im + n + 1] = true
      end
    end
    cnd += 1
  end
  return orb
end

function stab(I, C::ZLatAutoCtx{SS, T, U}) where {SS, T, U}
  V = C.V
#     	computes the orbit of fp.e[I] under the
#				generators in G->g[I]...G->g[n-1] and elements
#				stabilizing fp.e[I],
#				has some heuristic break conditions,
#				the generators in G->g[i] stabilize
#				fp.e[0]...fp.e[i-1] but not fp.e[i],
#				G->ng[i] is the number of generators in G->g[i],
#				the first G->nsg[i] of which are elements which
#				are obtained as stabilizer elements in
#				<G->g[0],...,G->g[i-1]>, G->ord[i] is the orbit
#				length of fp.e[i] under
#				<G->g[i],...,G->g[n-1]>	*****/
#group	*G;
#fpstruct fp;
#veclist	V;
#{
#	int	*orb, len, cnd, tmplen;
#	int	**w, *flag, ***H, ***Hj, **S, **tmp, ***Ggj;
#	int	i, j, k, l, dim, im, nH, nHj, fail;
#	int	Maxfail, Rest;
#
#/* some heuristic break conditions for the computation of stabilizer elements:
#   it would be too expensive to calculate all the stabilizer generators, which
#   are obtained from the orbit, since this is highly redundant,
#   on the other hand every new generator which enlarges the group is much
#   cheaper than one obtained from the backtrack,
#   after Maxfail subsequent stabilizer elements, that do not enlarge the group,
#   Rest more elements are calculated even if they leave the group unchanged,
#   since it turned out that this is often useful in the following steps,
#   increasing the parameters will possibly decrease the number of generators
#   for the group, but will increase the running time,
#   there is no magic behind this heuristic, tuning might be appropriate */
  dim = Hecke.dim(C)
  n = length(V)
  Rest = 0
  for i in I:dim
    if C.fp_diagonal[i] > 1 && C.orders[i] < C.fp_diagonal[i]
      Rest += 1
    end
  end

  Maxfail = Rest

  for i in 1:dim
    if C.fp_diagonal[i] > 1
      Maxfail += 1
    end
  end

  nH = 0
  for i in I:dim
    nH += C.ng[i]
  end

  Hj = Vector{T}(undef, nH + 1)
  H = Vector{T}(undef, nH)

#/* H are the generators of the group in which the stabilizer is computed */

  k = 0
  for i in I:dim
    for j in 1:C.ng[i]
      k += 1
      H[k] = C.g[i][j]
    end
  end

  w = Vector{Vector{Int}}(undef, 2 * n + 1)
  orb = zeros(Int, 2 * n)
  flag = zeros(Bool, 2 * n + 1)

#/* in w[V.n+i] an element is stored that maps fp.e[I] on v[i] */
#/* orb contains the orbit of fp.e[I] */
#/* if flag[i + V.n] = 1, then the point i is already in the orbit */
#/* S is a matrix to hold a stabilizer element temporarily */

  #@show I
  orb[1] = C.std_basis[I]
  flag[orb[1] + n + 1] = true
  #@show orb[1] + n + 1
  w[orb[1] + n + 1] = Int[ C.std_basis[i] for i in 1:dim]
  cnd = 1
  len = 1
  fail = 0

  while cnd <= len && fail < Maxfail + Rest
    for i in 1:nH
      if fail > Maxfail + Rest
        break
      end

      if fail >= Maxfail
      #/* there have already been Maxfail successive failures, now a random generator
      #   is applied to a random point of the orbit to get Rest more stabilizer
      #   elements */
        cnd = rand(1:len)
        i = rand(1:nH)
      end
      #@show orb, flag
      #@show cnd
      im = _operate(orb[cnd], H[i], V, C.operate_tmp)
      #@show im
      #@show w
      if !flag[im + n + 1]
#/* a new element is found, appended to the orbit and an element mapping
#   fp.e[I] to im is stored in w[im+V.n] */
        len += 1
        #@show orb, len
        orb[len] = im
        flag[im + n + 1] = true
        #@show w[orb[cnd] + n + 1]
        #@show H[i]
        #@show Int[_operate(w[orb[cnd] + n + 1][j], H[i], V) for j in 1:dim]
        w[im + n + 1] = Int[_operate(w[orb[cnd] + n + 1][j], H[i], V, C.operate_tmp) for j in 1:dim]
      else
#/* the image was already in the orbit */
        j = I
        while j <= dim
          if _operate(w[orb[cnd] + n + 1][j], H[i], V, C.operate_tmp) == w[im + n + 1][j]
            break
          end
          j += 1
        end
#/* j is the first index where the images of the old and the new element
#   mapping e[I] on im differ */
        if j <= dim && (C.orders[j] < C.fp_diagonal[j]  || fail >= Maxfail)
#/* new stabilizer element S = w[orb[cnd]+V.n] * H[i] * (w[im+V.n])^-1 */
          S = stabil(w[orb[cnd] + n + 1], w[im + n + 1], C.per, H[i], V, C)
          #@show S
          Hj[1] = S
          nHj = 1
          for k in j:dim
            for l in 1:C.ng[k]
              nHj += 1
              Hj[nHj] = C.g[k][l]
            end
          end
          tmplen = _orbitlen(C.std_basis[j], C.fp_diagonal[j], Hj, nHj, V, C)
          if tmplen > C.orders[j] || fail >= Maxfail
#/* the new stabilizer element S either enlarges the orbit of e[j]
#   or it is one of the additional elements after MAXFAIL failures */
            C.orders[j] = tmplen
            C.ng[j] = C.ng[j] + 1
            C.nsg[j] = C.nsg[j] + 1
            #@show C.g[j]
            #@show C.nsg[j]
            insert!(C.g[j], C.nsg[j], S)
#/* the new generator is inserted as stabilizer element nr. nsg[j]-1 */
            nH += 1
            push!(H, S)
            if fail < Maxfail
              fail = 0
            else
              fail += 1
            end
            resize!(Hj, nH + 1)
#/* the new generator is appended to H */
#/* the number of failures is reset to 0 */
          else
#/* the new stabilizer element S does not enlarge the orbit of e[j] */
            fail += 1
          end
        else
          if j <= dim && fail < Maxfail || (j == dim && fail >= Maxfail)
            fail += 1
          end
        end
      #/* if S is the identity and fail < Maxfail, nothing is done */
      end
    end
    if fail < Maxfail
      cnd += 1
    end
  end
end

#void stabil(S, x1, x2, per, G, V)	/*****	x1 corresponds to an element X1
#						mapping some vector e on p1,
#						x2 to an element X2 mapping e on
#						p2 and G is a generator mapping
#						p1 on p2, then S = X1*G*X2^-1
#						stabilizes e	*****/
function stabil(x1, x2, per, G, V, C)
  dim = length(x1)
  XG = zero_matrix(FlintZZ, dim, dim)
  X2 = zero_matrix(FlintZZ, dim, dim)
  x = Vector{Int}(undef, dim)
  for i in 1:dim
    x[i] = _operate(x1[i], G, V, C.operate_tmp) # fmpz case
  end

  XG = matgen(x, dim, per, V)
  X2 = matgen(x2, dim, per, V)

  b, S = can_solve_with_solution(X2, XG, side = :left)
  return S
end

function stabil(x1, x2, per, G::Matrix{Int}, V, C)
  #@show x1, x2
  dim = length(x1)
  x = Vector{Int}(undef, dim)
  for i in 1:dim
    x[i] = _operate(x1[i], G, V, C.operate_tmp)
  end

  XG = matgen(x, dim, per, V)
  X2 = matgen(x2, dim, per, V)

  @hassert :Lattice 1 begin XGG = deepcopy(XG); X22 = deepcopy(X2); true end
  SS = zeros(Int, dim, dim)
  _psolve(SS, X2, XG, dim, C.prime)
  @hassert :Lattice 1 SS * X22 == XGG

  return SS
end

fmpz_mat(M::Matrix{Int}) = matrix(FlintZZ, M)

zero_matrix(Int, r, c) = zeros(Int, r, c)

base_ring(::Vector{Int}) = Int

function _one(::Type{Matrix{Int}}, n::Int)
  z = zeros(Int, n, n)
  for i in 1:n
    z[i, i] = 1
  end
  return z
end

_one(::Type{fmpz_mat}, n::Int) = identity_matrix(FlintZZ, n)

_zero(::Type{Matrix{Int}}, n::Int, m::Int) = zeros(Int, n, m)

_zero(::Type{fmpz_mat}, n::Int, m::Int) = zero_matrix(FlintZZ, n, m)

function matgen(x, dim, per, v)
#/*****	generates the matrix X which has as row
#					per[i] the vector nr. x[i] from the
#					list v	*****
  X = zero_matrix(base_ring(v[1]), dim, dim)
  #@show x
  for i in 1:dim
    xi = x[i]
    if x[i] > 0
      for j in 1:dim
        X[per[i], j] = v[xi][j]
      end
    else
      for j in 1:dim
        X[per[i], j] = -v[-xi][j]
      end
    end
  end
  return X
end

# Isomorphism computation

function _try_iso_setup_small(Gi, Go)
  Ci = ZLatAutoCtx(Gi)
  fl, Cismall = try_init_small(Ci, false)
  if fl
    Co = ZLatAutoCtx(Go)
    fl2, Cosmall = try_init_small(Co, true, fmpz(Cismall.max))
    if fl2
      return true, Cismall, Cosmall
    end
  end

  return false, Cismall, Cismall
end

function _iso_setup(Gi, Go)
  Ci = ZLatAutoCtx(Gi)
  Co = ZLatAutoCtx(Go)
  init(Ci, true)
  init(Co, false, Ci.max)
  return Ci, Co
end

function isometry(Ci::ZLatAutoCtx{SS, T, U}, Co::ZLatAutoCtx{SS, T, U}) where {SS, T, U}
  d = dim(Co)
  C = Vector{Vector{Int}}(undef, d)
  # I could actually also test the minimum
  if length(Ci.V) != length(Co.V)
    return false, _zero(T, 0, 0)
  end
  for i in 1:d
    C[i] = zeros(Int, Ci.fp_diagonal[i])
  end
  x = zeros(Int, d)
  # compute the candidates for x[0]
  H = Vector{T}(undef, sum(length(gg) for gg in Co.g))
  k = 1
  for i in 1:length(Co.g)
    for j in 1:length(Co.g[i])
      H[k] = Co.g[i][j]
      k += 1
    end
  end
  isocand(C[1], 1, x, Ci, Co)

  found = iso(1, x, C, Ci, Co, H)
  if found
    ISO = matgen(x, d, Ci.per, Co.V)
    for k in 1:length(Ci.G)
      ISO * Co.G[k] * ISO' == Ci.G[k]
    end
    return true, ISO
  else
    return false, _zero(T, 0, 0)
  end
end

function iso(step, x, C, Ci, Co, G)
  d = dim(Ci)
  found = false
  @vprint :Lattice "Testing $(length(C[step])) many candidates\n"
  while !isempty(C[step]) && C[step][1] != 0 && !found
    @vprint :Lattice "Doing step $step\n"
    if step < d
      # choose the image of the base vector nr. step
      x[step] = C[step][1]
        # check whether x[1]..x[step]
      nbc = isocand(C[step + 1], step + 1, x, Ci, Co)
      if nbc == Ci.fp_diagonal[step + 1]
        # go deeper in the recursion
        Maxfail = 0
        # determine the heuristic value of Maxfail for the break condition in isostab
        for i in 1:step
          if Ci.fp_diagonal[i] > 1
            Maxfail += 1
          end
        end
        for i in (step + 1):d
          if Ci.fp_diagonal[i] > 1
            Maxfail += 2
          end
        end
        H = isostab(x[step], G, Co, Maxfail)
        found = iso(step + 1, x, C, Ci, Co, H)
      end
      if found
        return found
      end
      # This is horrible
      # this is remove orb from C[step], and then adding 0's at the end to make
      # it again as big as in the beginning. This can be done more efficiently.
      nc = length(C[step])
      orb = orbit(x[step], 1, G, length(G), Co.V, Co)
      no = length(orb)
      setdiff!(C[step], orb)
      newnc = length(C[step])
      resize!(C[step], nc)
      for i in newnc + 1:nc
        C[step][i] = 0
      end
    else
      x[d] = C[d][1]
      found = true
    end
  end
  return found
end

function isostab(pt, G, C::ZLatAutoCtx{S, T, U}, Maxfail) where {S, T, U}
# computes the orbit of V.v[pt]
# under the generators
#	G[0],...,G[nG-1] and elements
#	stabilizing V.v[pt], which are
#	stored in H, returns the number
#	of generators in H
# a heuristic break condition for the computation of stabilizer elements:
# it would be too expensive to calculate all the stabilizer generators, which
# are obtained from the orbit, since this is highly redundant,
# on the other hand every new generator which enlarges the group reduces the
# number of orbits and hence the number of candidates to be tested,
# after Maxfail subsequent stabilizer elements, that do not enlarge the group,
# the procedure stops,
# increasing Maxfail will possibly decrease the number of tests,
# but will increase the running time of the stabilizer computation
# there is no magic behind the heuristic, tuning might be appropriate */
#
  nG = length(G)
  d = dim(C)
#/* H are the generators of the stabilizer of C.V[pt] */
  V = C.V
  H = Vector{T}(undef, 1)
  nH = 0
  n = length(C.V)
  w = Vector{T}(undef, 2 * n + 2)
  orb = zeros(Int, 2 * n)
  orblen = 1
  flag = zeros(Bool, 2*n + 1)
#/* if flag[i + n + 1] = 1, then the point i is already in the orbit */
  orb[1] = pt
  flag[orb[1] + n + 1] = 1
#/* w[pt+V.n] is the Identity */
  w[orb[1] + n + 1] = _one(T, d)
  cnd = 1
  len = 1
  fail = 0
#/* fail is the number of successive failures */
  #A = zero_matrix(FlintZZ, d, d)
  #B = zero_matrix(FlintZZ, d, d)
  while (cnd <= len && fail < Maxfail)
    for i in 1:nG
      if fail >= Maxfail
        break
      end
      #@show G, i
      #@show orb[cnd]
      im = _operate(orb[cnd], G[i], V, C.operate_tmp)
      #@show im
      if !flag[im + n  + 1]
#/* a new element is found, appended to the orbit and an element mapping
        len += 1
        orb[len] = im
        flag[im + n + 1] = true
        w[im + n + 1] = w[orb[cnd] + n + 1] * G[i]
#   V.v[pt] to im is stored in w[im+V.n] */
      else
#/* the image was already in the orbit */
        B = w[orb[cnd] + n + 1] * G[i]
        #@show B
        #@show w[im + n + 1]
#/* check, whether the old and the new element mapping pt on im differ */
        if B != w[im + n + 1]
#/* new stabilizer element H[nH] = w[orb[cnd]+V.n] * G[i] * (w[im+V.n])^-1 */
          H[nH + 1] = w[orb[cnd] + n + 1] * G[i] * inv(w[im + n + 1])
          #	psolve((*H)[nH], A, B, dim, V.prime);
          rpt = rand(1:(n + 1))
          templen = _orbitlen(rpt, 2*n, H, nH + 1, V)
          while templen < orblen
#/* the orbit of this vector is shorter than a previous one, hence choose a new
#   random vector */
            rpt = rand(1:(n + 1))
            tmplen = _orbitlen(rpt, 2*n, H, nH + 1, V)
          end
          if tmplen > orblen
#/* the new stabilizer element H[nH] enlarges the group generated by H */
            orblen = tmplen
            nH += 1
            resize!(H, nH)
            fail = 0
          else
            fail += 1
          end
        end
        # if H[nH]is the identity, nothing is done
      end
    end
    cnd += 1
  end
  resize!(H, nH)
  return H
end

function isocand(CI, I, x, Ci::ZLatAutoCtx{S, T, U}, Co) where {S, T, U}
  d = dim(Ci)
  n = length(Ci.V)
  @assert n == length(Co.V)
  # Do something with bacher polynomials ...
  vec = Vector{S}(undef, d)
  vec2 = Vector{S}(undef, d)
  for i in 1:Ci.fp_diagonal[I]
    CI[i] = 0
  end
  nr = 0
  fail = false
  for j in 1:n
    if fail
      break
    end
    Vvj = Co.V[j]
    okp = 0
    okm = 0
    # do something with scpvec
    for i in 1:length(Co.G)
      # GiI = Ci.G[i][Ci.per[I]] this is the Ci.per[I]-th row of Ci.G[i]
      Fvi = Co.v[i]
      # vec is the vector of scalar products of Co.v[j] with the first I base vectors
      # x[1]...x[I-1]
      for k in 1:(I - 1)
        xk = x[k]
        # TODO: Check for symmetry as in the cand function
        # TODO: Use dot_product_with_row
        if xk > 0
          vec[k] = _dot_product(Vvj, Co.G[i], Co.V[xk])
          vec2[k] = _dot_product(Co.V[xk], Co.G[i], Vvj)
          #vec[k] = _dot_product(Vvj, Fvi, -xk)
        else
          vec[k] = -_dot_product(Vvj, Co.G[i], Co.V[-xk])
          vec2[k] = -_dot_product(Co.V[-xk], Co.G[i], Vvj)
          #vec[k] = -(Vvj * Co.G[i] * Co.V[-xk]')[1, 1]
          #vec2[k] = -(Co.V[-xk] * Co.G[i] * Vvj')[1, 1]
          #vec[k] = -_dot_product(Vvj, Fvi, -xk)
        end
      end
      good = true
      for k in 1:(I - 1)
        if vec[k] != Ci.G[i][Ci.per[I], Ci.per[k]] || vec2[k] != Ci.G[i][Ci.per[k], Ci.per[I]]
          good = false
          break
        end
      end
      if good && Co.V.lengths[j][i] == Ci.G[i][Ci.per[I], Ci.per[I]]
        okp += 1
      end
      good = true
      for k in 1:(I - 1)
        if vec[k] != -Ci.G[i][Ci.per[I], Ci.per[k]] || vec2[k] != -Ci.G[i][Ci.per[k], Ci.per[I]]

          good = false
          break
        end
      end
      if good && Co.V.lengths[j][i] == Ci.G[i][Ci.per[I], Ci.per[I]]
        okm += 1
      end
      # do something with scpvec
    end
    # do something with scpvec and DEP
    if okp == length(Ci.G)
      if nr < Ci.fp_diagonal[I]
        nr += 1
        CI[nr] = j
      else
        fail = true
      end
    end
    if okm == length(Ci.G)
      if nr < Ci.fp_diagonal[I]
        nr += 1
        CI[nr] = -j
      else
        fail = true
      end
    end
  end
  if fail
    nr = 0
  end

  #if nr == Ci.fp_diagonal[I] # DEP
  # update the blabla
  return nr
end

function assert_auto(C, order)
  G, o = _get_generators(C)
  if o != order
    throw(error("Order $o. Expected $order"))
  end

  for g in G
    for U in C.G
      if g * U * g' != U
        throw(error("Not an isometry.\nElement:\n $g\nGram matrix:\n$U"))
      end
    end
  end
  return true
end

################################################################################
#
#  Rewrite
#
################################################################################

################################################################################
#
#  Computational kernels
#
################################################################################

function _dot_product_with_column!(t::fmpz, v::fmpz_mat, A::fmpz_mat, k::Int, tmp::fmpz)
  mul!(t, v[1, 1], A[1, k])
  for i in 2:ncols(v)
    mul!(tmp, v[1, i], A[i, k])
    addeq!(t, tmp)
  end
  return t
end

function _dot_product_with_column(v::fmpz_mat, A::fmpz_mat, k::Int, tmp::fmpz = zero(FlintZZ))
  t = zero(FlintZZ)
  t = _dot_product_with_column!(t, v, A, k, tmp)
  return t
end

function _dot_product_with_row!(t::fmpz, v::fmpz_mat, A::fmpz_mat, k::Int, tmp::fmpz)
  mul!(t, v[1, 1], A[k, 1])
  for i in 2:ncols(v)
    mul!(tmp, v[1, i], A[k, i])
    addeq!(t, tmp)
  end
  return t
end

function _dot_product_with_row(v::fmpz_mat, A::fmpz_mat, k::Int, tmp::fmpz = zero(FlintZZ))
  t = zero(FlintZZ)
  t = _dot_product_with_row!(t, v, A, k, tmp)
  return t
end

function _dot_product_with_column!(t::Int, v::Vector{Int}, A::Matrix{Int}, k::Int, tmp::Int)
  @inbounds t = v[1] * A[1, k]
  @inbounds for i in 2:length(v)
    t = t + v[i] * A[i, k]
  end
  return t
end

function _dot_product_with_column(v::Vector{Int}, A::Matrix{Int}, k::Int, tmp::Int = zero(Int))
  t = zero(Int)
  t = _dot_product_with_column!(t, v, A, k, tmp)
  #@show size(A)
  #@assert (reshape(v, (1, length(v))) * A[1:size(A, 1), k:k])[1, 1] == t
  return t
end

function _dot_product_with_row!(t::Int, v::Vector{Int}, A::Matrix{Int}, k::Int, tmp::Int)
  @inbounds t = v[1] * A[k, 1]
  @inbounds for i in 2:length(v)
    t = t + v[i] * A[k, i]
  end
  return t
end

function _dot_product_with_row(v::Vector{Int}, A::Matrix{Int}, k::Int, tmp::Int = zero(Int))
  t = zero(Int)
  return _dot_product_with_row!(t, v, A, k, tmp)
end

# Generic vector times matrix

_vec_times_matrix!(w, v, A) = mul!(w, v, A)

_vec_times_matrix(v::Vector{Int}, A) = _vec_times_matrix!(zeros(Int, size(A, 2)), v , A)

function _vec_times_matrix!(w::Vector{Int}, v::Vector{Int}, A::Matrix{Int})
  c = size(A, 2)
  r = size(A, 1)
  for i in 1:c
    @inbounds w[i] = v[1] * A[1, i]
    for j in 2:r
      @inbounds w[i] += v[j] * A[j, i]
    end
  end
  #@assert reshape(w, 1, length(w)) == reshape(v, 1, length(v)) * A
  return w
end

function _dot_product(v::Vector{Int}, M::Matrix{Int}, w::Vector{Int})
  z = M * w
  zz = dot(v, z)
  #@assert zz == (reshape(v, 1, length(v)) * M * w)[1, 1]
  return zz
end

function _norm(v::Vector{Int}, M::Matrix{Int}, tmp::Vector{Int} = Vector{Int}(undef, length(v)))
  LinearAlgebra.mul!(tmp, M, v)
  return dot(v, tmp)
end

function _norm(v::fmpz_mat, M::fmpz_mat, tmp::fmpz_mat = zero_matrix(FlintZZ, 1, ncols(v)))
  mul!(tmp, v, M)
  return (v * tmp')[1, 1]
end

function _dot_product(v::fmpz_mat, M::fmpz_mat, w::fmpz_mat)
  return (v * M * w')[1, 1]
end

#

function _pgauss(r, A, B, n, p)
  ainv = invmod(A[r, r], p)
  for j in (r + 1):n
    if A[r, j] % p != 0
      f = A[r, j] * ainv % p
      for i in (r+1):n
        A[i, j] = (A[i, j] - f * A[i, r]) % p
      end
      for i in 1:n
        B[i, j] = (B[i, j] - f * B[i, r]) % p
      end
      A[r, j] = 0
    end
  end
end

function _psolve(X, A, B, n, p)
  for i in 1:(n - 1)
    j = i
    while A[i, j] % p == 0 && j <= n
      j += 1
    end
    if j == n + 1
      throw(error("Not possible"))
    end

    if j != i
      for k in i:n
        A[k, i], A[k, j] = A[k, j], A[k, i]
      end
      for k in 1:n
        B[k, i], B[k, j] = B[k, j], B[k, i]
      end
    end
    _pgauss(i, A, B, n, p)
  end
  for i in 1:n
    for j in n:-1:1
      sum = 0
      for k in n:-1:(j+1)
        sum = (sum + X[i, k] * A[k, j]) % p
      end
      ainv = invmod(A[j, j], p)
      X[i, j] = ((B[i, j] - sum) * ainv) % p
      c = 2 * X[i, j]
      if 2*c > p
        X[i, j] -= p
      elseif c <= -p
        X[i, j] += p
      end
    end
  end
end

###########################################
#
#  isless for fmpz_mat (vectors)
#
##########################################

function _isless(x::fmpz_mat, y::fmpz_mat)
  i = 0
  c = ncols(x)
  while i < c
    i += 1
    if x[i] == y[i]
      continue
    else
      return x[i] < y[i]
    end
  end
  return false
end

# should do this more C style
max_nbits(v::fmpz_mat) = maximum([nbits(v[1, i]) for i in 1:ncols(v)])

# Some tests that I need to add:
#
# G = matrix(FlintZZ, 8, 8, [4, -2, 0, 0, 0, 0, 0, 1, -2, 2, -1, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 0, 0, 0, 0, 0, -1, 2, 0, 1, 0, 0, 0, 0, 0, 0, 2])
#
# C = Hecke.ZLatAutoCtx([G]); Hecke.compute_short_vectors(C);
#
# Hecke.fingerprint(C)
# reduce(hcat, [C.fp[:, i] for i in 1:8][C.per]) == [240 240 2160 240 240 240 240 240; 0 56 126 126 126 126 126 126; 0 0 27 27 72 72 72 72; 0 0 0 10 40 16 40 40; 0 0 0 0 8 8 24 24; 0 0 0 0 0 4 6 12; 0 0 0 0 0 0 3 6; 0 0 0 0 0 0 0 2]


text:
abstract type AbstractKernel end

include("gauss_kernel.jl")
include("rq_kernel.jl")
include("iso_period_kernel.jl")
include("linear_kernel.jl")
include("neural_kernel_network.jl")
include("composite_kernel.jl")

ext: jl

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
function diffevol_minimise(f::Function, X::IntervalBox{N, T}, constraints::Vector{Constraint{T}}; ibc_chnl = RemoteChannel(()->Channel{Tuple{IntervalBox{N,T}, Float64}}(0)),
               diffevol_chnl = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(0)), np = 10*N, debug = false, de_ind = true, iterations = 20) where{N, T}

   pop = SVector{N, T}[]

   for i in 1:np
      indv = [X[j].lo + (1-rand())*(X[j].hi - X[j].lo) for j in 1:N]                        #Initialsing Population
      push!(pop, SVector{N, T}(indv))
   end

   global_min = Inf
   x_best = pop[np]

   num_of_iterations = 0

   while true

      fac = 2*rand()
      ind = rand(1:N)
      cr = rand()
      pop_new = SVector{N, T}[]

      temp = global_min

      if !de_ind
         if isready(diffevol_chnl)
            (x_best, temp, new_box) = take!(diffevol_chnl)  # Receiveing best individual from diffevol_minimise
            if debug
               println("Box recevied from IBC: ", (x_best, temp))
            end
            if temp == Inf
               break
            end
            if new_box != nothing
               X = new_box
            end
            push!(pop, x_best)
            np = np + 1
         end
      end

      for i in 1:np

         u = generate_random(1, np, i)
         v = generate_random(1, np, i, u)
         w = generate_random(1, np, i, u, v)   # Choosing index of three different individuals, different from the index of that individual whose mutant vector is going to form.

         m = bound_ensure(pop[u] + fac * (pop[v] - pop[w]), pop[u], X) # Mutation : Mutant Vector is created

         for j in 1:N                        # Recombination or CrossOver :  Mutant vector is itself is modified by Crossover rate (CR)
            if j != ind
               if rand() > cr
                  m = setindex(m, pop[i][j], j)
               end
            end
         end

         (c1, c2) = (0, 0)

         for constraint in constraints
            if invokelatest(constraint.C, m) ∈ constraint.bound
               c1 = c1 + 1
            end
            if invokelatest(constraint.C, pop[i]) ∈ constraint.bound
               c2 = c2 + 1
            end
         end

         if (c1, c2) == (length(constraints), length(constraints))
            if f(m) < f(pop[i])
               push!(pop_new, m)
            else
               push!(pop_new, pop[i])
            end
         else
            if c2 <= c1
               push!(pop_new, m)
            else
               push!(pop_new, pop[i])
            end
         end

         if c1 == length(constraints) || c2 == length(constraints)
            if f(pop_new[i]) < global_min
               global_min = f(pop_new[i])
               x_best = pop_new[i]
            end
         end
      end

      if global_min < temp
         best_box = (IntervalBox(Interval.(x_best)), global_min)
         if !de_ind
            put!(ibc_chnl, best_box)   #sending the best individual to ibc_minimise
         end
         if debug
            println("Box send to IBC: ", best_box)
         end
      end

      pop = pop_new

      num_of_iterations = num_of_iterations + 1

      if de_ind
         if num_of_iterations == iterations
            return (global_min, IntervalBox(Interval.(x_best)))
         end
      end

   end

end

function diffevol_maximise(f::Function, X::IntervalBox{N, T}, constraints::Vector{Constraint{T}}; ibc_chnl = RemoteChannel(()->Channel{Tuple{IntervalBox{N,T}, Float64}}(0)),
               diffevol_chnl = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(0)), np = 10*N, debug = false, de_ind = true, iterations = 20) where{N, T}

   diffevol_minimise(x -> -f(x), X, constraints, ibc_chnl = ibc_chnl, diffevol_chnl = diffevol_chnl, np = np, debug = debug, de_ind = de_ind, iterations = iterations)
end

max_issues_repo_name: yashcodes/Charibde

max_forks_repo_name: yashcodes/Charibde


text:
using CSV, DataFrames, DataFramesMeta, Plots, Dates, DelimitedFiles
# download and unzip the data file, from
# https://realized.oxford-man.ox.ac.uk/images/oxfordmanrealizedvolatilityindices.zip
# read the data into a DataFrame
data = CSV.read("oxfordmanrealizedvolatilityindices.csv",DataFrame)
# select the desired columns
data = data[:,[:Column1, :Symbol, :close_price, :rv10, :bv]]
# set desired names
rename!(data, :Column1 => :Date)
rename!(data, :rv10 => :rv)
# select desired index
sp500 = @subset(data, :Symbol .== ".SPX")
# drop times from date by taking the leading 10 characters (yyyy-mm-dd)
sp500.Date = [sp500.Date[i][1:10] for i = 1:size(sp500,1)]
# compute percentage returns
sp500.rets = vcat(0.0,[100.0 *(log(sp500.close_price[i]) - log(sp500.close_price[i-1])) for i=2:size(sp500,1)])
# scale volatility measures to be compatible with percentage returns (log price mult. by 100)
sp500.rv = 10000.0 .* sp500.rv
sp500.bv = 10000.0 .* sp500.bv
# select the date range
sp500 = @subset(sp500, :Date .>= "2013-12-17") # these dates give 1000 obs.
sp500 = @subset(sp500, :Date .<= "2017-12-05")
# write out variables to plain text file
data = sp500[:,[:rets, :rv, :bv]]
data = Matrix{Float64}(data)
CSV.write("sp500.csv", sp500)
#writedlm("sp500.txt", data)
# some plots
#plot(sp500.Date, sp500.rets, tickfontsize=5, legend=false)
#savefig("returns.png")
#plot(sp500.Date, [sp500.rv, sp500.bv], tickfontsize=5, label=["rv" "bv"])
#savefig("volatility.png")

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
using MetabolicFluxAnalysis
using Requests
using MAT

function import_cobra( filename::String , modelname::String )

    mf = MAT.matopen( filename )
    mv = MAT.matread( filename )

    model     = mv[modelname]
    modelKeys = keys(model)

    species_1   = model["metNames"]
    fluxes_pre  = model["rxns"]

    S    = model["S"]
    revs = model["rev"]

    # create species db
    species        = Vector{String}(length(species_1))
    species_db     = Dict{String,Species}()
    for zi=1:length(species_1)
        species_db[species_1[zi]] = Species(species_1[zi],0)
        species[zi]               = species_1[zi]
    end

    # create fluxes
    fluxes      = Vector{String}()
    flux_db     = Dict{String,Flux}()

    for zi=1:length(fluxes_pre)

        idx_neg = find( S[:,zi] .< 0 )
        idx_pos = find( S[:,zi] .> 0 )
        e::LC = LC{String}(); for i_n in idx_neg; e[species[i_n]] = -S[i_n,zi]; end
        p::LC = LC{String}(); for i_p in idx_pos; p[species[i_p]] =  S[i_p,zi]; end
        if(isempty(e)); e["*"]=1.; end
        if(isempty(p)); p["*"]=1.; end

        if(revs[zi]>0)
            push!(fluxes,string(fluxes_pre[zi],"_fw"))
            push!(fluxes,string(fluxes_pre[zi],"_bw"))
            e2 = copy(p)
            p2 = copy(e)
            flux_db[string(fluxes_pre[zi],"_fw")] = Flux(fluxes_pre[zi],e,p,Dict{RAS,PAS}())
            flux_db[string(fluxes_pre[zi],"_bw")] = Flux(fluxes_pre[zi],e2,p2,Dict{RAS,PAS}())
        else
            push!(fluxes,fluxes_pre[zi])
            flux_db[fluxes_pre[zi]] = Flux(fluxes_pre[zi],e,p,Dict{RAS,PAS}())
        end
    end

    return Network(fluxes,species,flux_db,species_db)
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _eqpole(l::T, b::T, southpole::Bool) where {T<:AbstractFloat}
    sgn = southpole ? -1 : 1
    l = deg2rad(sgn*l)
    b = deg2rad(sgn*b)
    r = 18 * sqrt(2 * (1 - sin(b))) * 3.53553391
    return r .* sincos(l)
end

"""
    eqpole(l, b[; southpole = false]) -> x, y

### Purpose ###

Convert right ascension \$l\$ and declination \$b\$ to coordinate \$(x, y)\$
using an equal-area polar projection.

### Explanation ###

The output \$x\$ and \$y\$ coordinates are scaled to be in the range \$[-90, 90]\$
and to go from equator to pole to equator.  Output map points can be
centered on the north pole or south pole.

### Arguments ###

* `l`: longitude, scalar or vector, in degrees
* `b`: latitude, same number of elements as right ascension, in degrees
* `southpole` (optional boolean keyword): keyword to indicate that the plot is
  to be centered on the south pole instead of the north pole.  Default is
  `false`.

### Output ###

The 2-tuple \$(x, y)\$:

* ``x`` coordinate, same number of elements as right ascension, normalized to be
  in the range ``[-90, 90]``.
* ``y`` coordinate, same number of elements as declination, normalized to be
  in the range ``[-90, 90]``.

### Example ###

```jldoctest
julia> using AstroLib

julia> eqpole(100, 35, southpole=true)
(-111.18287262822456, -19.604540237028665)

julia> eqpole(80, 19)
(72.78853915267848, 12.83458333897169)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
eqpole(l::Real, b::Real; southpole::Bool=false) =
    _eqpole(promote(float(l), float(b))..., southpole)

function eqpole(l::AbstractArray{L}, b::AbstractArray{<:Real};
                southpole::Bool=false) where {L<:Real}
    @assert length(l) == length(b)
    typel = float(L)
    x = similar(l, typel)
    y = similar(l, typel)
    for i in eachindex(l)
        x[i], y[i] = eqpole(l[i], b[i], southpole=southpole)
    end
    return x, y
end


text:
module TrainingTests

using Photon, Test

function simple_conv_model()
    model = Sequential(
        Conv2D(16, 3, relu),
        Conv2D(32, 3, relu),
        MaxPool2D(),
        Dense(32, relu),
        Dense(10)
    )
    return model
end

function test_save_load()
    m = simple_conv_model()
    X = KorA(randn(Float32,28,28,1,4))
    pred = m(X)
    workout = Workout(m, MSELoss())
    f = saveWorkout(workout)
    workout2 = loadWorkout(f)
    rm(f, force=true)
    pred2 = workout2.model(X)
    @assert pred == pred2
end

function getdata(s=28)
    [(
        KorA(randn(Float32,s,s,1,16)),
        KorA(randn(Float32,10,16))
    ) for i=1:10]
end

function test_train()
    model = simple_conv_model()
    workout = Workout(model, MSELoss())

    data = getdata()
    fit!(workout, data, epochs=2)

    @test workout.epochs == 2
    @test workout.steps == (2 * length(data))
    @test hasmetric(workout, :loss)
end

function test_train_valid()
    model = simple_conv_model()
    workout = Workout(model, MSELoss())

    data_tr = getdata()
    data_val = getdata()
    fit!(workout, data_tr, data_val, epochs=2)

    @test hasmetric(workout, :loss)
    @test hasmetric(workout, :val_loss)
end

function test_channel()
    model = simple_conv_model()
    workout = Workout(model, MSELoss())

    minibatch = (KorA(randn(Float32,30,30,3,4)), KorA(randn(Float32,10,4)))

    function randomdata()
        Channel() do channel
            for i in 1:4
                put!(channel, minibatch)
            end
        end
    end

    fit!(workout, randomdata, epochs=2)
    @test workout.epochs == 2
    @test hasmetric(workout, :loss)
end

@testset "Training" begin
    resetContext()
    test_train()
    test_train_valid()
    test_channel()
    test_save_load()
end

end


text:
@testset "mean" begin
    n = 9
    @testset "Basic" begin
        test_rrule(mean, randn(n))
    end
    @testset "with dims kwargs" begin
        test_rrule(mean, randn(n); fkwargs=(;dims=1))
        test_rrule(mean, randn(n,4); fkwargs=(;dims=2))
    end
end

@testset "variation: $var" for var in (std, var)
    test_rrule(var, randn(3))
    test_rrule(var, randn(4, 5); fkwargs=(; corrected=false))
    test_rrule(var, randn(ComplexF64, 6))
    test_rrule(var, Diagonal(randn(6)))

    test_rrule(var, randn(4, 5); fkwargs=(; dims=1))
    test_rrule(var, randn(ComplexF64, 4, 5); fkwargs=(; dims=2, corrected=false))
    test_rrule(var, UpperTriangular(randn(5, 5)); fkwargs=(; dims=1))

    x = PermutedDimsArray(randn(3, 4, 5), (3, 2, 1))
    xm = mean(x; dims=(1, 3))
    test_rrule(var, x; fkwargs=(; dims=(1, 3), mean=xm))
end


text:
#!/usr/bin/env julia

using Pkg

cfgribjl_dir = normpath(joinpath(@__DIR__, "..", ".."))

Pkg.activate(cfgribjl_dir); print("\n")

println(
    """
           Running pre-push hook
    -------------------------------------
    |  Tests run here while GRIB.jl is  |
    | missing from the registry as that |
    |  breaks CI testing functionality  |
    -------------------------------------
    """
)

include(joinpath(cfgribjl_dir, "./test/runtests.jl"))


text:
# Bogumił Kamiński, 2019-03-25
# a drop-in replacement of sudoku_solve from sudoku1.jl

# can also work under Julia 0.6, JuMP 0.18.5 and Cbc 0.4.3

const JULIA1 = VERSION >= v"1.0"

if JULIA1
    using LinearAlgebra
end

using JuMP
using Cbc

function sudoku_solve(lines, idx)
    t = [lines[10idx-j][k] - '0' for j in 8:-1:0, k in 1:9]
    m = if JULIA1
            Model(with_optimizer(Cbc.Optimizer, logLevel=0))
        else
            Model(solver=CbcSolver())
        end
    @variable(m, x[1:9, 1:9, 1:9], Bin)
    for i in 1:9, j in 1:9
        @constraint(m, sum(x[i, j, :]) == 1)
        @constraint(m, sum(x[i, :, j]) == 1)
        @constraint(m, sum(x[:, i, j]) == 1)
        t[i, j] > 0 && @constraint(m, x[i, j, t[i, j]] == 1)
    end
    for i in 0:2, j in 0:2, k in 1:9
        @constraint(m, sum(x[3i .+ (1:3), 3j .+ (1:3), k]) == 1)
    end
    JULIA1 ? optimize!(m) : solve(m)
    xv = JULIA1 ? value.(x) : getvalue(x)
    v = 100dot(1:9, xv[1,1,:]) + 10dot(1:9, xv[1,2,:]) + dot(1:9, xv[1,3,:])
    println(lines[10idx-9], ":\t", v)
    return v
end

max_stars_repo_path: 7. Sudoku/sudoku2.jl

max_issues_repo_path: 7. Sudoku/sudoku2.jl

max_forks_repo_path: 7. Sudoku/sudoku2.jl


text:
abstract type AbstractImmutableVector end

struct Vec2 <: AbstractImmutableVector
	x::Float32
	y::Float32

	Vec2() = new()
	Vec2(x, y) = new(x, y)
end

struct Vec3 <: AbstractImmutableVector
	x::Float32
	y::Float32
	z::Float32

	Vec3() = new()
	Vec3(x, y, z) = new(x, y, z)
end

struct Vec4 <: AbstractImmutableVector
	x::Float32
	y::Float32
	z::Float32
	w::Float32

	Vec4() = new()
	Vec4(x, y, z, w) = new(x, y, z, w)
end

vec2array(v::T) where T = [v.x, v.y, v.z]
position_func(positionField::Symbol) = v->vec2array(getfield(v, positionField))

struct MatrixColumn3
	e1::Float32
	e2::Float32
	e3::Float32
end

struct MatrixColumn4
	e1::Float32
	e2::Float32
	e3::Float32
	e4::Float32
end

abstract type AbstractImmutableMatrix end

struct Matrix3 <: AbstractImmutableMatrix
	c1::MatrixColumn3
	c2::MatrixColumn3
	c3::MatrixColumn3

	Matrix4() = new()
end

struct Matrix4 <: AbstractImmutableMatrix
	c1::MatrixColumn4
	c2::MatrixColumn4
	c3::MatrixColumn4
	c4::MatrixColumn4

	Matrix4() = new()
end

import Base: size, eltype

eltype(::Type{V}) where V<:AbstractImmutableVector = Float32
size(::Type{Vec2}) = (2,)
size(::Type{Vec3}) = (3,)
size(::Type{Vec4}) = (4,)

isvector(t::DataType) = t <: AbstractImmutableVector
size(::Type{T}, i) where T<:AbstractImmutableVector = size(T)[i]

eltype(::Type{Matrix3}) = Float32
size(::Type{Matrix3}) = (3, 3)

eltype(::Type{Matrix4}) = Float32
size(::Type{Matrix4}) = (4, 4)

size(::Type{T}, i) where T<:AbstractImmutableMatrix = size(T)[i]

ismatrix(t::DataType) = t <: AbstractImmutableMatrix

mutable struct SamplerType{GLTYPE}
end

get_sampler_type(::Type{SamplerType{GLTYPE}}) where GLTYPE = GLTYPE

const gl2jlTypes =
	Dict{UInt16, DataType}([
		(GL_FLOAT, Float32),
		(GL_FLOAT_VEC2, Vec2),
		(GL_FLOAT_VEC3, Vec3),
		(GL_FLOAT_VEC4, Vec4),
		(GL_FLOAT_MAT3, Matrix3),
		(GL_FLOAT_MAT4, Matrix4),
		(GL_SAMPLER_2D, SamplerType{Int64(GL_SAMPLER_2D)})
	])

gl2jltype(glType::Integer) = gl2jlTypes[glType]

const jl2glTypes =
	Dict{DataType, UInt16}([
		(Float32, GL_FLOAT),
		(Float64, GL_DOUBLE),
		(Float16, GL_HALF_FLOAT),
		(UInt16,  GL_UNSIGNED_SHORT),
		(Int16,   GL_SHORT),
		(UInt8,   GL_UNSIGNED_BYTE),
		(Int8,	GL_BYTE),
		(UInt32,  GL_UNSIGNED_INT),
		(Int32,   GL_INT)
	])

jl2gltype(jlType::DataType) = jl2glTypes[jlType]

function typeelements(dataType::DataType)
	fieldNames = fieldnames(dataType)
	if isempty(fieldNames)
		# simple type
		return (dataType, 1)
	else
		elType = Nothing
		for t in dataType.types
			if elType == Nothing
				elType = t
			elseif elType != t
				error("typeelements: A field with non-uniform types found")
			end
		end
		return (elType, length(fieldNames))
	end
end

function set_array_field(vert::Vector{T}, field::Symbol, values::Array) where T
	fieldInd = findfirst(isequal(field), fieldnames(T))
	elType, elCount = typeelements(fieldtype(T, fieldInd))
	offs = fieldoffset(T, fieldInd)
	dst = convert(Ptr{elType}, pointer(vert)) + offs
	valueStride = div(length(values), length(vert))
	srcBase = 0
	for vertInd = 0:length(vert)-1
		for elInd = 1:elCount
			unsafe_store!(dst, values[srcBase+elInd], elInd)
		end
		dst += sizeof(T)
		srcBase += valueStride
	end
end

function set_array_fields(vert::Vector{T}, fieldArrays::Dict{Symbol, Array}) where T
	for field in fieldnames(T)
		if haskey(fieldArrays, field)
			set_array_field(vert, field, fieldArrays[field])
		else
			@info("GRU.set_array_fields(): Missing data for field $field")
		end
	end
end


text:
@testset "Session" begin
    session = BLPData.Session(service_download_timeout_msecs=2000)
    BLPData.get_opened_services_names(session) == BLPData.DEFAULT_SERVICE_NAMES
    @test BLPData.is_service_open(session, "refdata")
    @test session.handle != C_NULL
    BLPData.stop(session)
    BLPData.destroy!(session)
    @test session.handle == C_NULL
    @test_throws ErrorException BLPData.Session("refdata", port=9000)
end

SESSION = BLPData.Session()

@testset "Subscribe" begin
    sublist = BLPData.subscribe(SESSION, "//blp/mktdata/ticker/PETR4 BS Equity?fields=BID,ASK")
    BLPData.unsubscribe(SESSION, sublist)
end

@testset "Service" begin
    service = SESSION["refdata"]

    println("operation names for $(service.name)")
    println(BLPData.list_operation_names(service))

    @test BLPData.has_operation(service, "HistoricalDataRequest")
    @test !BLPData.has_operation(service, "HistoricalDatarequest")

    @testset "HistoricalDataRequest" begin
        op = service["HistoricalDataRequest"]
        @test op.name == "HistoricalDataRequest"
        println(op)
    end
end

@testset "generic request/response" begin
    security = "PETR4 BS Equity"
    fields = [ "PX_LAST", "VOLUME" ]
    date_start = Date(2020, 4, 1)
    date_end = Date(2020, 4, 16)

    queue, corr_id = BLPData.send_request(SESSION, "refdata", "HistoricalDataRequest") do req
        push!(req["securities"], security)
        append!(req["fields"], fields)
        req["startDate"] = date_start
        req["endDate"] = date_end
    end

    resp = BLPData.parse_response_as(Dict, queue, corr_id)
    @test haskey(resp[1], :securityData)
    @test haskey(resp[1][:securityData], :fieldData)
    field_data = resp[1][:securityData][:fieldData]
    @test haskey(field_data[1], :PX_LAST)
    @test isa(field_data[1][:PX_LAST], Number)
    println(resp)

    BLPData.purge(queue)
end

@testset "bdp" begin
    @testset "single field" begin
        result = BLPData.bdp(SESSION, "PETR4 BS Equity", "PX_LAST")
        @test haskey(result, :PX_LAST)
        @test isa(result[:PX_LAST], Number)
    end

    @testset "many fields" begin
        result = BLPData.bdp(SESSION, "PETR4 BS Equity", ["PX_LAST", "VOLUME"])
        @test haskey(result, :PX_LAST)
        @test haskey(result, :VOLUME)
        @test isa(result[:PX_LAST], Number)
        @test isa(result[:VOLUME], Number)
    end

    @testset "many securities" begin
        result = BLPData.bdp(SESSION, ["PETR4 BS Equity", "VALE3 BS Equity"], ["PX_LAST", "VOLUME"])

        for (k, v) in result
            @test k == "PETR4 BS Equity" || k == "VALE3 BS Equity"
            @test haskey(v, :PX_LAST)
            @test isa(v[:PX_LAST], Number)
            @test isa(v[:VOLUME], Number)
        end
    end

    @testset "Invalid security" begin
        invalid_security = "INVALID SECURITY NAME"
        @test_throws BLPData.BLPResultErrException BLPData.bdp(SESSION, invalid_security, "TICKER")

        try
            BLPData.bdp(SESSION, invalid_security, "TICKER")
        catch err
            @test err.err.category == "BAD_SEC"
            @test err.err.subcategory == "INVALID_SECURITY"
        end
    end
end

@testset "bdh" begin
    result = BLPData.bdh(SESSION, "IBM US Equity", ["PX_LAST", "VWAP_VOLUME"], Date(2020, 1, 2), Date(2020, 1, 30))
    df = DataFrame(result)
    @test DataFrames.names(df) == [ "date", "PX_LAST", "VWAP_VOLUME" ]
    @test size(df) == (20, 3)
    show(df)

    @testset "periodicity" begin
        df = DataFrame(BLPData.bdh(SESSION, "PETR4 BS Equity", "PX_LAST", Date(2018, 2, 1), Date(2020, 2, 10), periodicity="YEARLY"))
        @test DataFrames.names(df) == [ "date", "PX_LAST" ]
        @test size(df) == (2, 2)
    end

    @testset "options" begin
        options = Dict("periodicitySelection" => "YEARLY", "periodicityAdjustment" => "CALENDAR")
        df = DataFrame(BLPData.bdh(SESSION, "PETR4 BS Equity", "PX_LAST", Date(2018, 2, 1), Date(2020, 2, 10), options=options))
        @test DataFrames.names(df) == [ "date", "PX_LAST" ]
        @test size(df) == (2, 2)
    end

    @testset "Exceptions" begin
        @test_throws BLPData.BLPResultErrException BLPData.bdh(SESSION, "INVALID SECURITY Equity", "PX_LAST", Date(2020, 4, 1), Date(2020, 4, 30))
        @test_throws BLPData.BLPResultErrException BLPData.bdh(SESSION, "PETR4 BS Equity", ["PX_LAST", "INVALID FIELD"], Date(2020, 4, 1), Date(2020, 4, 30))

        @testset "NoUnwrap" begin
            bdh_result = BLPData.bdh(SESSION, "PETR3 BS Equity", ["PX_LAST", "VOLUMExxx"], Date(2020, 4, 1), Date(2020, 4, 30), error_handling=BLPData.NoUnwrap())
            @test ismissing(bdh_result.field_data_vec[1].result.VOLUMExxx)
            @test isa(bdh_result.field_exceptions[:VOLUMExxx], BLPData.FieldErr)
        end
    end

    @testset "historical price" begin
        ticker = "PETR4 BS Equity"
        fields = [ "PX_LAST", "TURNOVER", "PX_BID", "PX_ASK", "EQY_WEIGHTED_AVG_PX", "EXCHANGE_VWAP" ]
        options = Dict(
            "periodicityAdjustment" => "CALENDAR",
            "periodicitySelection" => "DAILY",
            "currency" => "BRL",
            "pricingOption" => "PRICING_OPTION_PRICE",
            "nonTradingDayFillOption" => "ACTIVE_DAYS_ONLY",
            "nonTradingDayFillMethod" => "NIL_VALUE",
            "adjustmentFollowDPDF" => false,
            "adjustmentNormal" => false,
            "adjustmentAbnormal" => false,
            "adjustmentSplit" => false
        )

        df = DataFrame(BLPData.bdh(SESSION, ticker, fields, Date(2019, 1, 1), Date(2019, 2, 10), options=options))
        @test DataFrames.names(df) == [ "date", "PX_LAST", "TURNOVER", "PX_BID", "PX_ASK", "EQY_WEIGHTED_AVG_PX", "EXCHANGE_VWAP" ]
        @test size(df) == (27, 7)
        show(df)
    end

    @testset "adjusted price" begin
        ticker = "PETR4 BS Equity"
        field = "PX_LAST"
        options = Dict(
            "periodicityAdjustment" => "CALENDAR",
            "periodicitySelection" => "DAILY",
            "currency" => "BRL",
            "pricingOption" => "PRICING_OPTION_PRICE",
            "nonTradingDayFillOption" => "ACTIVE_DAYS_ONLY",
            "nonTradingDayFillMethod" => "NIL_VALUE",
            "adjustmentFollowDPDF" => false,
            "adjustmentNormal" => true,
            "adjustmentAbnormal" => true,
            "adjustmentSplit" => true
        )

        df = DataFrame(BLPData.bdh(SESSION, ticker, field, Date(2019, 1, 1), Date(2019, 2, 10), options=options))
        @test DataFrames.names(df) == [ "date", "PX_LAST" ]
        @test size(df) == (27, 2)
        show(df)
    end
end

@testset "bdh_intraday_ticks" begin
    @testset "Single security one field" begin
        d0 = DateTime(2020, 4, 27, 13)
        d1 = DateTime(2020, 4, 27, 13, 5)
        res = BLPData.bdh_intraday_ticks(SESSION, "PETR4 BS Equity", ["ASK"], d0, d1)
        @test length(res) > 10
        @test length(res) < 30 # 20 max
        for t in res
            @test t.type == :ASK
        end
    end

    @testset "Single security multiple fields" begin
        d0 = DateTime(2020, 4, 27, 13)
        d1 = DateTime(2020, 4, 27, 13, 5)
        res = BLPData.bdh_intraday_ticks(SESSION, "PETR4 BS Equity", ["TRADE", "BID", "ASK"], d0, d1)
        @test length(res) > 10
        @test length(res) < 50
        show(DataFrame(res))
    end

    @testset "Multiple securities" begin
        d0 = DateTime(2020, 4, 27, 13)
        d1 = DateTime(2020, 4, 27, 13, 5)
        res = BLPData.bdh_intraday_ticks(SESSION, ["PETR4 BS Equity", "VALE3 BS Equity"], ["TRADE", "BID", "ASK"], d0, d1)
        @test haskey(res, "PETR4 BS Equity")
        @test haskey(res, "VALE3 BS Equity")
        @test length(res["PETR4 BS Equity"]) > 10
        @test length(res["PETR4 BS Equity"]) < 50
        @test length(res["VALE3 BS Equity"]) > 10
        @test length(res["VALE3 BS Equity"]) < 40
    end
end

@testset "bds" begin
    @testset "COMPANY_ADDRESS" begin
        df = DataFrame(BLPData.bds(SESSION, "PETR4 BS Equity", "COMPANY_ADDRESS"))
        @test DataFrames.names(df) == [ "Address" ]
        @test df[end, :Address] == "Brazil"
        show(df)
    end

    @testset "DVD_HIST_GROSS_WITH_AMT_STAT" begin
        df = DataFrame(BLPData.bds(SESSION, "PETR4 BS Equity", "DVD_HIST_GROSS_WITH_AMT_STAT"))
        @test DataFrames.names(df) == [ "Declared Date", "Ex-Date", "Record Date", "Payable Date", "Dividend Amount", "Dividend Frequency", "Dividend Type", "Amount Status" ]
    end
end

@testset "benchmarks" begin
    @testset "single ticker" begin
        println("single ticker benchmark")
        @time result = BLPData.bdh(SESSION, "PETR4 BS Equity", ["PX_LAST", "VWAP_VOLUME"], Date(2020, 1, 2), Date(2020, 1, 30))
    end

    tickers = [ "PETR4 BS Equity",
                "PETR3 BS Equity",
                "VALE3 BS Equity",
                "BBDC4 BS Equity",
                "BBDC3 BS Equity",
                "ITUB4 BS Equity",
                "B3SA3 BS Equity",
                "ABEV3 BS Equity",
                "BBAS3 BS Equity",
                "ITSA4 BS Equity",
                "MGLU3 BS Equity",
                "LREN3 BS Equity",
                "VIVT4 BS Equity"
             ]

    @testset "Async bds" begin
        println("Async bds benchmark")
        @time result = BLPData.bds(SESSION, tickers, "DVD_HIST_GROSS_WITH_AMT_STAT")
        @test isa(result, Dict)

        for ticker in tickers
            @test haskey(result, ticker)
        end
    end

    @testset "Async bdh many fields" begin
        println("Async bdh benchmark")
        @time result = BLPData.bdh(SESSION, tickers, ["PX_LAST", "VWAP_VOLUME"], Date(2020, 1, 2), Date(2020, 1, 30))
        @test isa(result, Dict)

        for ticker in tickers
            @test haskey(result, ticker)
            @test isa(result[ticker], Vector)

            for row in result[ticker]
                @test isa(row, NamedTuple)
                @test haskey(row, :PX_LAST)
                @test haskey(row, :VWAP_VOLUME)
                @test isa(row[:PX_LAST], Number)
                @test isa(row[:VWAP_VOLUME], Number)
            end
        end
    end

    @testset "Async bdh many fields" begin
        println("Async bdh benchmark")
        @time result = BLPData.bdh(SESSION, tickers, "PX_LAST", Date(2020, 1, 2), Date(2020, 1, 30))
        @test isa(result, Dict)

        for ticker in tickers
            @test haskey(result, ticker)
            @test isa(result[ticker], Vector)

            for row in result[ticker]
                @test isa(row, NamedTuple)
                @test haskey(row, :PX_LAST)
                @test isa(row[:PX_LAST], Number)
            end
        end
    end

    @testset "bdp multiple securities" begin
        println("bdp multiple securities")
        @time result = BLPData.bdp(SESSION, tickers, ["PX_LAST", "VOLUME"])
        @test isa(result, Dict)

        for ticker in tickers
            @test haskey(result, ticker)
            @test isa(result[ticker], NamedTuple)
            @test isa(result[ticker].PX_LAST, Number)
            @test isa(result[ticker].VOLUME, Number)
        end
    end
end

BLPData.stop(SESSION)

#=
single ticker benchmark
  0.294486 seconds (1.99 k allocations: 99.656 KiB)
Async bds benchmark
  0.584568 seconds (1.05 M allocations: 53.660 MiB, 1.98% gc time)
Async bdh benchmark
  0.334141 seconds (31.69 k allocations: 1.558 MiB, 4.23% gc time)
Async bdh benchmark
  0.333756 seconds (20.68 k allocations: 1.030 MiB)
bdp multiple securities
  0.418478 seconds (1.70 k allocations: 93.312 KiB)
=#

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module FixedPointNumbers

using Base: reducedim_initarray

import Base: ==, <, <=, -, +, *, /, ~, isapprox,
             convert, promote_rule, show, isinteger, abs, decompose,
             isnan, isinf, isfinite,
             zero, oneunit, one, typemin, typemax, floatmin, floatmax, eps, sizeof, reinterpret,
             float, trunc, round, floor, ceil, bswap,
             div, fld, rem, mod, mod1, fld1, min, max, minmax,
             rand

using Base: @pure

# T => BaseType
# f => Number of Bytes reserved for fractional part
abstract type FixedPoint{T <: Integer, f} <: Real end

export
    FixedPoint,
    Fixed,
    Normed,
    floattype,
# "special" typealiases
    # Q and U typealiases are exported in separate source files
# Functions
    scaledual

reinterpret(x::FixedPoint) = x.i
reinterpret(::Type{T}, x::FixedPoint{T,f}) where {T,f} = x.i

# construction using the (approximate) intended value, i.e., N0f8
*(x::Real, ::Type{X}) where {X<:FixedPoint} = X(x)

# comparison
==(x::T, y::T) where {T <: FixedPoint} = x.i == y.i
 <(x::T, y::T) where {T <: FixedPoint} = x.i  < y.i
<=(x::T, y::T) where {T <: FixedPoint} = x.i <= y.i
"""
    isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))

For FixedPoint numbers, the default criterion is that `x` and `y` differ by no more than `eps`, the separation between adjacent fixed-point numbers.
"""
function isapprox(x::T, y::T; rtol=0, atol=max(eps(x), eps(y))) where {T <: FixedPoint}
    maxdiff = T(atol+rtol*max(abs(x), abs(y)))
    rx, ry, rd = reinterpret(x), reinterpret(y), reinterpret(maxdiff)
    abs(signed(widen1(rx))-signed(widen1(ry))) <= rd
end
function isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))
    isapprox(promote(x, y)...; rtol=rtol, atol=atol)
end

# predicates
isinteger(x::FixedPoint{T,f}) where {T,f} = (x.i&(1<<f-1)) == 0

# traits
typemax(::Type{T}) where {T <: FixedPoint} = T(typemax(rawtype(T)), 0)
typemin(::Type{T}) where {T <: FixedPoint} = T(typemin(rawtype(T)), 0)
floatmin(::Type{T}) where {T <: FixedPoint} = eps(T)
floatmax(::Type{T}) where {T <: FixedPoint} = typemax(T)

widen1(::Type{Int8})   = Int16
widen1(::Type{UInt8})  = UInt16
widen1(::Type{Int16})  = Int32
widen1(::Type{UInt16}) = UInt32
widen1(::Type{Int32})  = Int64
widen1(::Type{UInt32}) = UInt64
widen1(::Type{Int64})  = Int128
widen1(::Type{UInt64}) = UInt128
widen1(::Type{Int128}) = Int128
widen1(::Type{UInt128}) = UInt128
widen1(x::Integer) = x % widen1(typeof(x))

const ShortInts = Union{Int8,UInt8,Int16,UInt16}
const LongInts = Union{UInt64, UInt128, Int64, Int128, BigInt}

"""
    floattype(::Type{T})

Return the minimum float type that represents `T` without overflow to `Inf`.

# Example

A classic usage is to avoid overflow behavior by promoting `FixedPoint` to `AbstractFloat`

```julia
julia> x = N0f8(1.0)
1.0N0f8

julia> x + x # overflow
0.996N0f8

julia> float_x = floattype(eltype(x))(x)
1.0f0

julia> float_x + float_x
2.0f0
```
"""
floattype(::Type{T}) where {T <: Real} = T # fallback
floattype(::Type{T}) where {T <: Union{ShortInts, Bool}} = Float32
floattype(::Type{T}) where {T <: Integer} = Float64
floattype(::Type{T}) where {T <: LongInts} = BigFloat
floattype(::Type{FixedPoint{T,f}}) where {T <: ShortInts,f} = Float32
floattype(::Type{FixedPoint{T,f}}) where {T <: Integer,f} = Float64
floattype(::Type{FixedPoint{T,f}}) where {T <: LongInts,f} = BigFloat
floattype(::Type{F}) where {F <: FixedPoint} = floattype(supertype(F))
floattype(x::FixedPoint) = floattype(typeof(x))

nbitsfrac(::Type{FixedPoint{T,f}}) where {T <: Integer,f} = f
nbitsfrac(::Type{F}) where {F <: FixedPoint} = nbitsfrac(supertype(F))

rawtype(::Type{FixedPoint{T,f}}) where {T <: Integer,f} = T
rawtype(::Type{F}) where {F <: FixedPoint} = rawtype(supertype(F))
rawtype(x::FixedPoint) = rawtype(typeof(x))

# This IOBuffer is used during module definition to generate typealias names
_iotypealias = IOBuffer()

# Printing. These are used to generate type-symbols, so we need them
# before we include any files.
function showtype(io::IO, ::Type{X}) where {X <: FixedPoint}
    print(io, typechar(X))
    f = nbitsfrac(X)
    m = sizeof(X)*8-f-signbits(X)
    print(io, m, 'f', f)
    io
end
function show(io::IO, x::FixedPoint{T,f}) where {T,f}
    show(io, round(convert(Float64,x), digits=ceil(Int,f/_log2_10)))
    get(io, :compact, false) || showtype(io, typeof(x))
end
const _log2_10 = 3.321928094887362

function Base.showarg(io::IO, a::Array{T}, toplevel) where {T<:FixedPoint}
    toplevel || print(io, "::")
    print(io, "Array{")
    showtype(io, T)
    print(io, ",$(ndims(a))}")
    toplevel && print(io, " with eltype ", T)
end

include("fixed.jl")
include("normed.jl")
include("deprecations.jl")
const UF = (N0f8, N6f10, N4f12, N2f14, N0f16)

eps(::Type{T}) where {T <: FixedPoint} = T(oneunit(rawtype(T)),0)
eps(::T) where {T <: FixedPoint} = eps(T)
sizeof(::Type{T}) where {T <: FixedPoint} = sizeof(rawtype(T))

# Promotions for reductions
const Treduce = Float64
if isdefined(Base, :r_promote)
    # Julia v0.6
    Base.r_promote(::typeof(+), x::FixedPoint{T}) where {T} = Treduce(x)
    Base.r_promote(::typeof(*), x::FixedPoint{T}) where {T} = Treduce(x)
    Base.reducedim_init(f::typeof(identity),
                   op::typeof(+),
                   A::AbstractArray{T}, region) where {T <: FixedPoint} =
                       Base.reducedim_initarray(A, region, zero(Treduce))
    Base.reducedim_init(f::typeof(identity),
                   op::typeof(*),
                   A::AbstractArray{T}, region) where {T <: FixedPoint} =
                       Base.reducedim_initarray(A, region, oneunit(Treduce))
else
    # Julia v0.7
    Base.add_sum(x::FixedPoint, y::FixedPoint) = Treduce(x) + Treduce(y)
    Base.reduce_empty(::typeof(Base.add_sum), ::Type{F}) where {F<:FixedPoint}  = zero(Treduce)
    Base.reduce_first(::typeof(Base.add_sum), x::FixedPoint)   = Treduce(x)
    Base.mul_prod(x::FixedPoint, y::FixedPoint) = Treduce(x) * Treduce(y)
    Base.reduce_empty(::typeof(Base.mul_prod), ::Type{F}) where {F<:FixedPoint} = one(Treduce)
    Base.reduce_first(::typeof(Base.mul_prod), x::FixedPoint)  = Treduce(x)
end

for f in (:div, :fld, :fld1)
    @eval begin
        $f(x::T, y::T) where {T <: FixedPoint} = $f(reinterpret(x),reinterpret(y))
    end
end
for f in (:rem, :mod, :mod1, :min, :max)
    @eval begin
        $f(x::T, y::T) where {T <: FixedPoint} = T($f(reinterpret(x),reinterpret(y)),0)
    end
end

# When multiplying by a float, reduce two multiplies to one.
# Particularly useful for arrays.
scaledual(Tdual::Type, x) = oneunit(Tdual), x
scaledual(b::Tdual, x) where {Tdual <: Number} = b, x
scaledual(Tdual::Type, x::Union{T,AbstractArray{T}}) where {T <: FixedPoint} =
    convert(Tdual, 1/oneunit(T)), reinterpret(rawtype(T), x)
scaledual(b::Tdual, x::Union{T,AbstractArray{T}}) where {Tdual <: Number,T <: FixedPoint} =
    convert(Tdual, b/oneunit(T)), reinterpret(rawtype(T), x)

@noinline function throw_converterror(::Type{T}, x) where {T <: FixedPoint}
    n = 2^(8*sizeof(T))
    bitstring = sizeof(T) == 1 ? "an 8-bit" : "a $(8*sizeof(T))-bit"
    io = IOBuffer()
    show(IOContext(io, :compact=>true), typemin(T)); Tmin = String(take!(io))
    show(IOContext(io, :compact=>true), typemax(T)); Tmax = String(take!(io))
    throw(ArgumentError("$T is $bitstring type representing $n values from $Tmin to $Tmax; cannot represent $x"))
end

rand(::Type{T}) where {T <: FixedPoint} = reinterpret(T, rand(rawtype(T)))
rand(::Type{T}, sz::Dims) where {T <: FixedPoint} = reinterpret(T, rand(rawtype(T), sz))

end # module


text:
"""
    lp(N::T) where {T<:AbstractEcologicalNetwork}

Uses label propagation to generate a first approximation of the modular
structure of a network. This is usually followed by the BRIM (`brim`) method.
This method supposedly performs better for large graphs, but we rarely observed
any differences between it and variations of BRIM alone on smaller graphs.

#### References

Liu, X., Murata, T., 2009. Community Detection in Large-Scale Bipartite
Networks, in: 2009 IEEE/WIC/ACM International Joint Conference on Web
Intelligence and Intelligent Agent Technology. Institute of Electrical &
Electronics Engineers (IEEE). https://doi.org/10.1109/wi-iat.2009.15
"""
function lp(N::T) where {T<:AbstractEcologicalNetwork}
  L = Dict([species(N)[i]=>i for i in 1:richness(N)])

  neighbors_b = Dict([s=>N[:,s] for s in species(N, dims=2)])
  neighbors_t = Dict([s=>N[s,:] for s in species(N, dims=1)])

  # Initial modularity
  imod = Q(N, L)
  amod = imod
  improved = true

  while improved
    update_t = shuffle(species(N; dims=1))
    update_b = shuffle(species(N; dims=2))

    for s1 in update_t
      linked = neighbors_t[s1]
      labels = [L[s2] for s2 in linked]
      if length(labels) > 0
        counts = Dict()
        (i->counts[i] = get(counts, i, 0) + 1).(labels)
        cmax = maximum(values(counts))
        ok_keys = keys(filter(k -> k.second==cmax, counts))
        if length(ok_keys) > 0
          newlab = StatsBase.sample(collect(ok_keys))
          L[s1] = newlab
        end
      end
    end

    for s2 in update_b
      linked = neighbors_b[s2]
      labels = [L[s1] for s1 in linked]
      if length(labels) > 0
        counts = Dict()
        (i->counts[i] = get(counts, i, 0) + 1).(labels)
        cmax = maximum(values(counts))
        ok_keys = keys(filter(k -> k.second==cmax, counts))
        if length(ok_keys) > 0
          newlab = StatsBase.sample(collect(ok_keys))
          L[s2] = newlab
        end
      end
    end

    # Modularity improved?
    amod = Q(N, L)
    imod, improved = amod > imod ? (amod, true) : (amod, false)
  end
  tidy_modules!(L)
  return (N, L)
end

"""
    salp(N::T; θ::Float64=1.0, steps::Int64=10_000, λ::Float64=0.999, progress::Bool=false) where {T <: BipartiteNetwork}

Label-propagation using simulated annealing. This function uses simulated
annealing to propagate labels from neighboring nodes. It accepts a network as
input. The schedule of the simulated annealing is linear: at step k+1, the
temperature is θλᵏ. The initial temperature has been picked so that after 100
timesteps, using the default λ, a move decreasing modularity by 0.05 (20% of the
theoretical maximum) is picked with a probability of 0.1.

Optional arguments regulating the behavior of the simulated annealing
routine are:

- `λ=0.999`, the rate of temperature decay
- `θ=0.002`, the initial temperature
- `steps=10_000`, the number of annealing steps to perform
- `progress=false`, whether to display an info message every 100 timesteps

The θ parameter can be picked using the following method: if we want to allow a
maximal loss of modularity of δ, after timestep k, with a decay parameter λ,
with a probability P, then θ = -δ/[λᵏ×ln(P)]⁻¹. By beibg more or less
restrictive on these parameters, the user can pick a value of θ for every
problem.

This function can work as a first step (like `lp`), but in explorations during
the development of the package, we found that `brim` was rarely (if ever) able
to optmize the output further. It can therefore be used on its own.
"""
function salp(N::T; θ::Float64=0.002, steps::Int64=10_000, λ::Float64=0.999, progress::Bool=false) where {T <: BipartiteNetwork}
  Y = copy(N)
  m = last(each_species_its_module(N))

  Q0 = Q(Y, m)
  for step in 1:steps
    temperature = θ*λ^(step-1)
    update_side = rand() < 0.5 ? 1 : 2
    updated_species = rand(species(Y; dims=update_side))
    original_module = m[updated_species]
    neighbors = update_side == 1 ? N[updated_species,:] : N[:,updated_species]
    modules = [get(m, n, 0) for n in collect(neighbors)]
    m[updated_species] = rand(modules)
    QR = Q(Y, m)
    Δ = Q0 - QR
    if rand() ≤ exp(-Δ/temperature)
      Q0 = QR
    else
      m[updated_species] = original_module
    end
    if progress && (step % 100 == 0)
      info("t: $(lpad(step, 9)) \t θ: $(round(temperature, 2)) \t Q: $(round(Q0, 3)) \t m: $(length(unique(values(m))))")
    end
  end
  EcologicalNetworks.tidy_modules!(m)
  return Y, m
end