text:
################################################################################
#
#       GrpAb/GrpAbFinGen.jl : Finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

export AbelianGroup, DiagonalGroup, issnf, ngens, nrels, rels, snf, isfinite,
       isinfinite, rank, order, exponent, istrivial, isisomorphic,
       direct_product, istorsion, torsion_subgroup, sub, quo, iscyclic,
       psylow_subgroup, issubgroup, abelian_groups, flat, tensor_product,
       dual, chain_complex, isexact, homology, free_resolution

import Base.+, Nemo.snf, Nemo.parent, Base.rand, Nemo.issnf

################################################################################
#
#  Functions for some abstract interfaces
#
################################################################################

elem_type(G::GrpAbFinGen) = GrpAbFinGenElem

elem_type(::Type{GrpAbFinGen}) = GrpAbFinGenElem

parent_type(::Type{GrpAbFinGenElem}) = GrpAbFinGen

##############################################################################
#
#  Constructors
#
##############################################################################
# We do we have AbelianGroup and DiagonalGroup?
@doc Markdown.doc"""
    AbelianGroup(M::fmpz_mat) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
function AbelianGroup(M::fmpz_mat; name :: String = "")
  G = GrpAbFinGen(M)
  if name != ""
    set_name!(G, name)
  end
  return G
end

@doc Markdown.doc"""
    AbelianGroup(M::Array{fmpz, 2}) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
function AbelianGroup(M::Array{fmpz, 2}; name :: String = "")
  G = AbelianGroup(matrix(FlintZZ, M))
  if name != ""
    set_name!(G, name)
  end
  return G
end

@doc Markdown.doc"""
    AbelianGroup(M::Array{Integer, 2}) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `ncols(M)` generators and each row of `M` describes one relation.
"""
function AbelianGroup(M::Array{T, 2}; name :: String = "") where T <: Integer
  G = AbelianGroup(matrix(FlintZZ, M))
  if name != ""
    set_name!(G, name)
  end
  return G
end

@doc Markdown.doc"""
    AbelianGroup(M::Array{fmpz, 1}) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `length(M)` generators and one relation.
"""
function AbelianGroup(M::Array{fmpz, 1}; name :: String = "")
  G = AbelianGroup(matrix(FlintZZ, 1, length(M), M))
  if name != ""
    set_name!(G, name)
  end
  return G
end

@doc Markdown.doc"""
    AbelianGroup(M::Array{Integer, 1}) -> GrpAbFinGen

Creates the abelian group with relation matrix `M`. That is, the group will
have `length(M)` generators and one relation.
"""
function AbelianGroup(M::Array{T, 1}; name :: String = "") where T <: Integer
  G = AbelianGroup(matrix(FlintZZ, 1, length(M), M))
  if name != ""
    set_name!(G, name)
  end
  return G
end

@doc Markdown.doc"""
    DiagonalGroup(M::fmpz_mat) -> GrpAbFinGen

Assuming that $M$ has only one row, this function creates the direct product of
the cyclic groups $\mathbf{Z}/m_i$, where $m_i$ is the $i$th entry of `M`.
"""
function DiagonalGroup(M::fmpz_mat; name :: String = "")
  if nrows(M) != 1
    error("The argument must have only one row")
  end

  N = zero_matrix(FlintZZ, ncols(M), ncols(M))
  for i = 1:ncols(M)
    N[i,i] = M[1, i]
  end
  if issnf(N)
    G = GrpAbFinGen(fmpz[M[1, i] for i = 1:ncols(M)])
    if ncols(M) != 0 && !iszero(M[1, ncols(M)])
      G.exponent = M[1, ncols(M)]
    end
  else
    G = GrpAbFinGen(N)
    if ncols(M) != 0
      res = lcm(fmpz[M[1, i] for i = 1:ncols(M)])
      if !iszero(res)
        G.exponent = res
      end
    end
  end
  name == "" || set_name!(G, name)
  return G
end

@doc Markdown.doc"""
    DiagonalGroup(M::Array{Union{fmpz, Integer}, 1}) -> GrpAbFinGen

Creates the direct product of the cyclic groups $\mathbf{Z}/m_i$,
where $m_i$ is the $i$th entry of `M`.
"""
function DiagonalGroup(M::Array{T, 1}; name :: String = "") where T <: Union{Integer, fmpz}
  N = zero_matrix(FlintZZ, length(M), length(M))
  for i = 1:length(M)
    N[i,i] = M[i]
  end
  if issnf(N)
    G = GrpAbFinGen(M)
  else
    G = GrpAbFinGen(N)
  end
  if !isempty(M)
    res = lcm(M)
    if !iszero(res)
      G.exponent = res
    end
  end
  name == "" || set_name!(G, name)
  return G
end

################################################################################
#
#  String I/O
#
################################################################################
function show(io::IO, A::GrpAbFinGen)
  @show_name(io, A)
  @show_special(io, A)

  if issnf(A)
    show_snf(io, A)
  else
    show_gen(io, A)
  end
end

function show_hom(io::IO, G::GrpAbFinGen)
  D = get_special(G, :hom)
  D === nothing && error("only for hom")
  print(io, "hom of ")
  print(IOContext(io, :compact => true), D)
end

function show_direct_product(io::IO, G::GrpAbFinGen)
  D = get_special(G, :direct_product)
  D === nothing && error("only for direct products")
  print(io, "direct product of ")
  show(IOContext(io, :compact => true), D)
end

function show_tensor_product(io::IO, G::GrpAbFinGen)
  D = get_special(G, :tensor_product)
  D === nothing && error("only for tensor products")
  print(io, "tensor product of ")
  show(IOContext(io, :compact => true), D)
end

function show_gen(io::IO, A::GrpAbFinGen)
  print(io, "(General) abelian group with relation matrix\n$(A.rels)")
  if isdefined(A, :snf_map)
    println(io, "\nwith structure of ", domain(A.snf_map))
  end
end

function show_snf(io::IO, A::GrpAbFinGen)
  if get(io, :compact, false)
    print(io, "Abelian group with structure: ")
  else
    print(io, "GrpAb: ")
  end
  show_snf_structure(io, A)
end

function show_snf_structure(io::IO, A::GrpAbFinGen, mul = "x")
  @assert issnf(A)
  len = length(A.snf)
  if len == 0
    print(io, "Z/1")
    return
  end

  if A.snf[1] == 0
    if len == 1
      print(io, "Z")
    else
      print(io, "Z^$(len)")
    end
    return
  end

  i = 1
  while i <= len
    inv = A.snf[i]
    j = 1
    while i + j <= len && inv == A.snf[i + j]
      j += 1
    end
    if iszero(inv)
      print(io, "Z")
    else
      print(io, "Z/$(inv)")
    end
    if j > 1
      print(io, "^$(j)")
    end
    if i + j - 1 < len
      print(io, " $mul ")
    end
    i += j
  end
end

################################################################################
#
#  Field access
#
################################################################################

@doc Markdown.doc"""
    issnf(G::GrpAbFinGen) -> Bool

Returns whether the current relation matrix of the group $G$ is in Smith
normal form.
"""
issnf(A::GrpAbFinGen) = A.issnf

@doc Markdown.doc"""
    ngens(G::GrpAbFinGen) -> Int

Returns the number of generators of $G$ in the current representation.
"""
function ngens(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return ncols(A.rels)
  end
end

@doc Markdown.doc"""
    nrels(G::GrpAbFinGen) -> Int

Returns the number of relations of $G$ in the current representation.
"""
function nrels(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return nrows(A.rels)
  end
end

@doc Markdown.doc"""
    rels(A::GrpAbFinGen) -> fmpz_mat

Returns the currently used relations of $G$ as a single matrix.
"""
function rels(A::GrpAbFinGen)
  if issnf(A)
    return rels_snf(A)
  else
    return rels_gen(A)
  end
end

function rels_gen(A::GrpAbFinGen)
  return A.rels
end

function rels_snf(A::GrpAbFinGen)
  M = zero_matrix(FlintZZ, ngens(A), ngens(A))
  for i = 1:ngens(A)
    M[i,i] = A.snf[i]
  end
  return M
end

################################################################################
#
#  Hermite normal form
#
################################################################################

function assure_has_hnf(A::GrpAbFinGen)
  if isdefined(A, :hnf)
    return nothing
  end
  if isdefined(A, :exponent)
    A.hnf = hnf_modular_eldiv(A.rels, A.exponent)
  else
    A.hnf = hnf(A.rels)
  end
  return nothing
end

################################################################################
#
#  Smith normal form
#
################################################################################

@doc Markdown.doc"""
    snf(A::GrpAbFinGen) -> GrpAbFinGen, Map

Returns a pair $(G, f)$, where $G$ is an abelian group in canonical Smith
normal form isomorphic to $G$ and an isomorphism $f : G \to A$.
"""
function snf(G::GrpAbFinGen)
  if isdefined(G, :snf_map)
    return domain(G.snf_map)::GrpAbFinGen, G.snf_map::GrpAbFinGenMap
  end

  if issnf(G)
    G.snf_map = GrpAbFinGenMap(G) # identity
    return G, G.snf_map::GrpAbFinGenMap
  end

  S, _, T = snf_with_transform(G.rels, false, true)

  return _reduce_snf(G, S, T, inv(T))
end

# For S in SNF with G.rels = U*S*T and Ti = inv(T) this removes
# the ones at the diagonal of S and constructs the homomorphism.
function _reduce_snf(G::GrpAbFinGen, S::fmpz_mat, T::fmpz_mat, Ti::fmpz_mat)

  d = fmpz[S[i,i] for i = 1:min(nrows(S), ncols(S))]

  while length(d) < ngens(G)
    push!(d, 0)
  end

  #s = Array{fmpz, 1}()
  s = fmpz[ d[i] for i in 1:length(d) if d[i] !=  1]
  #for i = 1:length(d)
  #  if d[i] != 1
  #    push!(s, d[i])
  #  end
  #end
  TT = zero_matrix(FlintZZ, nrows(T), length(s))
  j = 1
  for i = 1:length(d)
    if d[i] != 1
      for k=1:nrows(T)
        TT[k, j] = T[k, i]
      end
      j += 1
    end
  end

  TTi = zero_matrix(FlintZZ, length(s), nrows(T))

  j = 1
  for i = 1:length(d)
    if d[i] != 1
      for k=1:nrows(T)
        TTi[j, k] = Ti[i, k]
      end
      j += 1
    end
  end

  H = GrpAbFinGen(s)
  if !isempty(s) && !iszero(s[end])
    H.exponent = s[end]
    G.exponent = s[end]
  end
  mp = hom(H, G, TTi, TT, check = false)
  G.snf_map = mp
  return H, mp::GrpAbFinGenMap
end

################################################################################
#
#  Finiteness
#
################################################################################

@doc Markdown.doc"""
    isfinite(A::GrpAbFinGen) -> Bool

Returns whether $A$ is finite.
"""
isfinite(A::GrpAbFinGen) = issnf(A) ? isfinite_snf(A) : isfinite_gen(A)

isfinite_snf(A::GrpAbFinGen) = length(A.snf) == 0 || !iszero(A.snf[end])

isfinite_gen(A::GrpAbFinGen) = isfinite(snf(A)[1])

@doc Markdown.doc"""
    isinfinite(A::GrpAbFinGen) -> Bool

Returns whether $A$ is infinite.
"""
isinfinite(A::GrpAbFinGen) = !isfinite(A)

################################################################################
#
#  Rank
#
################################################################################

@doc Markdown.doc"""
    rank(A::GrpAbFinGenGen) -> Int

Returns the rank of $A$, that is, the dimension of the
$\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
"""
rank(A::GrpAbFinGen) = issnf(A) ? rank_snf(A) : rank_gen(A)

rank_snf(A::GrpAbFinGen) = length(findall(x -> x == 0, A.snf))

rank_gen(A::GrpAbFinGen) = rank(snf(A)[1])

################################################################################
#
#  Order
#
################################################################################

@doc Markdown.doc"""
    order(A::GrpAbFinGen) -> fmpz

Returns the order of $A$. It is assumed that $A$ is finite.
"""
order(A::GrpAbFinGen) = issnf(A) ? order_snf(A) : order_gen(A)

function order_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  return prod(A.snf)
end

order_gen(A::GrpAbFinGen) = order(snf(A)[1])

################################################################################
#
#  Exponent
#
################################################################################

@doc Markdown.doc"""
    exponent(A::GrpAbFinGen) -> fmpz

Returns the exponent of $A$. It is assumed that $A$ is finite.
"""
function exponent(A::GrpAbFinGen)
  if isdefined(A, :exponent)
    return A.exponent
  end
  if issnf(A)
    res = exponent_snf(A)
    if !iszero(res)
      A.exponent = res
    end
    return res
  else
    res = exponent_gen(A)
    if !iszero(res)
      A.exponent = res
    end
    return res
  end
end

function exponent_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  ngens(A)==0 && return fmpz(1)
  return A.snf[end]
end

exponent_gen(A::GrpAbFinGen) = exponent(snf(A)[1])

################################################################################
#
#  Trivial
#
################################################################################

@doc Markdown.doc"""
    istrivial(A::GrpAbFinGen) -> Bool

Checks if $A$ is the trivial group.
"""
istrivial(A::GrpAbFinGen) = isfinite(A) && isone(order(A))

################################################################################
#
#  Isomorphism
#
################################################################################

@doc Markdown.doc"""
    isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen) -> Bool

Checks if $G$ and $H$ are isomorphic.
"""
function isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen)
  b = filter(x -> x != 1, snf(G)[1].snf) == filter(x -> x != 1, snf(H)[1].snf)
  return b
end

################################################################################
#
#  Direct product
#
################################################################################
#TODO: check the universal properties here!!!
@doc Markdown.doc"""
    direct_product(G::GrpAbFinGen...; task::Symbol = :sum) -> GrpAbFinGen, GrpAbFinGenMap, GrpAbFinGenMap

Returns the direct product $D$ of the abelian groups $G_i$. {{{task}}} can be
":sum", ":prod", ":both" or ":none" and determines which canonical maps
are computed as well: ":sum" for the injections, ":prod" for the
   projections.
"""
function direct_product(G::GrpAbFinGen...
             ; add_to_lattice::Bool = false, L::GrpAbLattice = GroupLattice, task::Symbol = :sum)
  @assert task in [:prod, :sum, :both, :none]

  Dp = AbelianGroup(cat([rels(x) for x = G]..., dims = (1,2)))

  set_special(Dp, :direct_product =>G, :show => show_direct_product)
  inj = GrpAbFinGenMap[]
  pro = GrpAbFinGenMap[]
  j = 0
  for g = G
    if task in [:sum, :both]
      m = hom(g, Dp, GrpAbFinGenElem[Dp[j+i] for i = 1:ngens(g)])
      append!(L, m)
      push!(inj, m)
    end
    if task in [:prod, :both]
      m = hom(Dp, g, vcat(GrpAbFinGenElem[g[0] for i = 1:j], gens(g), GrpAbFinGenElem[g[0] for i=j+ngens(g)+1:ngens(Dp)]))
      append!(L, m)
      push!(pro, m)
    end
    j += ngens(g)
  end
  if task == :none
    return Dp
  elseif task == :prod
    return Dp, pro
  elseif task == :sum
    return Dp, inj
  else
    return Dp, pro, inj
  end
end

⊕(A::GrpAbFinGen...) = direct_product(A..., task = :none)
export ⊕

@doc Markdown.doc"""
    canonical_injection(G::GrpAbFinGen, i::Int) -> Map
Given a group $G$ that was created as a direct product, return the
injection from the $i$th component.
"""
function canonical_injection(G::GrpAbFinGen, i::Int)
  D = get_special(G, :direct_product)
  D === nothing && error("1st argument must be a direct product")
  s = sum(ngens(D[j]) for j = 1:i-1)
  h = hom(D[i], G, [G[s+j] for j = 1:ngens(D[i])])
  return h
end

@doc Markdown.doc"""
    canonical_projection(G::GrpAbFinGen, i::Int) -> Map
Given a group $G$ that was created as a direct product, return the
projection onto the $i$th component.
"""
function canonical_projection(G::GrpAbFinGen, i::Int)
  D = get_special(G, :direct_product)
  D === nothing && error("1st argument must be a direct product")
  H = D[i]
  h = hom(G, H, vcat( [GrpAbFinGenElem[H[0] for j = 1:ngens(D[h])] for h = 1:i-1]...,
                         gens(H),
                      [GrpAbFinGenElem[H[0] for j = 1:ngens(D[h])] for h = i+1:length(D)]...))
  return h
end

function matrix(M::Map{GrpAbFinGen, GrpAbFinGen})
  if typeof(M) == GrpAbFinGenMap
    return M.map
  end
  G = domain(M)
  return vcat([M(g).coeff for g = gens(G)])
end

function matrix(M::Generic.IdentityMap{GrpAbFinGen})
  return identity_matrix(FlintZZ, ngens(domain(M)))
end

@doc Markdown.doc"""
    hom(G::GrpAbFinGen, H::GrpAbFinGen, A::Array{ <: Map{GrpAbFinGen, GrpAbFinGen}, 2}) -> Map
Given groups $G$ and $H$ that are created as direct products as well
as a matrix $A$ containing maps $A[i,j] : G_i \to H_j$, return
the induced homomorphism.
"""
function hom(G::GrpAbFinGen, H::GrpAbFinGen, A::Array{ <: Map{GrpAbFinGen, GrpAbFinGen}, 2})
  r, c = size(A)
  if c == 1
    dG = [G]
  else
    dG = get_special(G, :direct_product)
  end
  if r == 1
    dH = [H]
  else
    dH = get_special(H, :direct_product)
  end
  if dG === nothing || dH === nothing
    error("both groups need to be direct products")
  end
  @assert all(i -> domain(A[i[1], i[2]]) == dG[i[1]] && codomain(A[i[1], i[2]]) == dH[i[2]], Base.Iterators.ProductIterator((1:r, 1:c)))
  h = hom(G, H, vcat([hcat([matrix(A[i,j]) for j=1:c]) for i=1:r]))
  return h
end

function _flat(G::GrpAbFinGen)
  s = get_special(G, :direct_product)
  if s === nothing
    return [G]
  end
  return vcat([_flat(x) for x = s]...)
end

function _tensor_flat(G::GrpAbFinGen)
  s = get_special(G, :tensor_product)
  if s === nothing
    return [G]
  end
  return vcat([_tensor_flat(x) for x = s]...)
end

@doc Markdown.doc"""
    flat(G::GrpAbFinGen) -> GrpAbFinGen, Map
Given a group $G$ that is created using (iterated) direct products, or
(iterated) tensor product,
return a group that is a flat product: $(A \oplus B) \oplus C$
is returned as $A \oplus B \oplus C$, (resp. $\otimes$)
together with the  isomorphism.
"""
function flat(G::GrpAbFinGen)
  s = get_special(G, :direct_product)
  if get_special(G, :direct_product) !== nothing
    H = direct_product(_flat(G)..., task = :none)
  elseif get_special(G, :tensor_product) !== nothing
    H = tensor_product(_tensor_flat(G)..., task = :none)
  else
    H = G
  end
  return hom(G, H, identity_matrix(FlintZZ, ngens(G)), identity_matrix(FlintZZ, ngens(G)))
end
######################################################################
# Lift of homomorphisms
######################################################################
#=
  G
  | phi
  V
  F <- H
    psi
 and Im(phi) subset Im(psi), then G -> H can be constructed
=#

@doc Markdown.doc"""
    lift(phi::Map, psi::Map) -> Map
Given $\phi: G\to F$ and $\psi:H \to F$ s.th. $\Im(\phi) \subseteq \Im(\psi)$
return the map $G\to H$ to make the diagram commute.
"""
function lift(phi::Map, psi::Map)
  x = [haspreimage(psi, image(phi, g)) for g = gens(domain(phi))]
  @assert all(t -> t[1], x)
  return hom(domain(phi), domain(psi), [t[2] for t = x])
end

@doc Markdown.doc"""
    zero_map(G::GrpAbFinGen) -> Map
Create the map $G \to \{0\}$.
"""
function zero_map(G::GrpAbFinGen)
  Z = AbelianGroup([1])
  set_name!(Z, "Zero")
  return hom(G, Z, [Z[0] for i=1:ngens(G)])
end

######################################################################
# complex/ free resolution
######################################################################
function iszero(h::T) where {T <: Map{<:GrpAbFinGen, <:GrpAbFinGen}}
  return all(x -> iszero(h(x)), gens(domain(h)))
end

mutable struct ChainComplex{T}
  @declare_other
  maps::Array{<:Map{<:T, <:T}, 1}
  direction::Symbol
  exact::Array{Bool, 1}
  function ChainComplex(A::S; check::Bool = true, direction:: Symbol = :left) where {S <:Array{<:Map{<:T, <:T}, 1}} where {T}
    if check
      @assert all(i-> iszero(A[i]*A[i+1]), 1:length(A)-1)
    end
    r = new{T}()
    r.maps = A
    r.direction = direction
    return r
  end
end

length(C::ChainComplex) = length(C.maps)
Base.map(C::ChainComplex, i::Int) = C.maps[i]
obj(C::ChainComplex, i::Int) = (i==0 ? domain(C.maps[1]) : codomain(C.maps[i]))

function show(io::IO, C::ChainComplex)
  @show_name(io, C)
  @show_special(io, C)

  Cn = get_special(C, :name)
  if Cn === nothing
    Cn = "C"
  end
  name_mod = String[]
  name_map = String[]
  mis_map = Tuple{Int, <:Map}[]
  mis_mod = Tuple{Int, <:Any}[]

  if C.direction == :left
    rng = 0:length(C)
    arr = ("--", "-->")
    dir = 1
  else
    rng = length(C)+1:-1:1
    arr = ("<--", "--")
    dir = 0
  end

  for i=1:length(C)
    phi = map(C, i)
    if get_special(phi, :name) !== nothing
      push!(name_map, get_special(phi, :name))
    else
      push!(name_map, "")
      push!(mis_map, (i, phi))
    end
  end
  for i=0:length(C)
    M = obj(C, i)
    if get_special(M, :name) !== nothing
      push!(name_mod, get_special(M, :name))
    else
      push!(name_mod, "$(Cn)_$i")
      push!(mis_mod, (i, M))
    end
  end

  io = IOContext(io, :compact => true)
  for i=rng
    if i == first(rng)
      print(io, name_mod[i+dir])
      continue
    end
    if name_map[i] != ""
      print(io, " ", arr[1], " ", name_map[i], " ", arr[2], " ", name_mod[i+dir])
    else
      print(io, " ", arr[1], arr[2], " ", name_mod[i+dir])
    end
  end
  if length(mis_mod) > 0 # || length(mis_map) > 0
    print(io, "\nwhere:\n")
    for (i, M) = mis_mod
      print(io, "\t$(Cn)_$i = ", M, "\n")
    end
#    for (i, phi) = mis_map
#      print(io, "\tphi_$i = ", phi, "\n")
#    end
  end
end

@doc Markdown.doc"""
    chain_complex(A::Map{GrpAbFinGen, GrpAbFinGen, <:Any, <:Any}...) -> ChainComplex{GrpAbFinGen}
Given maps $A_i$ s.th. $\Im(A_i) \subseteq \Kern(A_{i+1})$, this creates
the chain complex.
"""
function chain_complex(A::Map{GrpAbFinGen, GrpAbFinGen, <:Any, <:Any}...)
  return ChainComplex(collect(A))
end

function chain_complex(A::Array{<:Map{GrpAbFinGen, GrpAbFinGen, <:Any, <:Any}, 1})
  return ChainComplex(A)
end

Base.lastindex(C::ChainComplex) = length(C)
getindex(C::ChainComplex, u::UnitRange) = ChainComplex(C.maps[u], check = false)

@doc Markdown.doc"""
    isexact(C::ChainComplex) -> Bool
Tests is the complex $A_i: G_i \to G_{i+1}$
is exact, ie. if $\Im(A_i) = \Kern(A_{i+1})$.
"""
function isexact(C::ChainComplex)
  return all(i->iseq(image(C.maps[i])[1], kernel(C.maps[i+1])[1]), 1:length(C)-1)
end

@doc Markdown.doc"""
    free_resolution(G::GrpAbFinGen) -> ChainComplex{GrpAbFinGen}
A free resultion for $G$, ie. a chain complex terminating in
$G \to \{0\}$ that is exact.
"""
function free_resolution(G::GrpAbFinGen)
  A = DiagonalGroup(zeros(FlintZZ, ngens(G)))
  R = rels(G)
  B = DiagonalGroup(zeros(FlintZZ, nrows(R)))
  h_A_G = hom(A, G, gens(G))
  h_B_A = hom(B, A, [A(R[i, :]) for i=1:ngens(B)])
  Z = AbelianGroup(Int[1])
  set_name!(Z, "Zero")
  return chain_complex(hom(Z, B, [B[0]]), h_B_A, h_A_G, hom(G, Z, [Z[0] for i = 1:ngens(G)]))
end

mutable struct ChainComplexMap{T} <: Map{ChainComplex{T}, ChainComplex{T}, HeckeMap, ChainComplexMap}
  header::MapHeader{ChainComplex{T}, ChainComplex{T}}
  maps::Array{<:Map{<:T, <:T}, 1}
  function ChainComplexMap(C::ChainComplex{T}, D::ChainComplex{T}, A::S; check::Bool = !true) where {S <: Array{<:Map{<:T, <:T}, 1}} where {T}
    r = new{T}()
    r.header = MapHeader(C, D)
    r.maps = A
    return r
  end
end

@doc Markdown.doc"""
    hom(C::ChainComplex{T}, D::ChainComplex{T}, phi::Map{<:T, <:T}) where {T} -> ChainComplexMap
Given chain complexes $C_i: G_i \to G_{i+1}$ and $D_i: H_i \to H_{i+1}$
as well as a map $\phi = \phi_n: G_n \to H_n$, lift $\phi$ to
the entire complex: $\phi_i: G_i \to H_i$ s.th. all squares commute.
"""
function hom(C::ChainComplex{T}, D::ChainComplex{T}, phi::Map{<:T, <:T}) where {T}
  @assert length(C) == length(D)
  @assert domain(C.maps[end]) == domain(phi)
  @assert domain(D.maps[end]) == codomain(phi)

  h = [phi]
  for i=length(C)-1:-1:1
    push!(h, lift(C.maps[i]*h[end], D.maps[i]))
  end
  return ChainComplexMap(C, D, reverse(h))
end

@doc Markdown.doc"""
    hom(C::ChainComplex{T}, G::T) -> ChainComplex{T}
Given a complex $A_i: G_i \to G_{i+1}$ and a module $G$,
compute the derived complex $\hom(G_i, G)$.
"""
function hom(C::ChainComplex{GrpAbFinGen}, G::GrpAbFinGen)
  A = GrpAbFinGenMap[]
  H = [hom(domain(C.maps[1]), G)]
  H = vcat(H, [hom(codomain(f), G) for f = C.maps])

  R = GrpAbFinGenMap[]
  for i=1:length(C)
    A = H[i+1][1] # hom(C_i+1, G)
    B = H[i][1]   # hom(C_i  , G)
    #need map from A -> B
    #   C.maps[i] : E -> D
    D = codomain(C.maps[i])
    E = domain(C.maps[i])
    #  H[2][i+1]: A -> Hom(D, G)
    #  H[2][i]  : B -> hom(E, G)
    g = GrpAbFinGenElem[]
    for h = gens(A)
      phi = H[i+1][2](h) # D -> G
      psi = C.maps[i] * phi
      push!(g, preimage(H[i][2], psi))
    end
    push!(R, hom(A, B, g))
  end
  return ChainComplex(reverse(R))
end

@doc Markdown.doc"""
    hom(C::ChainComplex{T}, G::T) -> ChainComplex{T}
Given a complex $A_i: G_i \to G_{i+1}$ and a module $G$,
compute the derived complex $\hom(G, G_i)$.
"""
function hom(G::GrpAbFinGen, C::ChainComplex)
  A = GrpAbFinGenMap[]
  H = [hom(G, domain(C.maps[1]))]
  H = vcat(H, [hom(G, codomain(f)) for f = C.maps])

  R = GrpAbFinGenMap[]
  for i=1:length(C)
    A = H[i+1][1] # hom(G, C_i+1)
    B = H[i][1]   # hom(G, C_i)
    #need map from A -> B
    #   C.maps[i] : E -> D
    D = codomain(C.maps[i])
    E = domain(C.maps[i])
    #  H[2][i+1]: A -> Hom(G, D)
    #  H[2][i]  : B -> hom(G, E)
    g = GrpAbFinGenElem[]
    for h = gens(B)
      phi = H[i][2](h) # G -> E
      psi = phi * C.maps[i]
      push!(g, preimage(H[i+1][2], psi))
    end
    push!(R, hom(B, A, g))
  end
  return ChainComplex(R)
end

@doc Markdown.doc"""
    homology(C::ChainComplex{GrpAbFinGen}) -> Array{GrpAbFinGen, 1}
Given a complex $A_i: G_i \to G_{i+1}$,
compute the homology, ie. the modules $H_i = \Kern A_{i+1}/\Im A_i$
"""
function homology(C::ChainComplex{GrpAbFinGen})
  H = GrpAbFinGen[]
  for i=1:length(C)-1
    push!(H, snf(quo(kernel(C.maps[i+1])[1], image(C.maps[i])[1])[1])[1])
  end
  return H
end

function snake_lemma(C::ChainComplex{T}, D::ChainComplex{T}, A::Array{<:Map{T, T}, 1}) where {T}
  @assert length(C) == length(D) == 3
  @assert length(A) == 3
  @assert domain(A[1]) == obj(C,0) && codomain(A[1]) == obj(D, 1)
  @assert domain(A[2]) == obj(C,1) && codomain(A[2]) == obj(D, 2)
  @assert domain(A[3]) == obj(C,2) && codomain(A[3]) == obj(D, 3)

  ka, mka = kernel(A[1])
  kb, mkb = kernel(A[2])
  kc, mkc = kernel(A[3])
  ca, mca = cokernel(A[1])
  cb, mcb = cokernel(A[2])
  cc, mcc = cokernel(A[3])

  res = GrpAbFinGenMap[]
  push!(res, GrpAbFinGenMap(mka * map(C, 1) * inv(mkb)))
  push!(res, GrpAbFinGenMap(mkb * map(C, 2) * inv(mkc)))
  #now the snake
  push!(res, GrpAbFinGenMap(mkc * inv(map(C, 2)) * A[2] * inv(map(D, 2)) * mca))
  #and the boring rest
  push!(res, GrpAbFinGenMap(inv(mca) * map(D, 2) * mcb))
  push!(res, GrpAbFinGenMap(inv(mcb) * map(D, 3) * mcc))
  return chain_complex(res...)
end

################################################################################
#Tensor product
################################################################################

function tensor_product2(G::GrpAbFinGen, H::GrpAbFinGen)
  RG = rels(G)
  RH = rels(H)
  R = vcat(kronecker_product(RG', identity_matrix(FlintZZ, ngens(H)))',
           kronecker_product(identity_matrix(FlintZZ, ngens(G)), RH')')
  G = AbelianGroup(R)
end

struct TupleParent{T <: Tuple}
  function TupleParent(t::T) where {T}
    return new{T}()
  end
end

function show(io::IO, P::TupleParent{T}) where {T}
  print(io, "parent of tuples of type $T")
end

elem_type(::Type{TupleParent{T}}) where {T} = T
elem_type(::TupleParent{T}) where {T} = T

parent(t::Tuple) = TupleParent(t)

@doc Markdown.doc"""
    tensor_product(G::GrpAbFinGen...; task::Symbol = :map) -> GrpAbFinGen, Map
Given groups $G_i$ compute the tensor product $G_1\otimes \cdots \otimes G_n$.
If {{{task}}} is set to ":map", a map $\phi$ is returned that
maps tuples in $G_1 \times \cdots \times G_n$ to pure tensors
$g_1 \otimes \cdots \otimes g_n$. The map admits a preimage as well.
"""
function tensor_product(G::GrpAbFinGen...; task::Symbol = :map)
  @assert task in [:map, :none]
  local T
  if length(G) == 1
    T = G[1]
  else
    T = tensor_product2(G[2], G[1])
    for i = 3:length(G)
      T = tensor_product2(G[i], T)
    end
  end
  set_special(T, :tensor_product => G, :show => show_tensor_product)
  if task == :none
    return T
  end

  g = vec(collect(Base.Iterators.ProductIterator(Tuple(gens(g) for g = reverse(G)))))

  function pure(g::GrpAbFinGenElem...)
    @assert length(g) == length(G)
    @assert all(i-> parent(g[i]) == G[i], 1:length(G))

    return T(vec(collect(prod(x) for x = Base.Iterators.product([h.coeff for h = reverse(g)]...))))
  end
  function pure(T::Tuple)
    return pure(T...)
  end
  function inv_pure(t::GrpAbFinGenElem)
    p = Base.findall(i -> !iszero(t[i]), 1:ngens(T))
    if length(p) == 0
      return Tuple(collect(g[0] for g = G))
    end
    @assert length(p) == 1
    @assert t[p[1]] == 1
    return reverse(g[p[1]])
  end

  return T, MapFromFunc(pure, inv_pure, TupleParent(Tuple([g[0] for g = G])), T)
end

⊗(G::GrpAbFinGen...) = tensor_product(G..., task = :none)
export ⊗

@doc Markdown.doc"""
    hom(G::GrpAbFinGen, H::GrpAbFinGen, A::Array{ <: Map{GrpAbFinGen, GrpAbFinGen}, 1}) -> Map
Given groups $G = G_1 \otimes \cdots \otimes G_n$ and
$H = H_1 \otimes \cdot \otimes H_n$ as well as maps
$\phi_i: G_i\to H_i$, compute the tensor product of the maps.
"""
function hom(G::GrpAbFinGen, H::GrpAbFinGen, A::Array{ <: Map{GrpAbFinGen, GrpAbFinGen}, 1})
  tG = get_special(G, :tensor_product)
  tG === nothing && error("both groups must be tensor products")
  tH = get_special(H, :tensor_product)
  tH === nothing && error("both groups must be tensor products")
  @assert length(tG) == length(tH) == length(A)
  @assert all(i-> domain(A[i]) == tG[i] && codomain(A[i]) == tH[i], 1:length(A))
  M = matrix(A[1])'
  for i=2:length(A)
    M = kronecker_product(matrix(A[i])', M)
  end
  return hom(G, H, M')
end

@doc Markdown.doc"""
    tensor_product(C::ChainComplex{T}, G::T) -> ChainComplex{T}
Given a complex $A_i: G_i \to G_{i+1}$ and a module $G$,
compute the derived complex $G_i \otimes G$.
"""
function tensor_product(C::ChainComplex, G::GrpAbFinGen)
  A = GrpAbFinGenMap[]
  H = [tensor_product(domain(C.maps[1]), G, task  = :none)]
  H = vcat(H, [tensor_product(codomain(f), G, task = :none) for f = C.maps])

  R = GrpAbFinGenMap[]
  I = identity_map(G)
  for i = 1:length(C)
    push!(R, hom(H[i], H[i+1], [C.maps[i], I]))
  end
  return chain_complex(R)
end

################################################################################
#
#  Torsion
#
################################################################################

istorsion(G::GrpAbFinGen) = isfinite(G)

function torsion_subgroup(G::GrpAbFinGen)
  S, mS = snf(G)
  subs = GrpAbFinGenElem[]
  for i in 1:ngens(S)
    if !iszero(S.snf[i])
      push!(subs, mS(S[i]))
    end
  end
  return sub(G, subs)
end

##############################################################################
#
#  Subgroups
#
##############################################################################

@doc Markdown.doc"""
    sub(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1}) -> GrpAbFinGen, Map

Create the subgroup $H$ of $G$ generated by the elements in `s` together
with the injection $\iota : H \to G$.
"""
function sub(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)

  if length(s) == 0
    S = GrpAbFinGen(fmpz[1])
    I = zero_matrix(FlintZZ, ngens(S), ngens(G))
    mp = hom(S, G, I, check = false)
    if add_to_lattice
      append!(L, mp)
    end
    return S, mp
  end

  p = s[1].parent
  @assert G == p
  @assert all(x->x.parent == G, s)
  m = zero_matrix(FlintZZ, length(s) + nrels(p), ngens(p) + length(s))
  for i = 1:length(s)
    for j = 1:ngens(p)
      m[i + nrels(p), j] = s[i][j]
    end
    m[i + nrels(p), i + ngens(p)] = 1
  end
  if issnf(p)
    for i = 1:ngens(p)
      m[i, i] = p.snf[i]
    end
  else
    for i = 1:nrels(p)
      for j = 1:ngens(p)
        m[i, j] = p.rels[i, j]
      end
    end
  end
  if isdefined(G, :exponent)
    h = hnf_modular_eldiv(m, G.exponent)
  else
    h = hnf(m)
  end
  fstWithoutOldGens = 1
  for i in nrows(h):-1:1, j in ngens(p):-1:1
    if !iszero(h[i,j])
      fstWithoutOldGens = i + 1
      break
    end
  end
  r = view(h, fstWithoutOldGens:nrows(h), ngens(p) + 1:ncols(h))
  S = AbelianGroup(r)

  mS = hom(S, p, view(m, (nrels(p) + 1):nrows(h), 1:ngens(p)), check = false)

  if add_to_lattice
    append!(L, mS)
  end
  return S, mS
end

@doc Markdown.doc"""
    sub(s::Array{GrpAbFinGenElem, 1}) -> GrpAbFinGen, Map

Assuming that the non-empty array `s` contains elements of an abelian group
$G$, this functions returns the subgroup $H$ of $G$ generated by the elements
in `s` together with the injection $\iota : H \to G$.
"""
function sub(s::Array{GrpAbFinGenElem, 1},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  length(s) == 0 && error("Array must be non-empty")
  return sub(parent(s[1]), s, add_to_lattice, L)
end

@doc Markdown.doc"""
    sub(G::GrpAbFinGen, M::fmpz_mat) -> GrpAbFinGen, Map

Create the subgroup $H$ of $G$ generated by the elements corresponding to the
rows of $M$ together with the injection $\iota : H \to G$.
"""
function sub(G::GrpAbFinGen, M::fmpz_mat,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  m = zero_matrix(FlintZZ, nrows(M) + nrels(G), ngens(G) + nrows(M))
  for i = 1:nrows(M)
    for j = 1:ngens(G)
      m[i + nrels(G), j] = M[i,j]
    end
    m[i + nrels(G), i + ngens(G)] = 1
  end
  if issnf(G)
    for i = 1:ngens(G)
      m[i, i] = G.snf[i]
    end
  else
    for i = 1:nrels(G)
      for j = 1:ngens(G)
        m[i, j] = G.rels[i, j]
      end
    end
  end

  if isdefined(G, :exponent) && !iszero(G.exponent)
    h = hnf_modular_eldiv(m, G.exponent)
  else
    h = hnf(m)
  end
  fstWithoutOldGens = 1

  for i in nrows(h):-1:1, j in ngens(G):-1:1
    if !iszero(h[i,j])
      fstWithoutOldGens = i + 1
      break
    end
  end
  r = view(h, fstWithoutOldGens:nrows(h), ngens(G) + 1:ncols(h))
  S = AbelianGroup(r)
  mS = hom(S, G, view(m, (nrels(G) + 1):nrows(h), 1:ngens(G)), check = false)

  if add_to_lattice
    append!(L, mS)
  end
  return S, mS
end

function _sub_integer_snf(G::GrpAbFinGen, n::fmpz, add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  ind = 1
  while gcd(n, G.snf[ind]) == G.snf[ind] || ind <= ngens(G)
    ind += 1
  end
  if ind == ngens(G) && gcd(n, G.snf[ind]) == G.snf[ind]
    Gnew = GrpAbFinGenElem(Int[])
    mp = hom(Gnew, G, GrpAbFinGenElem[])
    if add_to_lattice
      append!(L, mp)
    end
    return Gnew, mp
  end
  invariants = Vector{fmpz}(undef, ngens(G)-ind+1)
  for_map = Vector{fmpz}(undef, ngens(G)-ind+1)
  for i = ind:length(G.snf)
    if iszero(G.snf[i])
      invariants[i-ind+1] = 0
    else
      res = gcd(n, V.snf[i])
      invariants[i-ind+1] = divexact(G.snf[i], res)
    end
  end
  Gnew = DiagonalGroup(invariants)
  mat_map = zero_matrix(FlintZZ, length(invariants), ngens(G))
  for i = 1:ngens(Gnew)
    mat_map[i, ind+i-1] = n
  end
  mp = hom(Gnew, G, mat_map)
  if add_to_lattice
    append!(L, mp)
  end
  return Gnew, mp
end

@doc Markdown.doc"""
    sub(G::GrpAbFinGen, n::fmpz) -> GrpAbFinGen, Map

Create the subgroup $n \cdot G$ of $G$ together
with the injection $\iota : n\cdot G \to G$.
"""
function sub(G::GrpAbFinGen, n::fmpz,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  if issnf(G)
    return _sub_integer_snf(G, n, add_to_lattice, L)
  end
  H, mH = sub(G, scalar_matrix(FlintZZ, ngens(G), deepcopy(n)), add_to_lattice, L)
  if isdefined(G, :exponent)
    res = divexact(G.exponent, gcd(n, G.exponent))
    if !iszero(res)
      H.exponent = res
    end
  end
  return H, mH
end

@doc Markdown.doc"""
    sub(G::GrpAbFinGen, n::Integer) -> GrpAbFinGen, Map

Create the subgroup $n \cdot G$ of $G$ together
with the injection $\iota : n \cdot G \to G$.
"""
function sub(G::GrpAbFinGen, n::Integer,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  return sub(G, fmpz(n), add_to_lattice, L)
end

################################################################################
#
#  Quotients
#
################################################################################

@doc Markdown.doc"""
  quo(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1}) -> GrpAbFinGen, Map

Create the quotient $H$ of $G$ by the subgroup generated by the elements in
$s$, together with the projection $p : G \to H$.
"""
function quo(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  if length(s) == 0
    I = identity_matrix(FlintZZ, ngens(G))
    m = hom(G, G, I, I, check = false)
    if add_to_lattice
      append!(L, m)
    end
    return G, m
  end

  p = s[1].parent
  @assert G == p
  m = zero_matrix(FlintZZ, length(s)+nrels(p), ngens(p))
  for i = 1:length(s)
    for j = 1:ngens(p)
      m[i + nrels(p),j] = s[i][j]
    end
  end
  if issnf(p)
    for i = 1:ngens(p)
      m[i, i] = p.snf[i]
    end
  else
    for i = 1:nrels(p)
      for j = 1:ngens(p)
        m[i, j] = p.rels[i, j]
      end
    end
  end

  Q = AbelianGroup(m)
  I = identity_matrix(FlintZZ, ngens(p))
  m = hom(p, Q, I, I, check = false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

@doc Markdown.doc"""
  quo(G::GrpAbFinGen, M::fmpz_mat) -> GrpAbFinGen, Map

Create the quotient $H$ of $G$ by the subgroup generated by the elements
corresponding to the rows of $M$. together with the projection $p : G \to H$.
"""
function quo(G::GrpAbFinGen, M::fmpz_mat,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  m = vcat(rels(G), M)
  Q = AbelianGroup(m)
  I = identity_matrix(FlintZZ, ngens(G))
  m = hom(G, Q, I, I, check = false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

@doc Markdown.doc"""
    quo(G::GrpAbFinGen, n::Integer}) -> GrpAbFinGen, Map
    quo(G::GrpAbFinGen, n::fmpz}) -> GrpAbFinGen, Map

Returns the quotient $H = G/nG$ together with the projection $p : G \to H$.
"""
function quo(G::GrpAbFinGen, n::Union{fmpz, Integer},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  if issnf(G)
    return quo_snf(G, n, add_to_lattice, L)
  else
    return quo_gen(G, n, add_to_lattice, L)
  end
end

function quo_snf(G::GrpAbFinGen, n::Union{fmpz, Integer},
                 add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  r = [gcd(x, n) for x = G.snf]
  I = identity_matrix(FlintZZ, ngens(G))
  Q = DiagonalGroup(r)
  m = hom(G, Q, I, I, check = false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

function quo_gen(G::GrpAbFinGen, n::Union{fmpz, Integer},
                 add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  m = vcat(G.rels, n*identity_matrix(FlintZZ, ngens(G)))
  Q = AbelianGroup(m)
  if isdefined(G, :exponent)
    Q.exponent = gcd(n, G.exponent)
  end
  I = identity_matrix(FlintZZ, ngens(G))
  m = hom(G, Q, I, I, check = false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

################################################################################
#
#  use lattice...
#
################################################################################

function +(G::GrpAbFinGen, H::GrpAbFinGen, L::GrpAbLattice = GroupLattice)
  fl, GH, mG, mH = can_map_into_overstructure(L, G, H)
  if !fl
    error("no common overgroup known")
  end
  return sub(GH, vcat([GH(mG[i, :]) for i=1:nrows(mG)], [GH(mH[i, :]) for i=1:nrows(mH)]))[1]
end

function Base.intersect(G::GrpAbFinGen, H::GrpAbFinGen, L::GrpAbLattice = GroupLattice)
  fl, GH, mG, mH = can_map_into_overstructure(L, G, H)
  if !fl
    error("no common overgroup known")
  end
  #M = [ mG identity_matrix(FlintZZ, nrows(mG)); mH zero_matrix(FlintZZ, nrows(mH), nrows(mG)) ;
        #rels(GH) zero_matrix(FlintZZ, nrels(GH), nrows(mG))]
  M = vcat(vcat(hcat(mG, identity_matrix(FlintZZ, nrows(mG))), hcat(mH, zero_matrix(FlintZZ, nrows(mH), nrows(mG)))), hcat(rels(GH), zero_matrix(FlintZZ, nrels(GH), nrows(mG))))
  h = hnf(M)
  i = nrows(h)
  while i > 0 && iszero(sub(h, i:i, 1:ngens(GH)))
    i -= 1
  end
  return sub(G, [G(sub(h, j:j, ngens(GH)+1:ncols(h))) for j=i+1:nrows(h)])[1]
end

function Base.intersect(A::Array{GrpAbFinGen, 1})
  a = first(A)
  for b = A
    a = intersect(a, b, true)
  end
  return a
end

function Base.issubset(G::GrpAbFinGen, H::GrpAbFinGen, L::GrpAbLattice = GroupLattice)
  fl, GH, mG, mH = can_map_into_overstructure(L, G, H)
  if !fl
    error("no common overgroup known")
  end
  hH = hom(H, GH, mH)
  return all(x -> haspreimage(hH, GH(mG[x, :]))[1], 1:nrows(mG))
end

function issubgroup(G::GrpAbFinGen, H::GrpAbFinGen, L::GrpAbLattice = GroupLattice)
  fl, GH, mG, mH = can_map_into_overstructure(L, G, H)
  if !fl
    error("no common overgroup known")
  end
  hH = hom(H, GH, mH)
  n = matrix(FlintZZ, 0, ngens(H), fmpz[])
  for j=1:nrows(mG)
    fl, x = haspreimage(hH, GH(mG[j, :]))
    if !fl
      return false, hH
    end
    n = vcat(n, x.coeff)
  end
  return true, hom(G, H, n)
end

#cannot define == as this produces problems elsewhere... need some thought
function iseq(G::GrpAbFinGen, H::GrpAbFinGen, L::GrpAbLattice = GroupLattice)
  isfinite(G) && (order(G) == order(H) || return false)
  return issubgroup(G, H)[1] && issubgroup(H, G)[1]
end

function Base.isequal(G::GrpAbFinGen, H::GrpAbFinGen)
  return G === H
end

function quo(G::GrpAbFinGen, U::GrpAbFinGen)
  fl, m = issubgroup(U, G)
  fl || error("not a subgroup")
  return quo(G, m.map)
end
################################################################################
#
#  Make Smith normal form
#
################################################################################

function make_domain_snf(m::Map{GrpAbFinGen, T}) where T
  G = domain(m)
  S, mS = snf(G)
  return mS*m
end

################################################################################
#
#  Cyclic
#
################################################################################

@doc Markdown.doc"""
    iscyclic(G::GrpAbFinGen) -> Bool

Returns whether $G$ is cyclic.
"""
function iscyclic(G::GrpAbFinGen)
  if !issnf(G)
    S = snf(G)[1]
    return ngens(S) == 1
  else
    return ngens(G) == 1
  end
end

################################################################################
#
#  p-Sylow subgroup
#
################################################################################

function _psylow_subgroup_gens(G::GrpAbFinGen, p::Union{fmpz, Integer})
  @assert issnf(G)
  z = GrpAbFinGenElem[]
  for i in 1:ngens(G)
    k, m = remove(G.snf[i], p)
    if k != 0
      push!(z, m*G[i])
    end
  end
  return z
end

function psylow_subgroup(G::GrpAbFinGen, p::Union{fmpz, Integer},
                         to_lattice::Bool = true)
  S, mS = snf(G)
  z = _psylow_subgroup_gens(S, p)
  zz = [ image(mS, x) for x in z ]
  return sub(G, zz, to_lattice)
end

################################################################################
#
#  Some special abelian groups
#
################################################################################

# TH: Isn't this the same as UnitsModM.jl?
# TODO: remove this from here. It does not belong here
@doc Markdown.doc"""
    multgrp_of_cyclic_grp(n::fmpz) -> GrpAbFinGen

Returns the multiplicative group of the cyclic group with $n$ elements.
"""
function multgrp_of_cyclic_grp(n::fmpz)
  composition = Array{fmpz,1}()
  for (p,mp) in factor(n)
    if (p == 2) && (mp >= 2)
      push!(composition,2)
      push!(composition,fmpz(2)^(mp-2))
    else
      push!(composition,(p-1)*p^(mp-1))
    end
  end
  return DiagonalGroup(composition)
end

multgrp_of_cyclic_grp(n::Integer) = multgrp_of_cyclic_grp(fmpz(n))

################################################################################
#
#  Isomorphism to abelian groups
#
################################################################################

@doc Markdown.doc"""
    find_isomorphism(G, op, A::GrpAb) -> Dict, Dict

Given an abelian group $A$ and a collection $G$ which is an abelian group with
the operation `op`, this functions returns isomorphisms $G \to A$ and $A \to G$
encoded as dictionaries.

It is assumed that $G$ and $A$ are isomorphic.
"""
function find_isomorphism(G, op, A::GrpAb)
  H, GtoH, HtoG = find_isomorphism_with_abelian_group(G, op)
  @assert issnf(H)
  Asnf, AsnftoA = snf(A)
  s = length(H.snf)
  id = identity_matrix(FlintZZ, s)
  AsnftoH = hom(Asnf, H, id, id)
  GtoA = Dict{eltype(G), GrpAbFinGenElem}()
  AtoG = Dict{GrpAbFinGenElem, eltype(G)}()
  for g in G
    GtoA[g] = AsnftoA(AsnftoH\(GtoH[g]))
  end
  for a in A
    AtoG[a] = HtoG[AsnftoH(AsnftoA\a)]
  end

  #for g in G
  #  @assert AtoG[GtoA[g]] == g
  #end
  #for a in G
  #  @assert AtoG[GtoA[a]] == a
  #end
  return GtoA, AtoG
end

function find_isomorphism_with_abelian_group(G, op)
  id = find_identity(G, op)
  S = small_generating_set(G, op, id)
  list = push!(copy(S), id)
  n = length(G)
  elem_to_index = Dict{eltype(G), Int}()
  words = Dict{eltype(G), Array{Int}}()
  rels = Vector{Vector{Int}}()

  l = length(list)

  for i in 1:length(S)
    v = zeros(Int, length(S))
    v[i] = 1
    words[S[i]] = v
  end

  words[list[end]] = zeros(Int, length(S))

  for i in 1:length(G)
    elem_to_index[G[i]] = i
  end

  first_round = true # One has to do this once even if length(list) == n
  while first_round || length(list) != n
    first_round = false
    for g in list
      for i in 1:length(S)
        s = S[i]
        m = op(g, s)

        if m in list
          w = words[S[i]] .+ words[g] .- words[m]
          if !iszero(w)
            push!(rels, w)
          end
        end

        if !(m in list)
          push!(list, m)
          words[m] = words[s] .+ words[g]
        end
      end
    end
  end

  rel_mat = zero_matrix(FlintZZ, length(rels), length(S))
  for i in 1:length(rels)
    for j in 1:length(S)
      rel_mat[i, j] = rels[i][j]
    end
  end

  A = AbelianGroup(rel_mat)
  Asnf, mAsnf = snf(A)
  GtoAsnf = Dict{eltype(G), GrpAbFinGenElem}()
  AsnftoG = Dict{GrpAbFinGenElem, eltype(G)}()
  for g in G
    w = words[g]
    GtoAsnf[g] = sum(w[i] * (mAsnf \ A[i]) for i in 1:length(S))
  end
  for a in Asnf
    x = id
    v = Int[mAsnf(a).coeff[1, i] for i in 1:length(S)]
    for i in 1:length(S)
      x = op(x, pow(S[i], op, v[i], id))
    end
    AsnftoG[a] = x
  end

  for g in G
    @assert AsnftoG[GtoAsnf[g]] == g
  end
  for a in Asnf
    @assert GtoAsnf[AsnftoG[a]] == a
  end

  return Asnf, GtoAsnf, AsnftoG
end

function pow(x, op, n, id)
  if n == 0
    return id
  end
  y = x
  for i in 1:(n-1)
    y = op(y, x)
  end
  return y
end

################################################################################
#
#  All abelian groups
#
################################################################################

@doc Markdown.doc"""
    abelian_groups(n::Int) -> Vector{GrpAbFinGen}

Given a positive integer $n$, return a list of all abelian groups of order $n$.
"""
function abelian_groups(n::Int)
  if n == 1
    return GrpAbFinGen[DiagonalGroup(Int[])]
  end
  nn = fmpz(n)
  fac = factor(nn)
  sylow_lists = Vector{Vector{GrpAbFinGen}}()
  for (p, e) in fac
    push!(sylow_lists, GrpAbFinGen[DiagonalGroup(fmpz[p^i for i in reverse(t)]) for t in AllParts(e)])
  end
  C = Base.Iterators.product(sylow_lists...)
  grps = GrpAbFinGen[]
  for c in C
    G = c[1]
    for i in 2:length(fac)
      G = snf(direct_product(G, c[i], task = :none))[1]
    end
    push!(grps, G)
  end
  return grps
end

################################################################################
#
#  Sum of maps
#
################################################################################

function +(f::GrpAbFinGenMap, g::GrpAbFinGenMap)
  @assert domain(f) == domain(g)
  @assert codomain(f) == codomain(g)
  return hom(domain(f), codomain(g), GrpAbFinGenElem[f(x)+g(x) for x in gens(domain(f))])
end

function -(f::GrpAbFinGenMap, g::GrpAbFinGenMap)
  @assert domain(f) == domain(g)
  @assert codomain(f) == codomain(g)
  return hom(domain(f), codomain(g), GrpAbFinGenElem[f(x)-g(x) for x in gens(domain(f))])
end

################################################################################
#
#  Functions related to the action of a group of endomorphisms
#
################################################################################

function induce_action_on_subgroup(mS::GrpAbFinGenMap, acts::Vector{GrpAbFinGenMap})
  res = Vector{GrpAbFinGenMap}(undef, length(acts))
  S = domain(mS)
  for i = 1:length(acts)
    imgs = Vector{GrpAbFinGenElem}(undef, ngens(S))
    for j = 1:length(imgs)
      imgs[j] = mS\(acts[i](mS(S[j])))
    end
    res[i] = hom(S, S, imgs)
  end
  return res
end

function fixed_subgroup(f::GrpAbFinGenMap, to_lattice::Bool = true)
  @assert domain(f) == codomain(f)
  return kernel(f - id_hom(domain(f)), to_lattice)
end

function isfixed_point_free(act::Vector{GrpAbFinGenMap})
  G = domain(act[1])
  intersection_of_kernels = G
  minus_id = hom(G, G, GrpAbFinGenElem[-x for x in gens(G)])
  for i = 1:length(act)
    k, mk = fixed_subgroup(act[i])
    intersection_of_kernels = intersect(intersection_of_kernels, k)
    if order(intersection_of_kernels) == 1
      return true
    end
  end
  return false
end

################################################################################
#
#  Elementary divisors
#
################################################################################

@doc Markdown.doc"""
    elementary_divisors(G::GrpAbFinGen) -> Vector{fmpz}

Given $G$, returns the elementary divisors of $G$, that is, the unique positive
integers $e_1,\dotsc,e_k$ with $e_i \mid e_{i + 1}$ and
$G \cong \mathbf{Z}/e_1\mathbf{Z} \times \dotsb \times \mathbf{Z}/e_k\mathbf{Z}$.
"""
function elementary_divisors(G::GrpAbFinGen)
  if issnf(G)
    return copy(G.snf)
  else
    return elementary_divisors(snf(G)[1])
  end
end

################################################################################
#
#  Has quotient
#
################################################################################

@doc Markdown.doc"""
    has_quotient(G::GrpAbFinGen, invariant::Vector{Int}) -> Bool

Given an abelian group $G$, returns true if it has a quotient with given elementary
divisors and false otherwise.
"""
function has_quotient(G::GrpAbFinGen, invariants::Array{Int,1})
  H = DiagonalGroup(invariants)
  H = snf(H)[1]
  G1 = snf(G)[1]
  arr_snfG1 = filter(x -> x != 1, G1.snf)
  arr_snfH = filter(x -> x != 1, H.snf)
  if length(arr_snfG1) < length(arr_snfH)
    return false
  end
  for i = 0:length(arr_snfH)-1
    if !divisible(arr_snfG1[end-i], arr_snfH[end-i])
      return false
    end
  end
  return true
end

################################################################################
#
#  Find complement
#
################################################################################

#TODO: a better algorithm?
@doc Markdown.doc"""
    has_complement(f::GrpAbFinGenMap) -> Bool, GrpAbFinGenMap

Given a map representing a subgroup of a group $G$, returns either true and
an injection of a complement in $G$, or false.
"""
function has_complement(m::GrpAbFinGenMap)
  G = codomain(m)
  S, mS = snf(G)
  H, mH = image(m*inv(mS))
  els = GrpAbFinGenElem[]
  for i = 1:ngens(S)
    new_els = vcat(els, GrpAbFinGenElem[S[i]])
    new_sub, map_new_sub = sub(S, new_els)
    if isone(order(intersect(new_sub, H)))
      els = new_els
    end
  end
  ss, mss = sub(S, els)
  if order(ss)*order(H) != order(S)
    return false, mss
  else
    els_G = GrpAbFinGenElem[mS(x) for x in els]
    return true, sub(G, els_G)[2]
  end
end

################################################################################
#
#  Identity
#
################################################################################

id(G::GrpAbFinGen) = G(zeros(fmpz, ngens(G)))


text:
struct Events{T <: DiseaseStateSequence}
  exposure::Union{Nothing, Vector{Float64}}
  infection::Vector{Float64}
  removal::Union{Nothing, Vector{Float64}}

  function Events{T}(e::V, i::V, r::V) where {T <: SEIR, V <: Vector{Float64}}
    if length(unique((length.([e; i; r])))) != 1
      @error "Length of event time vectors must be equal"
    end
    return new{T}(e, i, r)
  end

  function Events{T}(n::Integer) where T <: SEIR
    return new{T}(fill(NaN, n), fill(NaN, n), fill(NaN, n))
  end

  function Events{T}(e::V, i::V) where {T <: SEI, V <: Vector{Float64}}
    if length(unique((length.([e; i])))) != 1
      @error "Length of event time vectors must be equal"
    end
    return new{T}(e, i, nothing)
  end

  function Events{T}(n::Integer) where T <: SEI
    return new{T}(fill(NaN, n), fill(NaN, n), nothing)
  end

  function Events{T}(i::V, r::V) where {T <: SIR, V <: Vector{Float64}}
    if length(unique((length.([i; r])))) != 1
      @error "Length of event time vectors must be equal"
    end
    return new{T}(nothing, i, r)
  end

  function Events{T}(n::Integer) where T <: SIR
    return new{T}(nothing, fill(NaN, n), fill(NaN, n))
  end

  function Events{T}(i::V) where {T <: SI, V <: Vector{Float64}}
    return new{T}(nothing, i, nothing)
  end

  function Events{T}(n::Integer) where T <: SI
    return new{T}(nothing, fill(NaN, n), nothing)
  end
end

function Events{T}(a::Array{Float64,2}) where T <: DiseaseStateSequence
  if size(a, 2) == 3
    return Events{T}(a[:,1], a[:,2], a[:,3])
  elseif size(a, 2) == 2
    return Events{T}(a[:,1], a[:,2])
  elseif size(a, 2) == 1
    return Events{T}(a[:,1])
  else
    @error "Invalid array size for construction of an $(Events{T}) object"
  end
end

function Events{T}(x::DiseaseStates) where T <: DiseaseStateSequence
  events = Events{T}(length(x))
  for i = 1:length(x)
    if x[i] in convert(DiseaseStates, T)
      for j = convert(DiseaseStates, T)[2:findfirst(Ref(x[i]) .== convert(DiseaseStates, T))]
        events[j][i] = -Inf
      end
    else
      @error "Invalid initial state for individual $i"
    end
  end
  return events
end

function individuals(x::Events{M}) where{
                     M <: DiseaseStateSequence}
  return length(x.infection)
end

function Base.show(io::IO, x::Events{T}) where T <: DiseaseStateSequence
  return print(io, "$T model event times (n=$(individuals(x)))")
end

function Base.getindex(x::Events{T}, new_state::DiseaseState) where T <: DiseaseStateSequence
  if new_state == State_E
    return x.exposure
  elseif new_state == State_I
    return x.infection
  elseif new_state == State_R
    return x.removal
  else
    error("Unsupported indexing disease state")
  end
end

function Base.getindex(x::Events{T}, states::DiseaseStates) where T <: DiseaseStateSequence
  y = x[states[1]]
  for i = 2:length(states)
    y = hcat(y, x[states[i]])
  end
  return y
end

function Base.convert(::Type{Array{Float64, 2}}, x::Events{T}) where T <: DiseaseStateSequence
  return x[convert(DiseaseStates, T)[2:end]]
end

function Base.convert(::Type{Vector{Float64}}, x::Events{T}) where T <: DiseaseStateSequence
  return x[convert(DiseaseStates, T)[2:end]][:]
end

function Base.convert(::Type{Array{Float64, 2}}, x::Array{Events{T}, 1}) where T <: DiseaseStateSequence
  y = convert(Vector{Float64}, x[1])'
  for i = 2:length(x)
    y = vcat(y, convert(Vector{Float64}, x[i])')
  end
  return y
end

function Base.minimum(x::Events{T}) where T <: DiseaseStateSequence
  y = convert(Array{Float64, 2}, x)
  return minimum(y[(y .!== NaN) .& (y .> -Inf)])
end

function Base.maximum(x::Events{T}) where T <: DiseaseStateSequence
  y = convert(Array{Float64, 2}, x)
  return maximum(y[(y .!== NaN) .& (y .> -Inf)])
end

function mean(x::Vector{Events{T}}) where T <: DiseaseStateSequence
  return Events{T}(mean([convert(Array{Float64, 2}, i) for i in x]))
end

ext: jl


text:
using Game2048
using BenchmarkTools
using StatsBase

board = initboard()
move!(board, left)

@time res = simulate_games(10_000);
res
@time mean(res)
@time extrema(res)

@benchmark simulate_game()


text:
using .lib

module execute
	export exec
	import ..Data,..clean,..exitWithError,..MMNone,..MMSerial,..MMSpawn,..MMThread
	using Printf,Base.Threads

	d=nothing
	o=nothing
	e=nothing
	r=nothing

	ec=0
	res=""

	function exec(rd::Data)
		global d,o,e,r
		d=rd
		o=co2f(d.out,stdout)
		e=co2f(d.err,stderr)
		r=ro2f(d.result)

		if d.multiple==MMNone single()
		elseif d.multiple==MMSerial serial()
		elseif d.multiple==MMSpawn spawn()
		elseif d.multiple==MMThread thread()
		end

		print(r,res)
		close(o)
		close(e)
		close(r)
		exit(ec)
	end

	function single()
		global res

		p=SP(`$(d.command)`,join(d.command," "))

		sec=@elapsed begin
			runSP(p)
		end

		res=clean("""
			time: $(descTime(sec))
			process id: $(p.pid)
			$(descEC(p))
		""")
		ec=p.ec
	end

	function serial()
		global res

		pl=multiple()
		lp=pl[length(pl)]

		sec=@elapsed begin
			for p in pl
				runSP(p)
				if p.ec!=0
					lp=p
					break
				end
			end
		end

		res=join([
			"time: $(descTime(sec))",
			map(p->"process$(p.order) id: $(p.pid<0 ? "N/A" : p.pid)",pl)...,
			descEC(lp),""
		],"\n")
		ec=lp.ec
	end

	function spawn()
		pl=multiple()

		sec=@elapsed begin
			for p in pl startSP(p) end
			for p in pl waitSP(p) end
		end

		collect(pl,sec)
	end

	function thread()
		pl=multiple()

		sec=@elapsed begin
			@threads for p in pl
				runSP(p)
			end
		end

		collect(pl,sec)
	end

	Base.@kwdef mutable struct SP
		cmd::Base.CmdRedirect
		description::String
		process=nothing
		order::Int = 0
		pid::Int32 = -1
		ec::Int = 0
		SP(cmdObj,desc) = new(pipeline(cmdObj,stdout=stdout,stderr=stderr,stdin=stdin),desc)
	end
	function startSP(p::SP)
		try
			p.process=run(p.cmd,wait=false)
			pid=getpid(p.process)
			p.pid=pid
		catch e
			exitWithError("実行に失敗しました: $(p.description)")
		end
	end
	function waitSP(p::SP)
		wait(p.process)
		p.ec=p.process.exitcode
	end
	function runSP(p::SP)
		startSP(p)
		waitSP(p)
	end
	function descEC(p::SP)::String
		return @sprintf("exit code: %.d",p.ec)
	end
	function multiple()::Array{SP,1}
		pl=fill(-1,length(d.command))
		n=1
		l::Array{SP,1}=[]
		for c in d.command
			if occursin(r"#\{\}\(\)\[\]<>\|\&\*\?~;",c)
				sh=get(ENV,"SHELL","sh")
				p=SP(`$sh -c $c`)
			else
				s=c
				s=replace(s,"\\"=>"\\\\")
				s=replace(s,"`"=>"\\`")
				s=replace(s,"\$"=>"\\\$")
				try
					p=SP(eval(Meta.parse("`$s`")),d.command)
				catch e
					exitWithError("パースに失敗しました: $c")
				end
			end
			p.order=n
			n+=1
			append!(l,[p])
		end
		return l
	end
	function collect(pl::Array{SP,1},sec)
		global res,ec

		res=join(Iterators.flatten([
			["time $(descTime(sec))"],
			map(function (p)
			if p.ec>ec ec=p.ec end
				return ["process$(p.order) id: $(p.pid)",descEC(p)]
			end,pl)...,[""]
		]),"\n")
	end

	function co2f(d,i)
		if d=="inherit" return i
		elseif d=="discard" return devnull
		else return fh(d)
		end
	end

	function ro2f(d)
		if d=="stdout" return stdout
		elseif d=="stderr" return stderr
		else return fh(d)
		end
	end

	opened=Dict()
	function fh(path)
		if haskey(opened,path) return opened[path]
		else
			try
				io=open(path,"a")
				opened[path]=io
				return io
			catch e
				exitWithError("指定したパスには書き込みできません: "+path)
			end
		end
	end

	function descTime(sec)::String
		t=""
		r=sec/3600
		v=floor(r)
		if v>=1 t*=@sprintf("%.0fh ",v) end
		r=(r-v)*60
		v=floor(r)
		if v>=1 t*=@sprintf("%.0fm ",v) end
		r=(r-v)*60
		v=floor(r)
		if v>=1 t*=@sprintf("%.0fs ",v) end
		r=(r-v)*1000
		t*=@sprintf("%07.3fms",r)
		return t
	end

end

max_stars_repo_path: Time Measurement/Julia/execute.jl

max_issues_repo_path: Time Measurement/Julia/execute.jl

max_forks_repo_path: Time Measurement/Julia/execute.jl


text:
#=
This code is loosely based on
Strid and Walentin (2008) "Block Kalman filtering for large-scale DSGE models"
=#

mutable struct CTBlockKalmanFilter{S<:AbstractFloat}
    T_powers::Dict{Int64, Matrix{S}}
    R_blocks::Dict{Int64, Matrix{S}}
    C::Vector{S}
    Q::Matrix{S}
    Z::Matrix{S}
    D::Vector{S}
    E::Matrix{S}
    RQRp::Dict{Tuple{Int64, Int64}, Matrix{S}} # RQR', indexed block-wise
    s_t::Vector{S}
    P_t::Dict{Tuple{Int64, Int64}, Matrix{S}}  # P_t indexed block-wise
    n_simulate_states::Int64
    loglh_t::S     # P(y_t | y_{1:t})
end

"""
```
KalmanFilter(T, R, C, Q, Z, D, E, [s_0, P_0])
```

Outer constructor for the `KalmanFilter` type.
"""
function CTBlockKalmanFilter(T::Matrix{S}, R::Matrix{S}, C::Vector{S}, Q::Matrix{S},
                             Z::Matrix{S}, D::Vector{S}, E::Matrix{S}, s_0::Vector{S} = Vector{S}(0),
                             P_0::Matrix{S} = Matrix{S}(0, 0),
                             n_simulate_states::Int64=1) where {S<:AbstractFloat}
    if isempty(s_0) || isempty(P_0)
        s_0, P_0 = init_stationary_states(T, R, C, Q)
    end

    # Create dictionaries of products of T matrices, plus create blocks of R matrices/extend C -> TT, RR, CC
    T_powers = Dict{Int64, Matrix{S}}()  # key is power of T
    R_blocks = Dict{Int64, Matrix{S}}() # key is power of T multiplying R
    R_blocks[0] = R
    T_powers[1] = T
    for i = 2:n_simulate_states
        T_powers[i] = T_powers[i - 1] * T
        R_blocks[i - 1] = T_powers[i - 1] * R
    end

    # Create RQR' dictionary
    RQRp[(1,1)] = R_blocks[0] * Q * R_blocks[0]' # initialize first row as base case
    if n_simulate_states > 1
        for i = 2:n_simulate_states
            RQRp[(1,i)] = R-blocks[0] * Q * R_blocks[i - 1]'
        end

        # Recursive construction of rest of RQR' dictionary
        for i = 2:n_simulate_states
            for j = 1:n_simulate_states
                if i <= j
                    RQRp[(i, j)] = R_blocks[i - 1] * Q * R_blocks[j - 1] + RQRp[(i-1,j-1)]
                else
                    # RQR' is symmetric
                    RQRp[(i, j)] = RQRp[(j, i)]
                end
            end
        end
    end

    # Create P_t dictionary
    P_dict = Dict{Tuple{Int64, Int64}, Matrix{S}}()
    T_size = size(T, 1)
    P_dict[(1, 1)] = P_0[1:T_size, 1:T_size]
    if n_simulate_states > 1
        for i = 2:n_simulate_states
            for j = i:n_simulate_states
                P_dict[(i, j)] = P_0[1 + T_size * (i-1):T_size * i, 1 + T_size * (j-1):T_size * j]
            end
        end
    end

    return CTBlockKalmanFilter(T_powers, R_blocks, C, Q, Z, D, E, RQRp,
                               s_0, P_dict, n_simulate_states, NaN)
end

"""
```
ct_block_kalman_filter(y, T, R, C, Q, Z, D, E, s_0 = Vector(), P_0 = Matrix();
    outputs = [:loglh, :pred, :filt], Nt0 = 0)
```
This function implements the Kalman filter for the following state-space model:

```
s_{t+1} = C + T*s_t + R*ϵ_t    (transition equation)
y_t     = D + Z*s_t + u_t      (measurement equation)

ϵ_t ∼ N(0, Q)
u_t ∼ N(0, E)
Cov(ϵ_t, u_t) = 0
```

### Inputs

- `y`: `Ny` x `Nt` matrix containing data `y_1, ... , y_T`
- `s_0`: optional `Ns` x 1 initial state vector
- `P_0`: optional `Ns` x `Ns` initial state covariance matrix
- `T`: `Ns` x `Ns` state transition matrix
- `R`: `Ns` x `Ne` matrix in the transition equation mapping shocks to states
- `C`: `Ns` x 1 constant vector in the transition equation
- `Q`: `Ne` x `Ne` matrix of shock covariances
- `Z`: `Ny` x `Ns` matrix in the measurement equation mapping states to
  observables
- `D`: `Ny` x 1 constant vector in the measurement equation
- `E`: `Ny` x `Ny` matrix of measurement error covariances

where:

- `Nt`: number of time periods for which we have data
- `Ns`: number of states
- `Ne`: number of shocks
- `Ny`: number of observables

### Keyword Arguments

- `outputs`: some subset of `[:loglh, :pred, :filt]` specifying which outputs to
  compute and return. There will always be the same number of return values,
  but, for example, `s_pred` and `P_pred` will be returned as empty arrays if
  `:pred` is not in `outputs`
- `Nt0`: number of presample periods to omit from all return values

### Outputs

- `loglh`: length `Nt` vector of conditional log-likelihoods P(y_t | y_{1:t-1})
- `s_pred`: `Ns` x `Nt` matrix of one-step predicted state vectors s_{t|t-1}
- `P_pred`: `Ns` x `Ns` x `Nt` array of mean squared errors P_{t|t-1} of
  predicted state vectors
- `s_filt`: `Ns` x `Nt` matrix of filtered state vectors s_{t|t}
- `P_filt`: `Ns` x `Ns` x `Nt` matrix containing mean squared errors P_{t|t} of
  filtered state vectors
- `s_0`: `Ns` x 1 initial state vector. This may have been reassigned to the
  last presample state vector if `Nt0 > 0`
- `P_0`: `Ns` x `Ns` initial state covariance matrix. This may have been
  reassigned to the last presample state covariance if `Nt0 > 0`
- `s_T`: `Ns` x 1 final filtered state `s_{T|T}`
- `P_T`: `Ns` x `Ns` final filtered state covariance matrix `P_{T|T}`

Compute the one-step-ahead states s_{t|t-1} and state covariances P_{t|t-1} and
assign to `k`.
"""
function ct_block_kalman_filter(y::Matrix{S}, T::Matrix{S}, R::Matrix{S}, C::Vector{S},
                          Q::Matrix{S}, Z::Matrix{S}, D::Vector{S}, E::Matrix{S};
                          n_simulate_states::Int = 1, s_0::Vector{S} = Vector{S}(undef, 0),
                          P_0::Matrix{S} = Matrix{S}(undef, 0, 0),
                          outputs::Vector{Symbol} = [:loglh, :pred, :filt],
                          Nt0::Int = 0) where {S<:AbstractFloat}
    # Check s_0, P_0 right dimensions
    if size(s_0, 1) != n_simulate_states * size(T, 1)
        error("Dimension error: s_0 does not stack intermediate simulated states.")
    elseif size(P_0, 1) != n_simulate_states * size(T, 1)
        error("Dimension error: P_0 does correspond to a state vector with stacked intermediate simulated states.")
    end

    # Determine outputs
    return_loglh = :loglh in outputs
    return_pred  = :pred  in outputs
    return_filt  = :filt  in outputs

    # Initialize inputs and outputs
    k = CTBlockKalmanFilter(T, R, C, Q, Z, D, E, s_0, P_0; n_simulate_states)

    # Dimensions
    Ns = size(T, 1) * n_simulate_states # number of states
    Nt = size(y, 2)                     # number of periods of data

    mynan = convert(S, NaN)
    loglh  = return_loglh ? fill(mynan, Nt)         : Vector{S}(0)
    s_pred = return_pred  ? fill(mynan, Ns, Nt)     : Matrix{S}(0, 0)
    P_pred = return_pred  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(0, 0, 0)
    s_filt = return_filt  ? fill(mynan, Ns, Nt)     : Matrix{S}(0, 0)
    P_filt = return_filt  ? fill(mynan, Ns, Ns, Nt) : Array{S, 3}(0, 0, 0)

    # Populate initial states
    s_0 = k.s_t
    P_0 = k.P_t # note that this is now a dict

    # Loop through periods t
    for t = 1:Nt
        # Forecast
        forecast!(k)
        if return_pred
            s_pred[:,    t] = k.s_t
            P_tmp = zeros(Ns, Ns)
            T_size = size(T, 1)
            for i = 1:n_simulate_states
                for j = 1:n_simulate_states
                    if i < j
                        P_tmp[1 + T_size * (i-1):T_size * i, 1 + T_size * (j-1):T_size * j] = k.P_t[(j,i)]'
                    else
                        P_tmp[1 + T_size * (i-1):T_size * i, 1 + T_size * (j-1):T_size * j] = k.P_t[(i,j)]
                    end
                end
            end
            P_pred[:, :, t] = P_tmp
        end

        # Update and compute log-likelihood
        update!(k, y[:, t]; return_loglh = return_loglh)
        if return_filt
            s_filt[:, t] = k.s_t
            P_tmp = zeros(Ns, Ns)
            T_size = size(T, 1)
            for i = 1:n_simulate_states
                for j = 1:n_simulate_states
                    if i < j # get transpose of upper triangle
                        P_tmp[1 + T_size * (i-1):T_size * i, 1 + T_size * (j-1):T_size * j] = k.P_t[(j,i)]'
                    else     # held in P_t as an element
                        P_tmp[1 + T_size * (i-1):T_size * i, 1 + T_size * (j-1):T_size * j] = k.P_t[(i,j)]
                    end
                end
            end
            P_filt[:, :, t] = P_tmp
        end
        if return_loglh
            loglh[t]        = k.loglh_t
        end

        # # Update s_0 and P_0 if Nt0 > 0
        # if t == Nt0
        #     s_0 = k.s_t
        #     P_0 = k.P_t
        # end
    end

    # Populate final states
    s_T = k.s_t
    P_T = k.P_t

    # Remove presample periods
    loglh, s_pred, P_pred, s_filt, P_filt =
        remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt; outputs = outputs)

    return loglh, s_pred, P_pred, s_filt, P_filt, s_0, P_0, s_T, P_T
end

# Computes the one-step-ahead states s_{t|t-1} and state covariances P_{t|t-1}
# and assign to 'k'
function forecast!(k::CTBlockKalmanFilter{Float64})
    T_powers, C, RQRp = k.T_powers, k.C, k.RQRp
    orig_dim = size(T_powers[1], 1)
    s_filt_prev = copy(k.s_t[end - orig_dim:orig_dim]) # don't need all of s_filt
    P_filt = k.P_t

    for i = 1:k.n_simulate_states
        k.s_t[1 + (i - 1) * orig_dim:orig_dim * i] = T_powers[i] * s_filt_prev + C
        for j = i:k.n_simulate_states
            k.P_t[(i,j)] = T_powers[i] * P_filt[(n_simulate_states, n_simulate_states)] * T_powers[j]' + RQRp[(i,j)]
        end
    end

    return nothing
end

function update!(k::CTBlockKalmanFilter{Float64}, y_obs::Vector{Float64}; return_loglh::Bool = true)
    # Keep rows of measurement equation corresponding to non-NaN observables
    nonnan = .!isnan.(y_obs)
    y_obs = y_obs[nonnan]
    Z = k.Z[nonnan, :]
    D = k.D[nonnan]
    E = k.E[nonnan, nonnan]
    Ny = length(y_obs)
    s_pred = k.s_t
    P_pred = k.P_t
    T_size = size(k.T_powers[1], 1)

    y_pred = Z * s_pred  + D
    V_pred = Z * P_pred * Z' + E
    V_pred = (V_pred + V_pred')/2
    V_pred_inv = inv(V_pred)
    dy = y_obs - y_pred
    PZV = P_pred' * Z' * V_pred_inv

    k.s_t = s_pred * PZV * dy
    k.P_t = P_pred - PZV * Z * P_pred

    if return_loglh
        k.loglh_t = -(Ny * log(2*π) + log(det(V_pred)) + dy'*V_pred_inv*dy)/2
    end

end

# """
# ```
# remove_presample!(Nt0, loglh, s_pred, P_pred, s_filt, P_filt)
# ```

# Remove the first `Nt0` periods from all other input arguments and return.
# """
# function remove_presample!(Nt0::Int, loglh::Vector{S},
#                            s_pred::Matrix{S}, P_pred::Array{S, 3},
#                            s_filt::Matrix{S}, P_filt::Array{S, 3};
#                            outputs::Vector{Symbol} = [:loglh, :pred, :filt]) where {S<:AbstractFloat}
#     if Nt0 > 0
#         if :loglh in outputs
#             loglh  = loglh[(Nt0+1):end]
#         end
#         if :pred in outputs
#             s_pred = s_pred[:,    (Nt0+1):end]
#             P_pred = P_pred[:, :, (Nt0+1):end]
#         end
#         if :filt in outputs
#             s_filt = s_filt[:,    (Nt0+1):end]
#             P_filt = P_filt[:, :, (Nt0+1):end]
#         end
#     end
#     return loglh, s_pred, P_pred, s_filt, P_filt
# end


text:
import RemoveLFS

const GIT_HOST = lowercase(strip(get(ENV, "GIT_HOST", "",)))

if GIT_HOST == "github"
    const GITHUB_DST_ORGANIZATION = "UnofficialJuliaMirrorSnapshots"
    const GITHUB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const GITHUB_BOT_PERSONAL_ACCESS_TOKEN =
        ENV["GITHUB_BOT_PERSONAL_ACCESS_TOKEN"]
    const dst_provider =
        RemoveLFS.Hosts.GitHubDstHost.new_github_session(
            ;
            github_organization = GITHUB_DST_ORGANIZATION,
            github_bot_username = GITHUB_BOT_USERNAME,
            github_bot_personal_access_token =
                GITHUB_BOT_PERSONAL_ACCESS_TOKEN,
            )
elseif GIT_HOST == "gitlab"
    const GITLAB_DST_GROUP = "UnofficialJuliaMirrorSnapshots"
    const GITLAB_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const GITLAB_BOT_PERSONAL_ACCESS_TOKEN =
        ENV["GITLAB_BOT_PERSONAL_ACCESS_TOKEN"]
    const dst_provider =
        RemoveLFS.Hosts.GitLabDstHost.new_gitlab_session(
            ;
            gitlab_group = GITLAB_DST_GROUP,
            gitlab_bot_username = GITLAB_BOT_USERNAME,
            gitlab_bot_personal_access_token =
                GITLAB_BOT_PERSONAL_ACCESS_TOKEN,
            )
elseif GIT_HOST == "bitbucket"
    const BITBUCKET_DST_TEAM = "UnofficialJuliaMirrorSnapshots"
    const BITBUCKET_BOT_USERNAME = "UnofficialJuliaMirrorBot"
    const BITBUCKET_BOT_APP_PASSWORD =
        ENV["BITBUCKET_BOT_APP_PASSWORD"]
    const dst_provider =
        RemoveLFS.Hosts.BitbucketDstHost.new_bitbucket_session(
            ;
            bitbucket_team = BITBUCKET_DST_TEAM,
            bitbucket_bot_username = BITBUCKET_BOT_USERNAME,
            bitbucket_bot_app_password = BITBUCKET_BOT_APP_PASSWORD,
            )
else
    error(
        string(
            "You must set the \"GIT_HOST\" environment variable, ",
            "and it must be set to one of the following values: ",
            "GITHUB GITLAB BITBUCKET",
            )
        )
end


text:
using CPLEX

include(joinpath(Pkg.dir("MathProgBase"),"test","linprog.jl"))
linprogtest(CplexSolver(CPX_PARAM_PREIND=0, CPX_PARAM_LPMETHOD=2, CPX_PARAM_SCRIND=0))

include(joinpath(Pkg.dir("MathProgBase"),"test","mixintprog.jl"))
mixintprogtest(CplexSolver(CPX_PARAM_SCRIND=0))

include(joinpath(Pkg.dir("MathProgBase"),"test","quadprog.jl"))
quadprogtest(CplexSolver(CPX_PARAM_SCRIND=0))
socptest(CplexSolver(CPX_PARAM_SCRIND=0))

solver = CplexSolver(CPX_PARAM_SCRIND=0)
MathProgBase.setparameters!(solver, Silent=true, TimeLimit=100.0)
include(joinpath(Pkg.dir("MathProgBase"),"test","linproginterface.jl"))
linprogsolvertest(solver)

include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
coniclineartest(CplexSolver(CPX_PARAM_SCRIND=0))
conicSOCtest(CplexSolver(CPX_PARAM_SCRIND=0))
conicSOCRotatedtest(CplexSolver(CPX_PARAM_SCRIND=0))
conicSOCINTtest(CplexSolver(CPX_PARAM_SCRIND=0))


text:
using LinearAlgebra

"""
    filter_cube(snap_file::String, corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real},
                          parttype::Integer)

Reads positions from `snap_file` and returns the indices of particles contained in a box defined by `corner_lowerleft` and `corner_upperright`.
"""
function filter_cube(snap_file::String, corner_lowerleft::Array{<:Real}, corner_upperright::Array{<:Real}; parttype::Integer=0)

    # read positions from file
    pos = read_snap(snap_file, "POS", parttype)

    # double check if the corners are set correctly
    cube_lower = Array{eltype(corner_lowerleft[1]),1}(undef, 3)
    cube_upper = Array{eltype(corner_lowerleft[1]),1}(undef, 3)
    @inbounds for dim = 1:3
        cube_lower[dim] = ((corner_lowerleft[dim] < corner_upperright[dim]) ? corner_lowerleft[dim] : corner_upperright[dim])
        cube_upper[dim] = ((corner_lowerleft[dim] > corner_upperright[dim]) ? corner_lowerleft[dim] : corner_upperright[dim])
    end

    return findall( @views ( ( cube_lower[1] .<= pos[1,:] .<= cube_upper[1] ) .&
                             ( cube_lower[2] .<= pos[2,:] .<= cube_upper[2] ) .&
                             ( cube_lower[3] .<= pos[3,:] .<= cube_upper[3] )) )
end

"""
    check_in_cylinder(x, pt1, pt2, r)

Checks if a 3D point `x` is in a cylinder defined by its endpoints `pt1` and `pt2` and radius `r`.
"""
function check_in_cylinder(x, pt1, pt2, r)
    # https://stackoverflow.com/questions/47932955/how-to-check-if-a-3d-point-is-inside-a-cylinder
    ( ( x .- pt1 ) ⋅ ( pt2 .- pt1 ) ≥ 0.0                             ) &
    ( ( x .- pt2 ) ⋅ ( pt2 .- pt1 ) ≤ 0.0                             ) &
    ( ( norm( (x .- pt1) × (pt2 .- pt1 ) ) / norm( pt2 .- pt1 ) ) ≤ r )
end

"""
    filter_cylinder(filename::String, pt1::Array{<:Real}, pt2::Array{<:Real}, r::Real)

Reads the positions contained in a file and returns the indices of particles contained in a cylinder defined by the endpoints `pt1` and `pt2` and radius `r`.
"""
function filter_cylinder(filename::String, pt1::Array{<:Real}, pt2::Array{<:Real}, r::Real; parttype::Integer=0)

    in_cube = filter_cube(filename, pt1, pt2, parttype=parttype )

    if size(in_cube,1) > 0
        check_in_cylinder_helper(x) = check_in_cylinder(x, pt1, pt2, r)
        # read positions from file
        pos = read_snap(filename, "POS", parttype)

        pos_cylinder = @views pos[:, in_cube]
        is_in_cylinder = mapslices(check_in_cylinder_helper, pos_cylinder, dims=1)
        return in_cube[findall( is_in_cylinder .== true )]
    else
        return in_cube
    end
end

"""
    filter_sphere(filename::String, center::Array{<:Real}, r::Real; parttype::Integer=0)

Reads the positions contained in a file and returns the indices of particles contained in a cylinder defined by the endpoints `pt1` and `pt2` and radius `r`.
"""
function filter_sphere(filename::String, center::Array{<:Real}, r::Real; parttype::Integer=0)

    in_cube = filter_cube(filename, center .- r, center .+ r, parttype=parttype)

    if size(in_cube,1) > 0

        distance = @views @. √( (pos[1, in_cube] - center[1])^2 +
                                (pos[2, in_cube] - center[2])^2 +
                                (pos[3, in_cube] - center[3])^2 )

        return in_cube[findall( distance .<= r )]
    else
        return in_cube
    end
end


text:
function handle_color!(uniform_dict, instance_dict)
    color, udict = if haskey(uniform_dict, :color)
        to_value(uniform_dict[:color]), uniform_dict
    elseif haskey(instance_dict, :color)
        to_value(instance_dict[:color]), instance_dict
    else
        nothing, uniform_dict
    end
    if color isa Colorant ||
              color isa AbstractVector{<:Colorant} ||
              color === nothing
        delete!(uniform_dict, :colormap)
    elseif color isa AbstractArray{<:Real}
        udict[:color] = lift(x -> convert(Vector{Float32}, x), udict[:color])
        # udict[:color] = lift(x -> convert(Vector{Float32}, x), udict[:color])
        uniform_dict[:color_getter] = """
            vec4 get_color(){
                vec2 norm = get_colorrange();
                float normed = (color - norm.x) / (norm.y - norm.x);
                return texture(colormap, vec2(normed, 0));
            }
        """
    end
end

const IGNORE_KEYS = Set([:shading, :overdraw, :rotation, :distancefield, :markerspace,
                         :fxaa, :visible, :transformation, :alpha, :linewidth,
                         :transparency, :marker, :lightposition])

function create_shader(scene::Scene, plot::MeshScatter)
    # Potentially per instance attributes
    per_instance_keys = (:rotations, :markersize, :color, :intensity)
    per_instance = filter(plot.attributes.attributes) do (k, v)
        return k in per_instance_keys && !(isscalar(v[]))
    end
    per_instance[:offset] = plot[1]

    for (k, v) in per_instance
        per_instance[k] = Buffer(lift_convert(k, v, plot))
    end

    uniforms = filter(plot.attributes.attributes) do (k, v)
        return (!haskey(per_instance, k)) && isscalar(v[])
    end

    uniform_dict = Dict{Symbol,Any}()
    for (k, v) in uniforms
        k in IGNORE_KEYS && continue
        uniform_dict[k] = lift_convert(k, v, plot)
    end

    handle_color!(uniform_dict, per_instance)

    instance = normal_mesh(plot.marker[])

    if !hasproperty(instance, :uv)
        uniform_dict[:uv] = Vec2f0(0)
    end

    return InstancedProgram(WebGL(), lasset("particles.vert"), lasset("particles.frag"),
                            instance, VertexArray(; per_instance...); uniform_dict...)
end

@enum Shape CIRCLE RECTANGLE ROUNDED_RECTANGLE DISTANCEFIELD TRIANGLE

primitive_shape(::Union{String,Char,Vector{Char}}) = Cint(DISTANCEFIELD)
primitive_shape(x::X) where {X} = Cint(primitive_shape(X))
primitive_shape(::Type{<:Circle}) = Cint(CIRCLE)
primitive_shape(::Type{<:Rect2D}) = Cint(RECTANGLE)
primitive_shape(::Type{T}) where {T} = error("Type $(T) not supported")
primitive_shape(x::Shape) = Cint(x)

using AbstractPlotting: to_spritemarker

function scatter_shader(scene::Scene, attributes)
    # Potentially per instance attributes
    per_instance_keys = (:offset, :rotations, :markersize, :color, :intensity,
                         :uv_offset_width, :marker_offset)
    uniform_dict = Dict{Symbol,Any}()

    if haskey(attributes, :marker) && attributes[:marker][] isa Union{Vector{Char},String}
        x = pop!(attributes, :marker)
        attributes[:uv_offset_width] = lift(x -> AbstractPlotting.glyph_uv_width!.(collect(x)),
                                            x)
        uniform_dict[:shape_type] = Cint(3)
    end

    per_instance = filter(attributes) do (k, v)
        return k in per_instance_keys && !(isscalar(v[]))
    end

    for (k, v) in per_instance
        per_instance[k] = Buffer(lift_convert(k, v, nothing))
    end

    uniforms = filter(attributes) do (k, v)
        return !haskey(per_instance, k)
    end

    for (k, v) in uniforms
        k in IGNORE_KEYS && continue
        uniform_dict[k] = lift_convert(k, v, nothing)
    end

    get!(uniform_dict, :shape_type) do
        return lift(x -> primitive_shape(to_spritemarker(x)), attributes[:marker])
    end
    if uniform_dict[:shape_type][] == 3
        atlas = AbstractPlotting.get_texture_atlas()
        uniform_dict[:distancefield] = Sampler(atlas.data, minfilter=:linear,
                                               magfilter=:linear, anisotropic=16f0)
        uniform_dict[:atlas_texture_size] = Float32(size(atlas.data, 1)) # Texture must be quadratic
    else
        uniform_dict[:atlas_texture_size] = 0f0
        uniform_dict[:distancefield] = Observable(false)
    end

    if !haskey(per_instance, :uv_offset_width)
        get!(uniform_dict, :uv_offset_width) do
            return if haskey(attributes, :marker) &&
                      to_spritemarker(attributes[:marker][]) isa Char
                lift(x -> AbstractPlotting.glyph_uv_width!(to_spritemarker(x)),
                     attributes[:marker])
            else
                Vec4f0(0)
            end
        end
    end

    space = get(uniforms, :markerspace, Observable(SceneSpace))
    uniform_dict[:use_pixel_marker] = map(space) do space
        return space == Pixel
    end
    handle_color!(uniform_dict, per_instance)

    instance = uv_mesh(Rect2D(-0.5f0, -0.5f0, 1f0, 1f0))

    uniform_dict[:resolution] = scene.camera.resolution
    return InstancedProgram(WebGL(), lasset("simple.vert"), lasset("sprites.frag"),
                            instance, VertexArray(; per_instance...); uniform_dict...)
end

function create_shader(scene::Scene, plot::Scatter)
    # Potentially per instance attributes
    per_instance_keys = (:offset, :rotations, :markersize, :color, :intensity,
                         :marker_offset)
    per_instance = filter(plot.attributes.attributes) do (k, v)
        return k in per_instance_keys && !(isscalar(v[]))
    end
    attributes = copy(plot.attributes.attributes)
    attributes[:offset] = plot[1]
    attributes[:billboard] = map(rot -> isa(rot, Billboard), plot.rotations)
    attributes[:pixelspace] = getfield(scene.camera, :pixel_space)
    attributes[:model] = plot.model
    attributes[:markerspace] = plot.markerspace
    delete!(attributes, :uv_offset_width)
    return scatter_shader(scene, attributes)
end

function to_gl_text(string, positions_per_char::AbstractVector{T}, textsize,
                    font, align, rot, model, j, l) where T <: VecTypes
    atlas = get_texture_atlas()
    N = length(T)
    positions, uv_offset_width, scale = Point{3, Float32}[], Vec4f0[], Vec2f0[]
    char_str_idx = iterate(string)
    offsets = Vec2f0[]
    broadcast_foreach(1:length(string), positions_per_char, textsize, font, align) do idx, pos, tsize, font, align
        char, str_idx = char_str_idx
        mpos = model * Vec4f0(to_ndim(Vec3f0, pos, 0f0)..., 1f0)
        push!(positions, to_ndim(Point{3, Float32}, mpos, 0))
        push!(uv_offset_width, glyph_uv_width!(atlas, char, font))
        glyph_bb, ext = FreeTypeAbstraction.metrics_bb(char, font, tsize)
        if isa(tsize, Vec2f0) # this needs better unit support
            push!(scale, tsize) # Vec2f0, we assume it's already in absolute size
        else
            push!(scale, widths(glyph_bb))
        end
        push!(offsets, minimum(glyph_bb))
        char_str_idx = iterate(string, str_idx)
    end
    return positions, offsets, uv_offset_width, scale
end

function to_gl_text(string, startpos::VecTypes{N, T}, textsize, font, aoffsetvec, rot, model, j, l) where {N, T}
    atlas = get_texture_atlas()
    positions = layout_text(string, startpos, textsize, font, aoffsetvec, rot, model, j, l)
    uv = Vec4f0[]
    scales = Vec2f0[]
    offsets = Vec2f0[]
    for (c, font, pixelsize) in zip(string, attribute_per_char(string, font), attribute_per_char(string, textsize))
        push!(uv, glyph_uv_width!(atlas, c, font))
        glyph_bb, extent = FreeTypeAbstraction.metrics_bb(c, font, pixelsize)
        push!(scales, widths(glyph_bb))
        push!(offsets, minimum(glyph_bb))
    end
    return positions, offsets, uv, scales
end

function create_shader(scene::Scene, plot::AbstractPlotting.Text)

    liftkeys = (:position, :textsize, :font, :align, :rotation, :model, :justification, :lineheight)
    args = getindex.(Ref(plot), liftkeys)
    gl_text = lift(plot[1], args...) do str, pos, tsize, font, align, rotation, model, j, l
        # For annotations, only str (x[1]) will get updated, but all others are updated too!
        args = @get_attribute plot (position, textsize, font, align, rotation)
        to_gl_text(str, args..., model, j, l)
    end

    # unpack values from the one signal:
    positions, offset, uv_offset_width, scale = map((1, 2, 3, 4)) do i
        return lift(getindex, gl_text, i)
    end
    # Sigh, we also allow inplace mutation without triggering
    # plot.rotation to update, so we need to update also on the string array (plot[1])
    rotation = lift(plot[1], plot.rotation) do str, rotation
        return to_rotation(rotation)
    end
    color = lift(plot[1], plot.color) do str, color
        return to_color(color)
    end
    uniforms = Dict(:model => Observable(Mat4f0(I)), :shape_type => Observable(Cint(3)),
                    :color => color, :rotations => rotation, :markersize => scale,
                    :markerspace => Observable(SceneSpace),
                    :marker_offset => offset, :offset => positions,
                    :uv_offset_width => uv_offset_width,
                    :transform_marker => Observable(false), :billboard => Observable(false),
                    :pixelspace => getfield(scene.camera, :pixel_space))
    return scatter_shader(scene, uniforms)
end


text:
# String Vector with the tables used in 'query'
function qtables(query::String)::Vector{String}

    # Normalize Query
    query = replace(query,"\n"=>" ")

    nq = split(query," ")
    nq = nq[nq .!= ""]

    # Extract tables in query
    tables = String[]
    nt = false
    for w in nq
        nt = nt ? push!(tables,w)==0 : lowercase(w) in ["from","join"]
    end
    tables

end


text:
using Test

module X1
    using Reexport
    @reexport module Y1
        const Z1 = 1
        export Z1
    end
end

@test union!(Set(), names(X1)) == union!(Set(), [:X1, :Y1, :Z1])
@test X1.Z1 == 1

module Y2
    const Z2 = 2
    export Z2
end
module X2
    using Reexport
    # Locally defined modules require a prefix for lookup in 0.7. It isn't necessary
    # on older versions but doesn't hurt.
    @reexport using Main.Y2
end
@test union!(Set(), names(X2)) == union!(Set(), [:X2, :Y2, :Z2])
@test X2.Z2 == 2

module X3
    using Reexport
    module Y3
        const Z3 = 3
        export Z3
    end
    module Y4
        const Z4 = 4
        export Z4
    end
    @reexport using .Y3, .Y4
end
@test union!(Set(), names(X3)) == union!(Set(), [:X3, :Y3, :Y4, :Z3, :Z4])
@test X3.Z3 == 3
@test X3.Z4 == 4

module X4
    using Reexport
    @reexport using Main.Y2
end
@test union!(Set(), names(X4)) == union!(Set(), [:X4, :Y2, :Z2])
@test X4.Z2 == 2

module X5
    using Reexport
    module Y3
        const Z3 = 3
        export Z3
    end
    module Y4
        const Z4 = 4
        export Z4
    end
    @reexport using .Y3, .Y4
end
@test union!(Set(), names(X3)) == union!(Set(), [:X3, :Y3, :Y4, :Z3, :Z4])
@test X3.Z3 == 3
@test X3.Z4 == 4

module X6
    using Reexport
    module Y5
        const Z5 = 5
        export Z5
        const Z6 = 6
        export Z6
    end
    @reexport using .Y5: Z5, Z6
end
@test union!(Set(), names(X6)) == union!(Set(), [:X6, :Z5, :Z6])
@test X6.Z5 == 5
@test X6.Z6 == 6

using .X6
@test Z5 == 5
@test Z6 == 6

module X7
    using Reexport
    module Y7
        struct T7 end
        export T7
        Base.@deprecate_binding S7 T7
    end
    @reexport using .Y7
end
using .X7
@test Base.isexported(X7, :S7)

#== Imports ==#
module X8
    using Reexport

    module InnerX8
        const a = 1
        export a
    end
    @reexport import .InnerX8: a
end

module X9
    using Reexport

    module InnerX9
        const a = 1
        export a
    end
    @reexport import .InnerX9.a
end

module X10
    using Reexport

    module InnerX10_1
        const a = 1
        export a
    end

    module InnerX10_2
        const b = 1
        export b
    end

    @reexport import .InnerX10_1.a, .InnerX10_2.b
end

module X11
    using Reexport

    module InnerX11
        const b = 1
        export b
    end
    @reexport import .InnerX11
end

@testset "import" begin
    @testset "Reexported colon-qualified single import" begin
        @test Set(names(X8)) == Set([:X8, :a])
    end

    @testset "Reexported dot-qualified single import" begin
        @test Set(names(X9)) == Set([:X9, :a])
    end

    @testset "Reexported qualified multiple import" begin
        @test Set(names(X10)) == Set([:X10, :a, :b])
    end

    @testset "Reexported module import" begin
        @test Set(names(X11)) == Set([:X11, :InnerX11])
    end
end

#== block ==#
module X12
    using Reexport
    @reexport begin
        using Main.X9
        using Main.X10
    end
end

module X13
    using Reexport
    @reexport begin
        import Main.X9
        import Main.X10
    end
end

module X14
    using Reexport
    module InnerX14
        const a = 1
        export a
    end
    @reexport begin
        import Main.X9
        using Main.X10
        using .InnerX14: a
    end
end

@testset "block" begin
    @testset "block of using" begin
        @test Set(names(X12)) == union(Set(names(X9)), Set(names(X10)), Set([:X12]))
    end
    @testset "block of import" begin
        @test Set(names(X13)) == Set([:X13, :X9, :X10])
    end
    @testset "mixed using and import" begin
        @test Set(names(X14)) == union(Set([:X14, :X9, :a]), Set(names(X10)))
    end
end

#== macroexpand ==#
module X15
    using Reexport

    macro identity_macro(ex::Expr)
        ex
    end

    module InnerX15
        const a = 1
        export a
    end

    @reexport @identity_macro using .InnerX15: a
end
@testset "macroexpand" begin
    @test Set(names(X15)) == Set([:X15, :a])
end

module X16
    using Reexport, Test
    baremodule A
        using Reexport
        @reexport using Test
    end
    baremodule B
        using Reexport
        @reexport using Test: @testset
    end
    baremodule C
        using Reexport
        @reexport import Test.@test
    end
    baremodule D
        using Reexport
        # Ensure that there are no module name conflicts with Core
        baremodule Core
        end
    end
    @testset "baremodule" begin
        @test Base.isexported(A, Symbol("@test"))
        @test Reexport.exported_names(B) == [Symbol("@testset"), :B]
        @test Reexport.exported_names(C) == [Symbol("@test"), :C]
        @test Reexport.exported_names(D) == [:D]
    end
end


text:
using JLD
using PyPlot
@load "/Users/milan/cams/error/linvslog_unstructured.jld"

nmethods,nvars,nstats = size(linerror)

# sort by groups
aero = Array(1:15)
ozone = [35,44,45]
methane = [25,26,41]
dynamics = [33,34,54,55,53]
clouds = [22,27,28,31,32,51]
hydro = vcat(Array(36:40),46)
nitro = [42,43,52]
oopp = Array(47:50)
ces = vcat(Array(17:21),[23,24,16])
co12 = [29,30]
chydro = vcat(ces,hydro)

grouped = vcat(aero,co12,clouds,methane,chydro,dynamics,nitro,ozone,oopp)
groups = [aero,co12,clouds,methane,chydro,dynamics,nitro,ozone,oopp]

allvars = allvars[reverse(grouped)]
meanerror = meanerror[:,reverse(grouped)]
linerror = linerror[:,reverse(grouped),:]
logerror = logerror[:,reverse(grouped),:]

allvars[allvars .== "c"] .= "ch4_c"

## fit histogram
ioff()
fig,(ax1,ax2,ax3) = subplots(1,3,figsize=(10,8),sharey=true)

lwmm = 1
lw10 = 2
lw25 = 3

colours = ["C0","k","C1","C3"]

for i in 1:nvars

    # MEAN ERROR
    for j in 1:nmethods
        me = meanerror[j,i]
        y = i-j/5
        ax1.set_xscale("log")
        ax1.scatter(me,y,30,colours[j],edgecolor="k",zorder=5)
        ax1.axhline(i-1,color="grey",lw=0.5)
    end

    # ABSOLUTE ERROR
    for j in 1:nmethods
        mi,p10,p25,p50,p75,p90,ma = linerror[j,i,:]
        y = i-j/5
        ax2.semilogx([p50,ma],[y,y],colours[j],lw=lwmm)
        ax2.semilogx([p50,p90],[y,y],colours[j],lw=lw10)
        ax2.scatter(p50,y,30,colours[j],edgecolor="k",zorder=5)
        ax2.axhline(i-1,color="grey",lw=0.5)
    end

    # DECIMAL ERROR
    for j in 1:nmethods
        mi,p10,p25,p50,p75,p90,ma = logerror[j,i,:]
        y = i-j/5
        ax3.semilogx([p50,ma],[y,y],colours[j],lw=lwmm)
        ax3.semilogx([p50,p90],[y,y],colours[j],lw=lw10)
        ax3.scatter(p50,y,30,colours[j],edgecolor="k",zorder=5)
        ax3.axhline(i-1,color="grey",lw=0.5)
    end

end

for ax in (ax1,ax2,ax3)
    for x in [4,7,10,15,32,38,40]
        ax.axhline(x,color="grey",lw=2)
    end
end

ax1.set_yticks(Array(1:nvars).-0.5)
ax1.yaxis.set_tick_params(length=0)
ax2.yaxis.set_tick_params(length=0)
ax3.yaxis.set_tick_params(length=0)
ax1.set_yticklabels(allvars)
ax1.set_ylim(-6.5,nvars+0.25)

ax1.set_xlim(1e-12,9e-1)
ax2.set_xlim(2e-22,1e1)
ax3.set_xlim(1.1e-8,1e2)

alfa=0.25
for i in 1:nvars
    if meanerror[1,i] >= meanerror[3,i] && meanerror[2,i] < meanerror[4,i]
        ax1.fill_between(ax1.get_xlim(),[i-1,i-1],[i,i],color="C9",alpha=alfa)
    end

    for (ax,e) in zip([ax2,ax3],[linerror,logerror])
        if e[1,i,4] >= e[3,i,4] && e[2,i,4] < e[4,i,4]
            ax.fill_between(ax.get_xlim(),[i-1,i-1],[i,i],color="C9",alpha=alfa)
        end
    end
end

for i in 1:nvars
    if meanerror[2,i] >= meanerror[4,i]
        ax1.fill_between(ax1.get_xlim(),[i-1,i-1],[i,i],color="C2",alpha=alfa)
    end

    for (ax,e) in zip([ax2,ax3],[linerror,logerror])
        if e[2,i,4] >= e[4,i,4]
            ax.fill_between(ax.get_xlim(),[i-1,i-1],[i,i],color="C2",alpha=alfa)
        end
    end
end

mi,p10,p25,p50,p75,p90,ma = linerror[4,1,:]
ax2.text(p50,-0.5,"median",rotation=90,ha="center",va="top",fontsize=9)
ax2.text(p90,-0.5,"90%",rotation=90,ha="center",va="top",fontsize=9)
ax2.text(ma,-0.5,"max",rotation=90,ha="center",va="top",fontsize=9)

ax1.scatter(0,0,30,"C0",edgecolor="k",label="LinQuant16")
ax1.scatter(0,0,30,"k",edgecolor="k",label="LinQuant24")
ax1.scatter(0,0,30,"C1",edgecolor="k",label="LogQuant16")
ax1.scatter(0,0,30,"C3",edgecolor="k",label="LogQuant24")
ax3.fill_between(ax3.get_xlim(),[-10,-10],[-9,-9],color="C9",alpha=alfa,label="at 16 bit")
ax3.fill_between(ax3.get_xlim(),[-10,-10],[-9,-9],color="C2",alpha=alfa,label="at both 16 and 24 bit")
ax1.legend(loc=4,ncol=2,fontsize=8,scatterpoints=3)
ax3.legend(loc=4,fontsize=9,title="Log better than LinQuantization")

ax1.set_xlabel("norm. mean error")
ax2.set_xlabel("norm. absolute error")
ax3.set_xlabel("decimal error")

ax1.set_title("Mean error",loc="left")
ax2.set_title("Absolute error",loc="left")
ax3.set_title("Decimal error",loc="left")

ax1.set_title("a",loc="right",fontweight="bold")
ax2.set_title("b",loc="right",fontweight="bold")
ax3.set_title("c",loc="right",fontweight="bold")

tight_layout()
savefig("/Users/milan/git/Elefridge.jl/plots/linvslog_sorted.png",dpi=200)
close(fig)


text:
function UncorrelatedRandomVectorGenerator(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::NormalizedRandomGenerator)
    return UncorrelatedRandomVectorGenerator((Vector{jdouble}, Vector{jdouble}, NormalizedRandomGenerator), arg0, arg1, arg2)
end

function UncorrelatedRandomVectorGenerator(arg0::jint, arg1::NormalizedRandomGenerator)
    return UncorrelatedRandomVectorGenerator((jint, NormalizedRandomGenerator), arg0, arg1)
end

function next_vector(obj::UncorrelatedRandomVectorGenerator)
    return jcall(obj, "nextVector", Vector{jdouble}, ())
end


text:
# Autogenerated wrapper script for HDF5_jll for aarch64-apple-darwin
export libhdf5, libhdf5_hl

using Zlib_jll
using OpenSSL_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("HDF5")
JLLWrappers.@declare_library_product(libhdf5, "@rpath/libhdf5.200.dylib")
JLLWrappers.@declare_library_product(libhdf5_hl, "@rpath/libhdf5_hl.200.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, OpenSSL_jll, LibCURL_jll)
    JLLWrappers.@init_library_product(
        libhdf5,
        "lib/libhdf5.200.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libhdf5_hl,
        "lib/libhdf5_hl.200.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ArtGallery

using Term
using Dates

include(joinpath(@__DIR__, "types.jl"))
include(joinpath(@__DIR__, "gallery.jl"))

export album

"""
    album(id::String)
This function will show the `id` picture in `REPL`.
id can be one of the following:
- "buddha"
- "alpaca"
- "dragon"
- "victory"
- "monster"
- "computer"
- "kitty"
"""
function album(id::String="computer")
    style = gallery()
    if id == "buddha"
        return highlight(style.s1)
    elseif id == "alpaca"
        return highlight(style.s2)
    elseif id == "dragon"
        return highlight(style.s3)
    elseif id == "victory"
        return highlight(style.s4)
    elseif id == "monster"
        return highlight(style.s5)
    elseif id == "computer"
        return highlight(style.s6)
    elseif id == "kitty"
        return highlight(style.s7)
    else
        return highlight(style.s6)
    end
end # function goodluck()

end # module