text:
export *, +, absolute_basis, absolute_basis_matrix, ambient_space, bad_primes,
       basis_matrix, can_scale_totally_positive, coefficient_ideal, degree, diagonal,
       discriminant, dual, fixed_field, generators, gram_matrix_of_rational_span,
       hasse_invariant, hermitian_lattice, intersect, involution, isdefinite,
       isintegral, islocally_isometric, ismodular, isnegative_definite,
       ispositive_definite, isrationally_equivalent, jordan_decomposition,
       local_basis_matrix, norm, pseudo_matrix, quadratic_lattice, rank,
       rational_span, rescale, scale, volume, witt_invariant, lattice,
       Zlattice, automorphism_group_generators, automorphism_group_order,
       isisometric, islocal_norm, normic_defect, issublattice, issublattice_with_relations

export HermLat, QuadLat

################################################################################
#
#  Verbose and assert scopes
#
################################################################################

add_verbose_scope(:Lattice)

add_assert_scope(:Lattice)

################################################################################
#
#  Ambient space
#
################################################################################

@doc Markdown.doc"""
    has_ambient_space(L::AbsLat) -> Bool

Returns whether the ambient space of $L$ is set.
"""
function has_ambient_space(L::AbsLat)
  return isdefined(L, :space)
end

@doc Markdown.doc"""
    ambient_space(L::AbsLat) -> AbsSpace

Returns the ambient space of $L$. If the ambient space is not known, an
error is raised.
"""
function ambient_space(L::AbsLat)
  if has_ambient_space(L)
    return L.space
  else
    error("Ambient space not defined")
  end
end

################################################################################
#
#  Rational span
#
################################################################################

@doc Markdown.doc"""
    rational_span(L::AbsLat) -> AbsSpace

Returns the rational span of $L$.
"""
rational_span(::AbsLat)

################################################################################
#
#  Diagonal
#
################################################################################

@doc Markdown.doc"""
    diagonal(L::AbsLat) -> Vector

Returns the diagonal of the rational span of $L$.
"""
function diagonal_of_rational_span(L::AbsLat)
  D, _ = _gram_schmidt(gram_matrix_of_rational_span(L), involution(L))
  return diagonal(D)
end

################################################################################
#
#  Module properties
#
################################################################################

@doc Markdown.doc"""
    pseudo_matrix(L::AbsLat) -> PMat

Returns the basis pseudo-matrix of $L$.
"""
pseudo_matrix(L::AbsLat) = L.pmat

@doc Markdown.doc"""
    pseudo_basis(L::AbsLat) -> Vector{Tuple{Vector, Ideal}}

Returns the pseudo-basis of $L$.
"""
function pseudo_basis(L::AbsLat)
  M = matrix(pseudo_matrix(L))
  LpM = pseudo_matrix(L)
  O = base_ring(LpM)
  z = Vector{Tuple{Vector{elem_type(nf(O))}, fractional_ideal_type(O)}}(undef, nrows(M))
  for i in 1:nrows(M)
    z[i] = (elem_type(base_ring(M))[ M[i, j] for j in 1:ncols(M) ],
            coefficient_ideals(LpM)[i])
  end
  return z
end

@doc Markdown.doc"""
    coefficient_ideals(L::AbsLat) -> Vector{NfOrdIdl}

Returns the coefficient ideals of the pseudo-basis of $L$.
"""
coefficient_ideals(L::AbsLat) = coefficient_ideals(pseudo_matrix(L))

@doc Markdown.doc"""
    basis_matrix_of_rational_span(L::AbsLat) -> MatElem

Returns the basis matrix of the rational span of $L$.
"""
basis_matrix_of_rational_span(L::AbsLat) = matrix(pseudo_matrix(L))

base_field(L::AbsLat) = L.base_algebra

@doc Markdown.doc"""
    base_ring(L::AbsLat) -> Ring

Return the ring over which the lattice is defined.
"""
base_ring(L::AbsLat) = base_ring(L.pmat)

@doc Markdown.doc"""
    fixed_field(L::AbsLat) -> NumField

Returns the fixed field of the involution of $L$.
"""
fixed_field(L::AbsLat) = fixed_field(rational_span(L))

@doc Markdown.doc"""
    involution(L::AbsLat) -> Map

Returns the involution of the rational span of $L$.
"""
involution(::AbsLat)

@doc Markdown.doc"""
    rank(L::AbsLat) -> Int

Returns the rank of $L$, that is the dimension of the rational span of $L$.
"""
rank(L::AbsLat) = rank(rational_span(L))

@doc Markdown.doc"""
    degree(L::AbsLat) -> Int

Returns the dimension of the ambient space of $L$.
"""
function degree(L::AbsLat)
  return dim(ambient_space(L))
end

@doc Markdown.doc"""
    issublattice(L::AbsLat, M::AbsLat) -> Bool

Returns whether $M$ is a sublattice of $L$.
"""
issublattice(L::AbsLat, M::AbsLat)

################################################################################
#
#  Gram matrix
#
################################################################################

@doc Markdown.doc"""
    gram_matrix_of_rational_span(L::AbsLat) -> MatElem

Returns the gram matrix of the rational span of $L$.
"""
function gram_matrix_of_rational_span(L::AbsLat)
  if isdefined(L, :gram)
    return L.gram
  else
    return gram_matrix(ambient_space(L), L.pmat.matrix)
  end
end

################################################################################
#
#  Generators
#
################################################################################

# Check if one really needs minimal
# Steinitz form is not pretty

@doc Markdown.doc"""
    generators(L::AbsLat; minimal = false) -> Vector{Vector}

Returns a set of generators of $L$ over the base ring of $L$.

If `minimal == true`, the number of generators is minimal. Note that computing
minimal generators is expensive.
"""
function generators(L::AbsLat; minimal::Bool = false)
  K = nf(base_ring(L))
  T = elem_type(K)
  if !minimal
    if isdefined(L, :generators)
      return L.generators::Vector{Vector{T}}
    end
    v = Vector{T}[]
    St = pseudo_matrix(L)
    d = ncols(St)
    for i in 1:nrows(St)
      if base_ring(L) isa NfOrd
        I = numerator(St.coeffs[i])
        den = denominator(St.coeffs[i])
        _assure_weakly_normal_presentation(I)
        push!(v, T[K(I.gen_one)//den * matrix(St)[i, j] for j in 1:d])
        push!(v, T[K(I.gen_two)//den * matrix(St)[i, j] for j in 1:d])
      else
        I = numerator(coefficient_ideals(St)[i])
        den = denominator(coefficient_ideals(St)[i])
        for g in absolute_basis(I)
          push!(v, T[K(g)//den * matrix(St)[i, j] for j in 1:d])
        end
      end
    end
    L.generators = v
    return v
  else # minimal
    if isdefined(L, :minimal_generators)
      return L.minimal_generators::Vector{Vector{T}}
    end
    St = _steinitz_form(pseudo_matrix(L), Val{false})
    d = nrows(St)
    n = degree(L)
    v = Vector{T}[]
    for i in 1:(d - 1)
      #@assert isprincipal(coefficient_ideals(St)[i])[1]
      push!(v, T[matrix(St)[i, j] for j in 1:d])
    end

    I = numerator(coefficient_ideals(St)[d])
    den = denominator(coefficient_ideals(St)[d])
    if minimal && base_ring(L) isa NfOrd
      b, a = isprincipal(I)
      if b
        push!(v, T[K(a)//den * matrix(St)[n, j] for j in 1:d])
      end
      return v
    end

    if base_ring(L) isa NfOrd
      _assure_weakly_normal_presentation(I)
      push!(v, T[K(I.gen_one)//den * matrix(St)[n, j] for j in 1:d])
      push!(v, T[K(I.gen_two)//den * matrix(St)[n, j] for j in 1:d])
    else
      for g in absolute_basis(I)
        push!(v, T[K(g)//den * matrix(St)[n, j] for j in 1:d])
      end
    end
  end

  L.minimal_generators = v

  return v
end

################################################################################
#
#  Gram matrix of generators
#
################################################################################

@doc Markdown.doc"""
    gram_matrix_of_generators(L::AbsLat; minimal::Bool = false) -> MatElem

Returns the Gram matrix of a generating set of $L$. If `minimal` is true,
then a minimal generating set is used.
"""
function gram_matrix_of_generators(L::AbsLat; minimal::Bool = false)
  m = generators(L, minimal = minimal)
  M = matrix(nf(base_ring(L)), m)
  return gram_matrix(ambient_space(L), M)
end

################################################################################
#
#  Discriminant
#
################################################################################

@doc Markdown.doc"""
    discriminant(L::AbsLat) -> NfOrdFracIdl

Returns the discriminant of $L$, that is, the generalized index ideal
$[L^\# : L]$.
"""
function discriminant(L::AbsLat)
  d = det(gram_matrix_of_rational_span(L))
  v = involution(L)
  C = coefficient_ideals(L)
  I = prod(J for J in C)
  return d * I * v(I)
end

################################################################################
#
#  Rational (local) equivalence
#
################################################################################

@doc Markdown.doc"""
    hasse_invariant(L::AbsLat, p::Union{InfPlc, NfOrdIdl}) -> Int

Returns the Hasse invariant of the rational span of $L$ at $p$. The lattice
must be quadratic.
"""
hasse_invariant(L::AbsLat, p)

@doc Markdown.doc"""
    witt_invariant(L::AbsLat, p::Union{InfPlc, NfOrdIdl}) -> Int

Returns the Witt invariant of the rational span of $L$ at $p$. The lattice
must be quadratic.
"""
witt_invariant(L::AbsLat, p)

################################################################################
#
#  Rational equivalent
#
################################################################################

@doc Markdown.doc"""
    isrationally_equivalent(L::AbsLat, M::AbsLat, p::Union{InfPlc, NfOrdIdl})
                                                                         -> Bool

Returns whether the rational spans of $L$ and $M$ are equivalent over the
completion at $\mathfrak p$.
"""
isrationally_equivalent(::AbsLat, ::AbsLat, ::NfAbsOrdIdl)

function isrationally_equivalent(L::AbsLat, M::AbsLat, p::NfAbsOrdIdl)
  return isequivalent(rational_span(L), rational_span(M), p)
end

function isrationally_equivalent(L::AbsLat, M::AbsLat, p::InfPlc)
  return isequivalent(rational_span(L), rational_span(M), p)
end

@doc Markdown.doc"""
    isrationally_equivalent(L::AbsLat, M::AbsLat)
                                            -> Bool

Returns whether the rational spans of $L$ and $M$ are equivalent.
"""
function isrationally_equivalent(L::AbsLat, M::AbsLat)
  return isequivalent(rational_span(L), rational_span(M))
end

################################################################################
#
#  Definiteness
#
################################################################################

@doc Markdown.doc"""
    ispositive_definite(L::AbsLat) -> Bool

Returns whether the rational span of $L$ is positive definite.
"""
ispositive_definite(L::AbsLat) = ispositive_definite(rational_span(L))

@doc Markdown.doc"""
    isnegative_definite(L::AbsLat) -> Bool

Returns whether the rational span of $L$ is negative definite.
"""
isnegative_definite(L::AbsLat) = isnegative_definite(rational_span(L))

@doc Markdown.doc"""
    isdefinite(L::AbsLat) -> Bool

Returns whether the rational span of $L$ is definite.
"""
isdefinite(L::AbsLat) = isdefinite(rational_span(L))

@doc Markdown.doc"""
    can_scale_totally_positive(L::AbsLat) -> Bool, NumFieldElem

Returns whether there is a totally positive rescaled lattice of $L$. If so, the
second return value is an element $a$ such that $L^a$ is totally positive.
"""
function can_scale_totally_positive(L::AbsLat)
  a = _isdefinite(rational_span(L))
  if iszero(a)
    return false, a
  else
    return true, a
  end
end

################################################################################
#
#  Addition
#
################################################################################

# Some of these assertions can be relaxed, in particular in the scaling

@doc Markdown.doc"""
    +(L::AbsLat, M::AbsLat) -> AbsLat

Returns the sum of $L$ and $M$.

The lattices $L$ and $M$ must have the same ambient space.
"""
function Base.:(+)(L::T, M::T) where {T <: AbsLat}
  @assert has_ambient_space(L) && has_ambient_space(M)
  @assert ambient_space(L) === ambient_space(M)
  V = ambient_space(L)
  fr = nrows(pseudo_matrix(L)) == dim(V) || nrows(pseudo_matrix(M)) == dim(V)
  m = _sum_modules(L, pseudo_matrix(L), pseudo_matrix(M), fr)
  return lattice_in_same_ambient_space(L, m)
end

################################################################################
#
#  Intersection
#
################################################################################

@doc Markdown.doc"""
    intersect(L::AbsLat, M::AbsLat) -> AbsLat

Returns the intersection of $L$ and $M$.

The lattices $L$ and $M$ must have the same ambient space.
"""
function intersect(L::T, M::T) where {T <: AbsLat}
  @assert has_ambient_space(L) && has_ambient_space(M)
  @assert ambient_space(L) === ambient_space(M)
  V = ambient_space(L)
  fr = nrows(pseudo_matrix(L)) == dim(V) || nrows(pseudo_matrix(M)) == dim(V)
  if !fr
    error("Only implemented for lattices of full rank")
  end
  m = _intersect_modules(L, pseudo_matrix(L), pseudo_matrix(M), fr)
  return lattice_in_same_ambient_space(L, m)
end

################################################################################
#
#  Scalar multiplication
#
################################################################################

@doc Markdown.doc"""
    *(a::NumFieldElem, M::AbsLat) -> AbsLat

Returns the lattice $aM$ inside the ambient space of $M$.
"""
function Base.:(*)(a::NumFieldElem, L::AbsLat)
  @assert has_ambient_space(L)
  m = _module_scale_element(a, pseudo_matrix(L))
  return lattice_in_same_ambient_space(L, m)
end

function Base.:(*)(L::QuadLat, a)
  return a * L
end

@doc Markdown.doc"""
    *(a::NfOrdIdl, M::AbsLat) -> AbsLat

Returns the lattice $aM$ inside the ambient space of $M$.
"""
function Base.:(*)(a::Union{NfRelOrdIdl, NfAbsOrdIdl}, L::AbsLat)
  @assert has_ambient_space(L)
  m = _module_scale_ideal(a, pseudo_matrix(L))
  return lattice_in_same_ambient_space(L, m)
end

@doc Markdown.doc"""
    *(a::NfOrdFracIdl, M::AbsLat) -> AbsLat

Returns the lattice $aM$ inside the ambient space of $M$.
"""
function Base.:(*)(a::Union{NfRelOrdFracIdl, NfAbsOrdFracIdl}, L::AbsLat)
  @assert has_ambient_space(L)
  m = _module_scale_ideal(a, pseudo_matrix(L))
  return lattice_in_same_ambient_space(L, m)
end

################################################################################
#
#  Absolute basis
#
################################################################################

@doc Markdown.doc"""
    absolute_basis(L::AbsLat) -> Vector

Returns a $\mathbf{Z}$-basis of $L$.
"""
function absolute_basis(L::AbsLat)
  pb = pseudo_basis(L)
  z = Vector{Vector{elem_type(base_field(L))}}()
  for (v, a) in pb
    for w in absolute_basis(a)
      push!(z, w .* v)
    end
  end
  @assert length(z) == absolute_degree(base_field(L)) * rank(L)
  return z
end

################################################################################
#
#  Absolute basis matrix
#
################################################################################

@doc Markdown.doc"""
    absolute_basis_matrix(L::AbsLat) -> MatElem

Returns a $\mathbf{Z}$-basis matrix of $L$.
"""
function absolute_basis_matrix(L::AbsLat)
  pb = pseudo_basis(L)
  E = base_field(L)
  c = ncols(matrix(pseudo_matrix(L)))
  z = zero_matrix(E, rank(L) * absolute_degree(E), c)
  k = 1
  for (v, a) in pb
    for w in absolute_basis(a)
      for j in 1:c
        z[k, j] = w * v[j]
      end
      k += 1
    end
  end
  return z
end

################################################################################
#
#  Norm
#
################################################################################

# cache this
@doc Markdown.doc"""
    norm(L::AbsLat) -> NfOrdFracIdl

Returns the norm of $L$. This is a fractional ideal of the fixed field of $L$.
"""
norm(::AbsLat)

################################################################################
#
#  Scale
#
################################################################################

@doc Markdown.doc"""
    scale(L::AbsLat) -> NfOrdFracIdl

Returns the scale of $L$.
"""
scale(L::AbsLat)

################################################################################
#
#  Rescale
#
################################################################################

@doc Markdown.doc"""
    rescale(L::AbsLat, a::NumFieldElem) -> AbsLat

Returns the rescaled lattice $L^a$. Note that this has a different ambient
space than $L$.
"""
rescale(::AbsLat, ::NumFieldElem)

Base.:(^)(L::AbsLat, a::RingElement) = rescale(L, a)

################################################################################
#
#  Is integral
#
################################################################################

@doc Markdown.doc"""
    isintegral(L::AbsLat) -> Bool

Returns whether the lattice $L$ is integral.
"""
function isintegral(L::AbsLat)
  return isintegral(scale(L))
end

################################################################################
#
#  Dual lattice
#
################################################################################

@doc Markdown.doc"""
    dual(L::AbsLat) -> AbsLat

Returns the dual lattice of $L$.
"""
dual(::AbsLat)

################################################################################
#
#  Volume
#
################################################################################

@doc Markdown.doc"""
    volume(L::AbsLat) -> NfOrdFracIdl

Returns the volume of $L$.
"""
function volume(L::AbsLat)
  return discriminant(L)
end

################################################################################
#
#  Modularity
#
################################################################################

@doc Markdown.doc"""
    ismodular(L::AbsLat) -> Bool, NfOrdFracIdl

Returns whether $L$ is modular. In this case, the second return value is a
fractional ideal $\mathfrak a$ such that $\mathfrak a L^\# = L$, where
$L^\#$ is the dual of $L$.
"""
function ismodular(L::AbsLat)
  a = scale(L)
  if volume(L) == a^rank(L)
    return true, a
  else
    return false, a
  end
end

@doc Markdown.doc"""
    ismodular(L::AbsLat, p::NfOrdIdl) -> Bool, Int

Returns whether $L_{\mathfrak{p}}$ is modular. In this case the second return
value is an integer $v$ such that $L_{\mathfrak{p}}$ is
$\mathfrak{p}^v$-modular.
"""
function ismodular(L::AbsLat{<: NumField}, p)
  a = scale(L)
  if base_ring(L) == order(p)
    v = valuation(a, p)
    if v * rank(L) == valuation(volume(L), p)
      return true, v
    else
      return false, 0
    end
  else
    @assert base_ring(base_ring(L)) == order(p)
    q = prime_decomposition(base_ring(L), p)[1][1]
    v = valuation(a, q)
    if v * rank(L) == valuation(volume(L), q)
      return true, v
    else
      return false, 0
    end
  end
end

################################################################################
#
#  Local basis matrix
#
################################################################################

# The docstring is confusing.
# If p is a prime ideal of base_ring(L), then it actually does
# local_basis_matrix(L, minimum(p),...)
@doc Markdown.doc"""
    local_basis_matrix(L::AbsLat, p::NfOrdIdl; type = :any) -> MatElem

Given a prime ideal $\mathfrak p$ and a lattice $L$, this function returns
a basis matrix of lattice $M$ such that $M_{\mathfrak{p}} = L_{\mathfrak{p}}$.

- If `type == :submodule`, the lattice $L$ will be a sublattice of $M$.
- If `type == :supermodule`, the lattice $L$ will be a superlattice of $M$.
- If `type == :any`, there may not be any containment relation between $M$ and
  $L$.
"""
function local_basis_matrix(L::AbsLat, p; type::Symbol = :any)
  R = base_ring(L)
  S = order(p)
  if R === S
    return local_basis_matrix(L, minimum(p), type = type)
    #if type == :any
    #  return _local_basis_matrix(pseudo_matrix(L), p)
    #elseif type == :submodule
    #  return _local_basis_submodule_matrix(pseudo_matrix(L), p)
    #elseif type == :supermodule
    #  return _local_basis_supermodule_matrix(pseudo_matrix(L), p)
    #else
    #  throw(error("""Invalid :type keyword :$(type).
    #                 Must be either :any, :submodule, or :supermodule"""))
    #end
  elseif S === base_ring(R)
    if type == :any
      return _local_basis_matrix_prime_below(pseudo_matrix(L), p)
    elseif type == :submodule
      return _local_basis_matrix_prime_below_submodule(pseudo_matrix(L), p)
    elseif type == :supermodule
      throw(NotImplemented())
    else
      error("""Invalid :type keyword :$(type).
               Must be either :any, :submodule, or :supermodule""")
    end
  else
    error("Something wrong")
  end
end

################################################################################
#
#  Jordan decomposition
#
################################################################################

@doc Markdown.doc"""
    jordan_decomposition(L::AbsLat, p::NfOrdIdl)
                                -> Vector{MatElem}, Vector{MatElem}, Vector{Int}

Returns a Jordan decomposition of the completion of $L$ at $\mathfrak p$.

The return value consists of three lists $(M_i)_i$, $(G_i)_i$ and $(s_i)_i$ of
the same length $r$. The completions of the row spans of the matrices $M_i$
yield a Jordan decomposition of $L_{\mathfrak{p}}$ into modular sublattices
$L_i$ with gram matrices $G_i$ and scale of $\mathfrak{p}$-adic valuation $s_i$.
"""
jordan_decomposition(L::AbsLat, p::NfOrdIdl)

################################################################################
#
#  Local isometry
#
################################################################################

@doc Markdown.doc"""
    islocally_isometric(L::AbsLat, M::AbsLat, p::NfOrdIdl) -> Bool

Returns whether the completions of $L$ and $M$ at the prime ideal
$\mathfrak{p}$ are locally isometric.
"""
islocally_isometric(::AbsLat, ::AbsLat, ::NfOrdIdl)

################################################################################
#
#  Isotropy
#
################################################################################

@doc Markdown.doc"""
    isisotropic(L::AbsLat, p) -> Bool

Returns whether the completion of $L$ at $p$ is isotropic.
"""
isisotropic(L::AbsLat, p) = isisotropic(rational_span(L), p)

################################################################################
#
#  Restrict scalars
#
################################################################################

function restrict_scalars(L::AbsLat)
  V = ambient_space(L)
  Vabs, f = restrict_scalars(V, FlintQQ)
  Babs = absolute_basis(L)
  Mabs = zero_matrix(FlintQQ, length(Babs), rank(Vabs))
  for i in 1:length(Babs)
    v = f\(Babs[i])
    for j in 1:length(v)
      Mabs[i, j] = v[j]
    end
  end
  return ZLat(Vabs, Mabs)
end

################################################################################
#
#  Automorphism group
#
################################################################################

# Determine the gram matrices of the bilinear forms
# V x V -> Q, (x, y) -> Tr_K/Q(a*B(x, y))
# with respect to an absolute basis of L, for all a in generators.
function Zforms(L::AbsLat{<: NumField}, generators)
  return _Zforms(L, generators)
end

function Zforms(L::AbsLat{<: NumField})
  E = base_ring(ambient_space(L))
  if degree(E) > 1
    generators = elem_type(E)[E(1), absolute_primitive_element(E)]
  else
    generators = elem_type(E)[E(1)]
  end
  return _Zforms(L, generators)
end

function _Zforms(L::AbsLat{<: NumField}, generators::Vector)
  V = ambient_space(L)
  E = base_ring(V)
  Babs = absolute_basis(L)
  Babsmat = matrix(E, Babs)
  forms = fmpz_mat[]
  scalars = fmpq[]
  for b in generators
    Vres, VresToV = restrict_scalars(V, FlintQQ, b)
    G = gram_matrix(Vres, map(t -> preimage(VresToV, t), Babs))
    d = denominator(G)
    Gint = change_base_ring(FlintZZ, d * G)
    c = content(Gint)
    G = divexact(Gint, c)
    push!(forms, G)
    push!(scalars, d//c)
  end
  return forms, scalars, Babsmat, generators
end

# Compute the automorphism group of the lattice
# per default, the are given with respect to the basis of the ambient space
# if ambient_representation = true, they are given with respect to the coordinate
# space/ambient space
function assert_has_automorphisms(L::AbsLat{<: NumField}; redo::Bool = false)

  if !redo && isdefined(L, :automorphism_group_generators)
    return nothing
  end

  E = base_ring(ambient_space(L))

  ZgramL, scalarsL, BabsmatL, generatorsL = Zforms(L)
  @assert isone(generatorsL[1])

  # So the first one is either positive definite or negative definite
  # Make it positive definite. This does not change the automorphisms.
  if ZgramL[1][1, 1] < 0
    ZgramL[1] = -ZgramL[1]
  end

  # Make the Gram matrix small
  Glll, T = lll_gram_with_transform(ZgramL[1])
  Ttr = transpose(T)
  ZgramLorig = ZgramL
  ZgramL = copy(ZgramL)
  for i in 1:length(ZgramL)
    ZgramL[i] = T * ZgramL[i] * Ttr
  end

  # Create the automorphism context and compute generators as well as orders

  C = ZLatAutoCtx(ZgramL)
  fl, Csmall = try_init_small(C)
  if fl
    auto(Csmall)
    _gens, order = _get_generators(Csmall)
    gens = fmpz_mat[matrix(ZZ, g) for g in _gens]
  else
    init(C)
    auto(C)
    gens, order = _get_generators(C)
  end

  @hassert :Lattice 1 begin
    flag = true
    for g in gens
      gt = g'
      for i in 1:length(ZgramL)
        if g * ZgramL[i] * g' != ZgramL[i]
          flag = false
        end
      end
    end
    flag
  end

  # Now undo the LLL transformation
  Tinv = inv(T)
  for i in 1:length(gens)
    gens[i] = Tinv * gens[i] * T
  end

  # Now gens are with respect to the absolute basis of L
  @hassert :Lattice 1 begin
    flag = true
    for j in 1:length(ZgramL)
      for i in 1:length(gens)
        if gens[i] * ZgramLorig[j] * transpose(gens[i]) != ZgramLorig[j]
          flag = false
        end
      end
    end
    flag
  end

  # Now translate to get the automorphisms with respect to basis_matrix(L)
  BmatL = basis_matrix_of_rational_span(L)

  b1, s1 = can_solve_with_solution(BabsmatL, BmatL, side = :left)
  b2, s2 = can_solve_with_solution(BmatL, BabsmatL, side = :left)

  t_gens = Vector{typeof(BmatL)}(undef, length(gens))

  for i in 1:length(gens)
    t_gens[i] = s1 * change_base_ring(E, gens[i]) * s2
  end

  G = gram_matrix_of_rational_span(L)
  @hassert :Lattice 1 all(g * G * _map(transpose(g), involution(L)) == G
                            for g in t_gens)

  pm = pseudo_matrix(L)
  C = coefficient_ideals(pm)

  for g in t_gens
    @hassert :Lattice 1 all(g[i, j] in C[j] * inv(C[i])
                              for i in 1:nrows(g), j in 1:nrows(g))
  end

  # Now set the generators and the order

  L.automorphism_group_generators = t_gens
  L.automorphism_group_order = order
  return nothing
end

################################################################################
#
#  Automorphism group generators
#
################################################################################

@doc Markdown.doc"""
    automorphism_group_generators(L::AbsLat; ambient_representation = true)

Given a definite lattice $L$ returns generators for the automorphism group of $L$.
If `ambient_representation` is `true` (the default), the transformations are represented
with respect to the ambient space of $L$. Otherwise, the transformations are represented
with respect to the (pseudo-)basis of $L$.
"""
automorphism_group_generators(L::AbsLat; ambient_representation::Bool = true)

function automorphism_group_generators(L::AbsLat; ambient_representation::Bool = true)

  assert_has_automorphisms(L)

  gens = L.automorphism_group_generators

  if !ambient_representation
    if check
      Grel = gram_matrix(rational_span(L))
      for g in gens
        @assert g * Grel * _map(transpose(g), involution(L)) == Grel
      end
    end
    return copy(gens)
  else
    bm = basis_matrix_of_rational_span(L)
    bminv = inv(bm)
    gens = typeof(bm)[bminv * g * bm for g in gens]
    @hassert :Lattice 1 begin
      flag = true
      Gamb = gram_matrix(ambient_space(L))
      for g in gens
        if g * Gamb * _map(transpose(g), involution(L)) != Gamb
          flag = false
        end
      end
      flag
    end
    return gens
  end
end

################################################################################
#
#  Automorphism group order
#
################################################################################

@doc Markdown.doc"""
    automorphism_group_order(L::AbsLat)

Given a definite lattice $L$ return the order of the automorphism group of $L$.
"""
automorphism_group_order(L::AbsLat; redo::Bool = false)

function automorphism_group_order(L::AbsLat; redo::Bool = false)
  assert_has_automorphisms(L, redo = redo)
  return L.automorphism_group_order
end

################################################################################
#
#  Isometry
#
################################################################################

@doc Markdown.doc"""
    isisometric(L::ZLat, M::ZLat; ambient_representation::Bool = true)
                                                              -> (Bool, MatElem)

Tests if $L$ and $M$ are isometric. If this is the case, the second return value
is an isometry $T$ from $L$ to $M$.

By default, that isometry is represented with respect to the bases of the
ambient spaces, that is, $T V_M T^t = V_L$ where $V_L$ and $V_M$ are the gram
matrices of the ambient spaces of $L$ and $M$ respectively. If
`ambient_representation = false`, then the isometry is represented with respect
to the (pseudo-)bases of $L$ and $M$, that is, $T G_M T^t = G_L$ where $G_M$
and $G_L$ are the gram matrices of the (pseudo-)bases of $L$ and $M$
respectively.
"""
isisometric(L::AbsLat, M::AbsLat; ambient_representation::Bool = true)

function isisometric(L::AbsLat{<: NumField}, M::AbsLat{<: NumField};
                                            ambient_representation::Bool = true)
  V = ambient_space(L)
  W = ambient_space(M)
  E = base_ring(V)
  K = base_field(E)
  @assert base_ring(V) == base_ring(W)
  @assert base_ring(L) == base_ring(M)

  ZgramL, scalarsL, BabsmatL, generatorsL = Zforms(L)
  ZgramM, scalarsM, BabsmatM, generatorsM = Zforms(M, generatorsL)
  @assert generatorsL == generatorsM
  if scalarsL != scalarsM
    return false, zero_matrix(E, 0, 0)
  end

  # So the first one is either positive definite or negative definite
  # Make it positive definite. This does not change the automorphisms.
  if ZgramL[1][1, 1] < 0
    ZgramL[1] = -ZgramL[1]
    ZgramM[1] = -ZgramM[1]
  end

  ZgramLsmall = copy(ZgramL)
  ZgramMsmall = copy(ZgramM)

  # Make the Gram matrix small
  _, TL = lll_gram_with_transform(ZgramL[1])
  _, TM = lll_gram_with_transform(ZgramM[1])
  TLtr = transpose(TL)
  TMtr = transpose(TM)
  for i in 1:length(ZgramL)
    ZgramLsmall[i] = TL * ZgramL[i] * TLtr
    ZgramMsmall[i] = TM * ZgramM[i] * TMtr
  end

  fl, CLsmall, CMsmall = _try_iso_setup_small(ZgramLsmall, ZgramMsmall)
  if fl
    b, _T = isometry(CLsmall, CMsmall)
    T = matrix(FlintZZ, _T)
  else
    CL, CM = _iso_setup(ZgramLsmall, ZgramMsmall)
    b, T = isometry(CL, CM)
  end

  if b
    T = change_base_ring(FlintQQ, inv(TL)*T*TM)
    fl, s1 = can_solve_with_solution(BabsmatL, basis_matrix_of_rational_span(L), side = :left)
    fl, s2 = can_solve_with_solution(basis_matrix_of_rational_span(M), BabsmatM, side = :left)
    T = s1 * change_base_ring(E, T) * s2
    @hassert :Lattice 1 T * gram_matrix(rational_span(M)) *
                            _map(transpose(T), involution(L)) ==
                                gram_matrix(rational_span(L))
    if !ambient_representation
      return true, T
    else
      T = inv(basis_matrix_of_rational_span(L)) * T *
                 basis_matrix_of_rational_span(M)

      @hassert :Lattice 1 T * gram_matrix(ambient_space(M)) *
                              _map(transpose(T), involution(L)) ==
                                  gram_matrix(ambient_space(L))
      return true, T
    end
  else
    return false, zero_matrix(E, 0, 0)
  end
end

################################################################################
#
#  Maximal sublattices
#
################################################################################

function maximal_sublattices(L::AbsLat, p; use_auto::Bool = false,
                                           callback = false, max = inf)
  @req base_ring(L) == order(p) "asdsd"

  B = local_basis_matrix(L, p, type = :submodule)
  n = nrows(B)
  R = base_ring(L)
  K = nf(R)
  k, h = ResidueField(R, p)
  hext = extend(h, K)

  if use_auto
    throw(NotImplemented())
  end

  Ls = maximal_subspaces(k, n)
  pML = _module_scale_ideal(pseudo_matrix(L), p)
  result = typeof(L)[]
  keep = true
  cont = true
  E = Int[]
  for i in 1:length(Ls)
    m = map_entries(y -> hext\y, Ls[i])
    LL = lattice(ambient_space(L), _sum_modules(L, pseudo_matrix(m * B), pML))
    if !(callback isa Bool)
      keep, cont = callback(result, LL)::Tuple{Bool, Bool}
    end
    if keep
      push!(result, LL)
      push!(E, 1)
    end
    if !cont
      break
    end
    if length(result) >= max
      break
    end
  end
  return result, E
end

################################################################################
#
#  Minimal superlattices
#
################################################################################

function minimal_superlattices(L::AbsLat, p; use_auto::Bool = false,
                                             callback = false, max = inf)
  @req base_ring(L) == order(p) "asdsd"

  B = local_basis_matrix(L, p, type = :submodule)
  n = nrows(B)
  R = base_ring(L)
  K = nf(R)
  k, h = ResidueField(R, p)
  hext = extend(h, K)

  if use_auto
    throw(NotImplemented())
  end

  Ls = enumerate_lines(k, n)

  pinv = inv(p)
  ML = pseudo_matrix(L)
  result = typeof(L)[]
  keep = true
  cont = true
  E = Int[]
  for v in Ls
    # don't need to make a copy
    m = matrix(K, 1, n, map(y -> hext\y, v))
    ppm = pseudo_matrix(m*B, [pinv])
    LL = lattice(ambient_space(L), _sum_modules(L, ML, ppm))
    if !(callback isa Bool)
      keep, cont = callback(result, LL)
    end
    if keep
      push!(result, LL)
      push!(E, 1)
    end
    if !cont
      break
    end
    if length(result) >= max
      break
    end
  end
  return result, E
end

################################################################################
#
#  Orthogonal sum
#
################################################################################

# TODO: Make this a proper coproduct with injections?
function orthogonal_sum(M::T, N::T) where {T <: AbsLat}
  @req base_ring(M) === base_ring(N) "Base rings must be equal"
  U = ambient_space(M)
  V = ambient_space(N)
  W, f1, f2 = orthogonal_sum(U, V)
  rU = rank(U)
  rV = rank(V)
  rW = rank(W)
  pM = pseudo_matrix(M)
  pN = pseudo_matrix(N)
  MpM = matrix(pM)
  MpN = matrix(pN)
  H = pseudo_matrix(diagonal_matrix(MpM, MpN),
                    vcat(coefficient_ideals(pM), coefficient_ideals(pN)))
  return lattice(W, H), f1, f2
end


text:
"""
**AbstractMolecularContainer** is the parent-type for
all other container types.
"""
abstract AbstractMolecularContainer

ext: jl


text:
mutable struct CmdProgram <: Program
	name::String
	id_file::String
	info_before::String
	info_after::String
    cmd_dependencies::Vector{CmdDependency}
	inputs::Vector{String}
	input_types::Vector{DataType}
	default_inputs::Vector
	validate_inputs::Function
	prerequisites::Function
    cmd::Base.AbstractCmd
	infer_outputs::Function
	outputs::Vector{String}
	output_types::Vector{DataType}
	default_outputs::Vector
	validate_outputs::Function
	wrap_up::Function
end

"""
	CmdProgram <: Program

	CmdProgram(;
		name::String               = "Command Program",
		id_file::String            = "",
		info_before::String        = "auto",
		info_after::String         = "auto",
		cmd_dependencies::Vector{CmdDependency} = Vector{CmdDependency}(),
		inputs                     = Vector{String}(),
		validate_inputs::Function  = do_nothing,  # positional arguments: inputs::Dict{String, ValidInputTypes}
		prerequisites::Function    = do_nothing,  # positional arguments: inputs, outputs::Dict{String, ValidInputTypes}
		cmd::Base.AbstractCmd      = ``,
		infer_outputs::Function    = do_nothing,  # positional arguments: inputs::Dict{String, ValidInputTypes}
		outputs                    = Vector{String}(),
		validate_outputs::Function = do_nothing  # positional arguments: outputs::Dict{String, ValidInputTypes},
		wrap_up::Function          = do_nothing  # positional arguments: inputs, outputs::Dict{String, ValidInputTypes}
	) -> CmdProgram

Command program template. To run a `CmdProgram`, use `run(::CmdProgram; kwargs...).`

# Arguments

- `name::String`: Program name.

- `id_file::String`: The prefix of *run ID file*. To prevent from running the program with the same inputs and outputs twice, it will generate a unique *run ID file* after a successful run.

- `info_before::String`: Print it when starting the program.

- `info_after::String`: Print it when finishing the program.

- `cmd_dependencies::Vector{CmdDependency}`: Any command dependencies used in the program.

- `inputs` and `outputs`: Elements (or vectors containing elements) in the following format: (1) `keyword` (2) `keyword => data_type` (3) `keyword => default_value` (4) `keyword => default_value => data_type`.

  `keyword` is an argument name, can be `String` or `Symbol`.

  `default_value` is optional. If set, users may not provide this argument when running. Elsewise, users have to provide it. Caution: `nothing` is preserved and means default value not set. If `String`, it can contain other keywords, but need to quote using '<>', such as `"<arg>.txt"`

  `data_type` is optional. If set, the value provided have to be this data type, or an error will throw.

  *HOW DOES THIS WORK?*

  > `CmdProgram` stores a command template. In the template, replaceable portions are occupied by *keywords*, and all keywords are set in `inputs` and `outputs`.
  > `keyword`s will be replaced before running the program. Users need to provide a dictionary of `keyword::String => value` in `run(::Program, inputs::Dict{String}, outputs::Dict{String})`.

- `validate_inputs::Function`: A function to validate inputs. It takes *one* argument `Dict{String, ValidInputTypes}` whose keys are the same as `inputs`. If validation fail, throw error or return false.

- `prerequisites`: A function to run just before the main command. It prepares necessary things, such as creating directories. It takes *two* arguments `Dict{String, ValidInputTypes}` whose keys are the same as `inputs` and `outputs`, respectively.

- `cmd::AbstractCmd`: The main command template. In the template, keywords in `inputs::Vector{String}` and `outputs::Vector{String}` will be replaced when envoking `run(::CmdProgram, inputs::Dict{String, ValidInputTypes}, outputs::Dict{String, ValidInputTypes})`.

- `infer_outputs::Function`: A function to infer outputs from inputs. It takes *one* argument `Dict{String, ValidInputTypes}` whose keys are the same as `inputs`.

- `validate_outputs::Function`: A function to validate outputs. It takes *one* argument `Dict{String, ValidInputTypes}` whose keys are the same as `outputs`. If validation fail, throw error or return false.

- `wrap_up::Function`: the last function to run. It takes *two* arguments `Dict{String, ValidInputTypes}` whose keys are the same as `inputs` and `outputs`, respectively.

# Example

	p = CmdProgram(
		id_file = "id_file",
		inputs = [
			"input",
			"input2" => Int,
			"optional_arg" => 5,
			"optional_arg2" => 0.5 => Number
		],
		outputs =
			"output" => "<input>.output"
		,
		cmd = `echo input input2 optional_arg optional_arg2 output`
	)

	inputs = Dict(
		"input" => `in1`,
		"input2" => 2
	)

	outputs = Dict(
		"output" => "out"
	)

	run(p, inputs, outputs;
		touch_run_id_file = false
	)
"""
function CmdProgram(;
	name::String               = "Command Program",
	id_file::String            = "",
	info_before::String        = "auto",
	info_after::String         = "auto",
	cmd_dependencies::Vector{CmdDependency} = Vector{CmdDependency}(),
	inputs                     = Vector{String}(),
	validate_inputs::Function  = do_nothing,  # positional arguments: inputs::Dict{String}
	prerequisites::Function    = do_nothing,  # positional arguments: inputs, outputs::Dict{String}
	cmd::Base.AbstractCmd      = ``,
	infer_outputs::Function    = do_nothing,  # positional arguments: inputs::Dict{String}
	outputs                    = Vector{String}(),
	validate_outputs::Function = do_nothing,  # positional arguments: outputs::Dict{String}
	wrap_up::Function          = do_nothing  # positional arguments: inputs, outputs::Dict{String}
)
	inputs, input_types, default_inputs = parse_default(inputs)
	outputs, output_types, default_outputs = parse_default(outputs)

	CmdProgram(
		name,
		id_file,
		info_before,
		info_after,
		cmd_dependencies,
		inputs,
		input_types,
		default_inputs,
		validate_inputs,
		prerequisites,
		cmd,
		infer_outputs,
		outputs,
		output_types,
		default_outputs,
		validate_outputs,
		wrap_up
	)
end

"""
	run(
		p::CmdProgram;
		inputs=Dict{String}(),
		outputs=Dict{String}(),
		dir::AbstractString="",
		check_dependencies::Bool=true,
		skip_when_done::Bool=true,
		touch_run_id_file::Bool=true,
		verbose::Bool=true,
		dry_run::Bool=false,
		stdout=nothing,
		stderr=nothing,
		stdlog=nothing,
		append::Bool=false
	) -> (success::Bool, outputs::Dict{String})

	run(p::CmdProgram, inputs, outputs; kwargs...)

	run(p::CmdProgram, inputs; kwargs...)
	)  # only usable when `p.infer_outputs` is defined, or default outputs are set in `p`.

Run Command Program (CmdProgram) using specified `inputs` and `outputs`.

Return `(success::Bool, outputs::Dict{String})`

- `p::CmdProgram`: the command program template.

- `inputs::Dict{String}` and `outputs::Dict{String}`: `p::CmdProgram` stores a command template. In the template, replaceable portions are occupied by *keywords*, and all keywords can be found at `p.inputs` and `p.outputs` string vectors. Here, `inputs` and `outputs` are `Dict(keyword::String => replacement)`. The replacements do not have a length limit, unless a *keyword* refers to a filename (length == 1).

  > If data types of `inputs` and `outputs` are not `Dict{String}`, they will be converted as far as possible. If the conversion fails, program will throw an error.

- `dir::AbstractString = ""`: working directory to run the program.

- `check_dependencies::Bool = true`: check dependencies for `p` (`p.cmd_dependencies`).

- `skip_when_done::Bool = true`: Skip running the program and return `true` if it has been done before (the `run_id_file` exists and `p.validate_outputs(outputs)` passes.)

- `touch_run_id_file::Bool = true`: If `true`, touch a unique run ID file, which indicate the program is successfully run with given inputs and outputs. If `false`, the next time running the program, `skip_when_done=true` will not take effect.

- `verbose::Bool = true`: If `true`, print info and error messages. If `false`, print error messages only.

- `dry_run::Bool = false`: do not run the command, return `(command::AbstractCmd, run_id_file::String)`.

- `stdout`, `stderr`, `stdlog` and `append`: Redirect the program outputs to files. `stdlog` is the Julia logging of `@info`, `@warn`, `@error`, etc. Caution: If the original command (`p.cmd`) has redirection, arguments defined here might not be effective for the command.

!!! note
    Redirecting in Julia are not thread safe, so unexpected redirection might be happen if you are running programs in different `Tasks` or multi-thread mode.

### Workflow

1. Go to the working directory. Establish redirection. (`dir`, `stdout`, `stderr`, `stdlog`, `append`).

2. Validate compatibility between `p` and `inputs/outputs`.

3. Check whether the program has run before. (`skip_when_done`, `p.validate_outputs(outputs)`)

4. Check command dependencies. (`check_dependencies`, `p.cmd_dependencies`)

5. Validate `inputs`. (`p.validate_inputs(inputs)`)

6. Generate runnable command from `p` and `inputs/outputs`. (`stdout`, `stderr`, `append`)

7. Preparing before running main command. (`p.prerequisites(inputs, outputs)`)

8. Run command generated in #5.

9. Validate `outputs`. (`p.validate_outputs(outputs)`)

10. Wrap up. (`p.wrap_up(inputs, outputs)`)

11. Success, touch run id file, and return `(success::Bool, outputs::Dict{String})`. (`touch_run_id_file::Bool`)

# Example

	p = CmdProgram(
		id_file = "id_file",
		inputs = ["input", "input2"],
		outputs = ["output"],
		cmd = `echo input input2 output`
	)

	inputs = Dict(
		"input" => `in1`,
		"input2" => 2
	)

	outputs = Dict(
		"output" => "out"
	)

	run(p, inputs, outputs;
		touch_run_id_file = false
	)
"""
function _run(
	p::CmdProgram;
	inputs=Dict{String, Any}(),
	outputs=Dict{String, Any}(),
	skip_when_done::Bool=true,
	check_dependencies::Bool=true,
	# stdout=nothing,
	# stderr=nothing,
	# append::Bool=false,
	verbose::Bool=true,
	touch_run_id_file::Bool=true,
	dry_run::Bool=false
)
	# input/output completion
	inputs, outputs = xxputs_completion_and_check(p, inputs, outputs)

	# run id based on inputs and outputs
	run_id = generate_run_uuid(inputs, outputs)
	run_id_file = p.id_file * "." * string(run_id)

	# whether dry run
	if dry_run
		@goto dry_run_start
	end

	if verbose
		if p.info_before == "auto" || p.info_before == ""
			@info timestamp() * "Starting program: $(p.name)" command_template=p.cmd run_id inputs outputs
		else
			@info timestamp() * p.info_before command_template=p.cmd run_id inputs outputs
		end
	end

	# skip when done
	if skip_when_done && isfile(run_id_file) && isok(p.validate_outputs(outputs))
		@warn timestamp() * "Skipping finished program: $(p.name)" command_template=p.cmd run_id inputs outputs
		return true, outputs
	end

	# check dependencies
	if check_dependencies
		foreach(check_dependency, p.cmd_dependencies)
	end

	# preparation: remove run id file
	isfile(run_id_file) && rm(run_id_file)

	# preparation: validate inputs
	try
		isok(p.validate_inputs(inputs)) || error("ProgramInputValidationError: $(p.name): the prerequisites function returns false.")
	catch e
		@error timestamp() * "ProgramInputValidationError: $(p.name): fail to validate inputs (before running the main command). See error messages above." validation_function=p.validate_inputs command_template=p.cmd run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramInputValidationError")
		return false, outputs
	end

	@label dry_run_start
	# preparation: replace inputs and outputs in cmd, including redirecting files
	cmd = prepare_cmd(p, inputs, outputs)

	if dry_run
		return (cmd, run_id_file)
		# dry_run stops
	end

	# preparation: run function prerequisites before the main command
	try
		isok(p.prerequisites(inputs, outputs)) || error("ProgramPrerequisitesError: $(p.name): the prerequisites function returns false.")
	catch e
		@error timestamp() * "ProgramPrerequisitesError: $(p.name): fail to run the prerequisites function (before running the main command). See error messages above." prerequisites=p.prerequisites command_template=p.cmd run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramPrerequisitesError")
		return false, outputs
	end

	# run the main command
	try
		run(cmd)
		# run(pipeline(cmd, stdout=stdout, stderr=stderr, append=append))
	catch e
		@error timestamp() * "ProgramRunningError: $(p.name): fail to run the main command. See error messages above." prerequisites=p.prerequisites command_running=cmd run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramRunningError")
		return false, outputs
	end

	# confirmation: validate outputs
	try
		isok(p.validate_outputs(outputs)) || error("ProgramOutputValidationError: $(p.name): the validation function returns false.")
	catch e
		@error timestamp() * "ProgramOutputValidationError: $(p.name): fail to validate outputs (after running the main command). See error messages above." validation_function=p.validate_outputs command_running=cmd run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramOutputValidationError")
		return false, outputs
	end

	try
		isok(p.wrap_up(inputs, outputs)) || error("ProgramWrapUpError: $(p.name): the wrap_up function returns false.")
	catch e
		@error timestamp() * "ProgramWrapUpError: $(p.name): fail to run the wrap_up function. See error messages above." wrap_up=p.wrap_up command_running=cmd run_id inputs outputs exception=(e, catch_backtrace())
		error("ProgramWrapUpError")
		return false, outputs
	end

	# touch the run_id_file
	touch_run_id_file && touch(run_id_file)

	if verbose
		if p.info_after == "auto" || p.info_after == ""
			@info timestamp() * "Finishing program: $(p.name)" command_running=cmd run_id inputs outputs
		else
			@info timestamp() * p.info_after command_running=cmd run_id inputs outputs
		end
	end
	return true, outputs
end

function prepare_cmd(c::Cmd, inputs::Dict{String}, outputs::Dict{String})
	cmd = deepcopy(c)
	n = length(cmd.exec)
	for i = n:-1:1
		replacement = get(inputs, cmd.exec[i], nothing) |> to_cmd
		if isnothing(replacement)
			replacement = get(outputs, cmd.exec[i], nothing) |> to_cmd
			if !isnothing(replacement)
				splice!(cmd.exec, i, replacement)
			end
		else
			splice!(cmd.exec, i, replacement)
		end
	end
	cmd
end

function prepare_cmd(c::Base.CmdRedirect, inputs::Dict{String}, outputs::Dict{String})
	Base.CmdRedirect(
		prepare_cmd(c.cmd, inputs, outputs),
		prepare_cmd(c.handle, inputs, outputs),
		c.stream_no,
		c.readable
	)
end
function prepare_cmd(c::T, inputs::Dict{String}, outputs::Dict{String}) where T <: Union{Base.OrCmds, Base.ErrOrCmds, Base.AndCmds}
	T(
		prepare_cmd(c.a, inputs, outputs),
		prepare_cmd(c.b, inputs, outputs)
	)
end
function prepare_cmd(h::Base.TTY, inputs::Dict{String}, outputs::Dict{String})  # h: handle property of CmdRedirect
	h
end
function prepare_cmd(h::Base.FileRedirect, inputs::Dict{String}, outputs::Dict{String})  # h: handle property of CmdRedirect
	replacement = get(outputs, h.filename, nothing) |> to_cmd
	if isnothing(replacement)
		replacement = get(inputs, h.filename, nothing) |> to_cmd
		if !isnothing(replacement)
			@goto to_replace
		end
	else
		@goto to_replace
	end
	return h  # nothing change

	@label to_replace
	if length(replacement.exec) == 1
		return Base.FileRedirect(replacement.exec[1], h.append)
	else
		error("ProgramRedirectError: a key ($(h.filename)) in inputs or outputs found in redirecting filename, but the provided value contain multiple elements. Please only provide one element. For example, quote the filename when contain special characters.")
	end
end

"""
	prepare_cmd(p::CmdProgram, inputs, outputs)

Prepare the runable command. Keywords in CmdProgram will be given to values of inputs/outputs.
"""
function prepare_cmd(p::CmdProgram, inputs::Dict{String}, outputs::Dict{String})
	prepare_cmd(p.cmd, inputs, outputs)
end

function prepare_cmd(p::CmdProgram, inputs, outputs)
	prepare_cmd(p.cmd, to_xxput_dict(inputs), to_xxput_dict(outputs))
end


text:
@testset "Subgroup enumeration" begin
  @testset "p-subgroups" begin
    G = abelian_group([9, 3, 3])
    H = abelian_group([9, 3, 3, 5, 7])

    GG = abelian_group([8])
    @test 4 == length(collect(subgroups(GG)))

    @testset "All subgroups" begin
      @test 50 == length(collect(psubgroups(G, 3)))
    end

    @testset "Given subgroup type" begin

      T = psubgroups(G, 3, subtype = [3])
      @test 0 == length(collect(T))

      T = psubgroups(G, 3, subtype = [1, 1, 1, 1])
      @test 0 == length(collect(T))

      T = psubgroups(G, 3, subtype = [2])
      @test 9 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[9] for t in T])

      T = psubgroups(G, 3, subtype = [2])
      @test 9 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[9] for t in T])

      T = psubgroups(G, 3, subtype = [1, 1])
      @test 13 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[3, 3] for t in T])

      T = psubgroups(G, 3, subtype = Int[])
      @test 1 == length(collect(T))
      @test all([order(t[1]) == 1 for t in T])

      T = psubgroups(G, 3, subtype = Int[2, 1, 1])
      @test 1 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[3, 3, 9] for t in T])

      TH = psubgroups(G, 3, subtype = [3])
      @test 0 == length(collect(TH))

      TH = psubgroups(G, 3, subtype = [1, 1, 1, 1])
      @test 0 == length(collect(TH))

      TH = psubgroups(G, 3, subtype = [2])
      @test 9 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[9] for t in TH])

      TH = psubgroups(G, 3, subtype = [2])
      @test 9 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[9] for t in TH])

      TH = psubgroups(G, 3, subtype = [1, 1])
      @test 13 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[3, 3] for t in TH])

      TH = psubgroups(G, 3, subtype = Int[])
      @test 1 == length(collect(TH))
      @test all([order(t[1]) == 1 for t in TH])

      TH = psubgroups(G, 3, subtype = Int[2, 1, 1])
      @test 1 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[3, 3, 9] for t in TH])
    end

    @testset "Given quotient type" begin
      T = psubgroups(G, 3, quotype = [2], fun = quo)
      @test 9 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[9] for t in T])

      T = psubgroups(G, 3, quotype = [1, 1], fun = quo)
      @test 13 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[3, 3] for t in T])

      T = psubgroups(G, 3, quotype = [2, 1, 1], fun = quo)
      @test 1 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[3, 3, 9] for t in T])

      T = psubgroups(G, 3, quotype = Int[], fun = quo)
      @test 1 == length(collect(T))
      @test all([order(t[1]) == 1 for t in T])

      TH = psubgroups(G, 3, quotype = [2], fun = quo)
      @test 9 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[9] for t in TH])

      TH = psubgroups(G, 3, quotype = [1, 1], fun = quo)
      @test 13 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[3, 3] for t in TH])

      TH = psubgroups(G, 3, quotype = [2, 1, 1], fun = quo)
      @test 1 == length(collect(TH))
      @test all([snf(t[1])[1].snf == fmpz[3, 3, 9] for t in TH])

      TH = psubgroups(G, 3, quotype = Int[], fun = quo)
      @test 1 == length(collect(TH))
      @test all([order(t[1]) == 1 for t in TH])
    end

    @testset "Given order" begin

      T = psubgroups(G, 3, order = 1)
      @test 1 == length(collect(T))
      @test all([order(t[1]) == 1 for t in T])

      T = psubgroups(G, 3, order = 3)
      @test 13 == length(collect(T))
      @test all([order(t[1]) == 3 for t in T])

      T = psubgroups(G, 3, order = 9)
      @test 22 == length(collect(T))
      @test all([order(t[1]) == 9 for t in T])

      T = psubgroups(G, 3, order = 81)
      @test 1 == length(collect(T))
      @test all([order(t[1]) == 81 for t in T])

      TH = psubgroups(G, 3, order = 1)
      @test 1 == length(collect(TH))
      @test all([order(t[1]) == 1 for t in TH])

      TH = psubgroups(G, 3, order = 3)
      @test 13 == length(collect(TH))
      @test all([order(t[1]) == 3 for t in TH])

      TH = psubgroups(G, 3, order = 9)
      @test 22 == length(collect(TH))
      @test all([order(t[1]) == 9 for t in TH])

      TH = psubgroups(G, 3, order = 81)
      @test 1 == length(collect(TH))
      @test all([order(t[1]) == 81 for t in TH])
    end

    @testset "Given index" begin

      T = psubgroups(G, 3, index = 1, fun = quo)
      @test 1 == length(collect(T))
      @test all([order(t[1]) == 1 for t in T])

      T = psubgroups(G, 3, index = 3, fun = quo)
      @test 13 == length(collect(T))
      @test all([order(t[1]) == 3 for t in T])

      T = psubgroups(G, 3, index = 9, fun = quo)
      @test 22 == length(collect(T))
      @test all([order(t[1]) == 9 for t in T])

      T = psubgroups(G, 3, index = 81, fun = quo)
      @test 1 == length(collect(T))
      @test all([order(t[1]) == 81 for t in T])

      TH = psubgroups(G, 3, index = 1, fun = quo)
      @test 1 == length(collect(TH))
      @test all([order(t[1]) == 1 for t in TH])

      TH = psubgroups(G, 3, index = 3, fun = quo)
      @test 13 == length(collect(TH))
      @test all([order(t[1]) == 3 for t in TH])

      TH = psubgroups(G, 3, index = 9, fun = quo)
      @test 22 == length(collect(TH))
      @test all([order(t[1]) == 9 for t in TH])

      TH = psubgroups(G, 3, index = 81, fun = quo)
      @test 1 == length(collect(TH))
      @test all([order(t[1]) == 81 for t in TH])
    end
  end

  @testset "Arbitrary groups" begin
    G = abelian_group([3, 5, 7, 9, 25])

    @test_throws ErrorException subgroups(G, index = 3, subtype = [3, 1])

    @testset "All subgroups" begin
      @test 280 == length(collect(subgroups(G)))
    end

    @test_throws ErrorException subgroups(G, index = 3, order = 2)

    @testset "Given subtype" begin
      T = subgroups(G, subtype = [25, 27])
      @test length(collect(T)) == 0

      T = subgroups(G, subtype = [25, 9, 3])
      @test 5 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[3, 225] for t in T])

      T = subgroups(G, subtype = [5, 5, 9, 3])
      @test 1 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[15, 45] for t in T])

      T = subgroups(G, subtype = [3, 5, 9, 5])
      @test 1 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[15, 45] for t in T])

    end

    @testset "Given quotype" begin
      T = subgroups(G, quotype = [25, 27], fun = quo)
      @test 0 == length(collect(T))

      T = subgroups(G, quotype = [5, 7], fun = quo)
      @test 6 == length(collect(T))
      @test all([snf(t[1])[1].snf == fmpz[35] for t in T])
    end

    @testset "Given order" begin
      T = subgroups(G, order = 5*7*3)
      @test 24 == length(collect(T))
      @test all([order(t[1]) == 5*7*3 for t in T])
    end

    @testset "Given index" begin
      T = subgroups(G, index = 7*3, fun = quo)
      @test 4 == length(collect(T))
      @test all([order(t[1]) == 7*3 for t in T])
    end
  end
end


text:
include("Init.jl")
include("ProjectInvestment.jl")

##########################################################################################################
#Investment decision making algorithm"
#GEP = "cost-minimizing future investments" price projection method
#EMLab = "myopic agents" price projection method
#EMCAS = "exogenous scenarios for future investments" price projection method
investmentModelSelectionList = [:GEP,:EMLab,:EMCAS]
#Selecting the look-ahead horizon
genCoHorizonList = [5,10,15,20]

#Store the simulation history
evaluationHistory_List = Array{Any,2}(undef,length(investmentModelSelectionList),9)
systemHistory_List = Array{Any,2}(undef,length(investmentModelSelectionList),9)
allGenCos_List = Array{Any,2}(undef,length(investmentModelSelectionList),9)

#Looping through the cases tested in the study
for modelCount in 1:length(investmentModelSelectionList)
    systemSettings.simulationSettings.investmentModelSelection = investmentModelSelectionList[modelCount]
    #For "exogenous scenarios for future investments" method,
        #we test different scenario tree settings
    if modelCount == 3
        #First layer of the scenario tree
        newExpansionList = [0.95,0.9,0.85]
        #Second layer of the scenario tree
        newDivisionList = [(0.2,0.5,0.3),(0.3,0.2,0.5),(0.5,0.3,0.2)]
        caseCount = 0
        #Looping through different scenario tree settings
        #Outer loop: first layer
        for newExpansion in newExpansionList
            #Inner layer: second layer
            for newDivision in newDivisionList
                caseCount += 1
                #Create a new case and serve as input to the "Init.jl" function where GenCos are initialized
                #The function "NewCase" can be found in the file "ScenarioTree.jl"
                gapCoverPercentage,divisionTypeShare = NewCase(newExpansion,newDivision)
                #Initialize the system
                #MO = Market Operator
                #allGenCos = all Generation Companies
                #systemLoad = load data
                #allPowerFlows = all power flows (connection between market operator and generators)
                #technologyPool = all candidate technologies
                MO,allGenCos,systemLoad,allPowerFlows,technologyPool,systemState,systemHistory = InitSystem(systemSettings,genCoHorizon;gapCoverPercentage = gapCoverPercentage,divisionTypeShare = divisionTypeShare)
                #Run the simulation
                SimulateInvestment(systemSettings,MO,allGenCos,systemLoad,allPowerFlows,technologyPool,systemState,systemHistory)
                #Record the history
                systemHistory_List[modelCount,caseCount] = systemHistory
                allGenCos_List[modelCount,caseCount] = allGenCos
            end
        end
    else
        #For "cost-minimizing future investments" and "myopic agents" price projection method, we test the sensitivity of look-ahead horizon
        for genCoHorizonCount in 1:length(genCoHorizonList)
            #Similar operation as above
            MO,allGenCos,systemLoad,allPowerFlows,technologyPool,systemState,systemHistory = InitSystem(systemSettings,genCoHorizonList[genCoHorizonCount])
            SimulateInvestment(systemSettings,MO,allGenCos,systemLoad,allPowerFlows,technologyPool,systemState,systemHistory)
            systemHistory_List[modelCount,genCoHorizonCount] = systemHistory
            allGenCos_List[modelCount,genCoHorizonCount] = allGenCos
        end
    end
end


text:
# GaussQuadrature.jl
mutable struct DNq
    n::Array{Float64, 2}
    w::Array{Float64, 1}
end

"""
    initializenode(ndims::Int, quadpts::Int)
Initialize the nodes and weights of Gauss-Hermite quadrature (Multidimentional version).
Generate `ndims` dimensional quadrature nodes and weights.

# Arguments

# Optional arguments

- `dentype`
    - `:DN` Discrete normal
    - `:GH` Gauss-Hermite quadrature

# Example
```julia
ghn = initializenode(3, 21)
```
"""
function initializenode(ndims, quadpts; θrange = (-4, 4))
    nodes = range(θrange[1], stop = θrange[2], length = quadpts)
    weights = pdf.(Normal(), nodes) ./ sum(pdf.(Normal(), nodes))
    ghn = Dict(zip(1:quadpts, nodes))
    ghw = Dict(zip(1:quadpts, weights))
    # https://stackoverflow.com/questions/54256481/converting-array-of-cartesianindex-to-2d-matrix-in-julia
    as_ints(a::AbstractArray{CartesianIndex{L}}) where L = reshape(reinterpret(Int, a), (L, size(a)...))
    idx = convert(Matrix, as_ints(vec(CartesianIndices(Array{Float64}(undef, Tuple(fill(quadpts, ndims))))))')
    n = [ghn[idx[i, j]] for i in axes(idx, 1), j in axes(idx, 2)]
    w = prod([ghw[idx[i, j]] for i in axes(idx, 1), j in axes(idx, 2)], dims = 2)[:]
    return DNq(n, w)
end


text:
using SMC
using Distributions

type LinregState
    # these change
    intercept::Float64
    slope::Float64

    # these don't change they are parameters
    prior_std::Float64
    noise_std::Float64
end

# JUST USE INDEPENDENT MH PROPOSAL FOR REJUVENATION

immutable Datum
    x::Float64
    y::Float64
end

function log_likelihood(datum::Datum, state::LinregState)
    pred_y = state.intercept + state.slope * datum.x
    dist = Normal(pred_y, state.noise_std)
    logpdf(dist, datum.y)
end

immutable LinregInitializer
    datum::Datum
    prior_std::Float64
    noise_std::Float64
end

function SMC.forward(init::LinregInitializer)
    # sample the intercept and slope from the prior, weight by the likeliood
    intercept = rand(Normal(0.0, init.prior_std))
    slope = rand(Normal(0.0, init.noise_std))
    state = LinregState(intercept, slope, init.prior_std, init.noise_std)
    log_weight = log_likelihood(init.datum, state)
    (state, log_weight)
end

function SMC.backward(init::LinregInitializer, state::LinregState)
    log_likelihood(init.datum, state)
end

immutable LinregIncrementer
    data::Array{Datum,1} # this is the full list of all data, the incrementer only uses a prefix of it
    datum_index::Int
    rejuvenation_iters::Int
    function LinregIncrementer(data::Array{Datum,1}, datum_index::Int,
                               rejuvenation_iters::Int)
        @assert datum_index > 1
        new(data, datum_index, rejuvenation_iters)
    end
end

function independent_mh_step!(state::LinregState, data::Array{Datum,1})
    # propose from the prior
    intercept_old = state.intercept
    slope_old = state.slope
    old_ll = 0.0
    for datum in data
        old_ll += log_likelihood(datum, state)
    end
    state.intercept = rand(Normal(0.0, state.prior_std))
    state.slope = rand(Normal(0.0, state.prior_std))
    new_ll = 0.0
    for datum in data
        new_ll += log_likelihood(datum, state)
    end
    if log(rand()) >= new_ll - old_ll
        # reject
        state.intercept = intercept_old
        state.slope = slope_old
    end
end

function SMC.forward(incr::LinregIncrementer, state::LinregState)
    new_state = deepcopy(state)
    # do rejuvenation over already incorporate data points
    incorporated_data = incr.data[1:incr.datum_index-1]
    for i=1:incr.rejuvenation_iters
        independent_mh_step!(new_state, incorporated_data)
    end
    # incorporate new datapoint and compute log likelihood as weight
    log_weight = log_likelihood(incr.data[incr.datum_index], new_state)
    (new_state, log_weight)
end

function SMC.backward(incr::LinregIncrementer, new_state::LinregState)
    # compute log likelihood for weight, before rejuvenation
    log_weight = log_likelihood(incr.data[incr.datum_index], new_state)

    state = deepcopy(new_state)
    # do rejuvenation sweep over data
    incorporated_data = incr.data[1:incr.datum_index-1]
    for i=1:incr.rejuvenation_iters
        independent_mh_step!(state, incorporated_data)
    end
    (state, log_weight)
end

function LinregSMCScheme(x::Array{Float64,1}, y::Array{Float64,1}, prior_std::Float64, noise_std::Float64,
                         rejuvenation_iters::Int, num_particles::Int)
    @assert length(x) == length(y)
    data = Array{Datum,1}()
    for (xi, yi) in zip(x, y)
        push!(data, Datum(xi, yi))
    end
    initializer = LinregInitializer(data[1], prior_std, noise_std)
    incrementers = Array{LinregIncrementer,1}(length(data)-1)
    for i=2:length(data)
        incrementers[i-1] = LinregIncrementer(data, i, rejuvenation_iters)
    end
    SMCScheme(initializer, incrementers, num_particles)
end

function log_joint_density(state::LinregState, x::Array{Float64,1}, y::Array{Float64,1})
    log_prior = logpdf(Normal(0., state.prior_std), state.intercept)
    log_prior += logpdf(Normal(0., state.prior_std), state.slope)
    log_likelihood = 0.0
    for (xi, yi) in zip(x, y)
        y_expected = state.intercept + state.slope * xi
        log_likelihood += logpdf(Normal(y_expected, state.noise_std), yi)
    end
    log_prior + log_likelihood
end

immutable LinregSMCSampler <: Sampler
    inner_module::SMCSampler
    prior_std::Float64
    noise_std::Float64
    function LinregSMCSampler(x::Array{Float64,1}, y::Array{Float64,1}, prior_std::Float64, noise_std::Float64,
                             rejuvenation_iters::Int, num_particles::Int)
        smc_scheme = LinregSMCScheme(x, y, prior_std, noise_std, rejuvenation_iters, num_particles)
        inner_module = SMCSampler(smc_scheme, (sample) -> log_joint_density(sample, x, y))
        new(inner_module, prior_std, noise_std)
    end
end

function simulate(sampler::LinregSMCSampler)
    (sample, log_weight) = simulate(sampler.inner_module)
    ([sample.intercept, sample.slope], log_weight)
end

function regenerate(sampler::LinregSMCSampler, intercept_and_slope::Array{Float64,1})
    @assert length(intercept_and_slope) == 2
    (intercept, slope) = intercept_and_slope
    state = LinregState(intercept, slope, sampler.prior_std, sampler.noise_std)
    regenerate(sampler.inner_module, state)
end


text:
abstract type LCVMRuleModel <: MixingRule end

struct LCVMRule{γ} <: LCVMRuleModel
    components::Array{String,1}
    activity::γ
end

has_sites(::Type{<:LCVMRuleModel}) = false
has_groups(::Type{<:LCVMRuleModel}) = false
built_by_macro(::Type{<:LCVMRuleModel}) = false

function Base.show(io::IO, mime::MIME"text/plain", model::LCVMRule)
    return eosshow(io, mime, model)
end

function Base.show(io::IO, model::LCVMRule)
    return eosshow(io, model)
end

export LCVMRule
function LCVMRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    model = LCVMRule(components, init_activity)
    return model
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::LCVMRuleModel,α,a,b)
    n = sum(z)
    x = z./n
    invn2 = (one(n)/n)^2
    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
    b̄ = dot(z,Symmetric(b),z) * invn2

    ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)

    λ  = 0.7
    AV = -0.52
    AM = -0.623
    C1 = λ/AV+(1-λ)/AM

    ā = b̄*R̄*T*(C1*(g_E/(R̄*T)-0.3*sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps))+sum(x[i]*ᾱ[i,i] for i ∈ @comps))
    return ā,b̄
end

is_splittable(::LCVMRule) = true


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

name = "Xerces"
version = v"3.2.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://apache.mirror.digionline.de//xerces/c/3/sources/xerces-c-3.2.3.tar.gz", "fb96fc49b1fb892d1e64e53a6ada8accf6f0e6d30ce0937956ec68d39bd72c7e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd xerces-c-3.2.3/
atomic_patch -p1 "${WORKSPACE}/srcdir/ThreadTest.patch"
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -Dnetwork:BOOL=OFF
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/xerces-c-3.2.3/LICENSE
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libxerces-c", :libxerces)
]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@doc doc"""
Representation of a pair copula contruction. A `PCC` stores an `Array`
of copulas together with a `Vine` instance representing the actual
density decomposition. Within the `Vine` instance, the density
decomposition is stored in parent notation for reasons of memory
efficiency.

## Copulas

Involved copulas are stored as `Array` for memory efficiency. However,
upper level access functions treat the array of copulas as if they
were stored in a upper triangular matrix, with copula C_{ij|klm}
accessible through indices i and j. Thereby, i always is smaller than
j! In order to switch between triangular notation and array
implementation, function `triang2sub` converts linear triangular
matrix indexing into row and column indices, and function `sub2triang`
working the other direction.

### Simplifying assumption

In order to allow for conditional copulas, the array of copulas is of
abstract type, enabling static `PairCop` copulas as well as
`VarPairCop` as varying conditional copulas.

As a means of characterizing a pair copula construction, it suffices
to store varying conditional copulas as type `VarPairCop`, with one
field capturing the copula family, and a second field capturing an
anonymous function.

When it comes to evaluating functions on the conditional copula,
however, the conditional functional needs to be resolved. For example,
`hfun(varCop, u1, u2)` does not work, as there would be no way to
calculate the actual parameters of `varCop` inside of the function:
`varCop` might depend on u3, u4,...
"""->
type PCC{T} <: Copula
    pccCops::Array{T, 1}
    vine::Vine

    function PCC(cops::Array{T, 1}, vn::Vine)
        ## test sizes
        nVars = size(vn.trees, 2)
        nCops = sum((nVars-1):-1:1)
        if !(size(cops, 1) == nCops)
            error("Vine size must coincide with number of copulas. ")
        end
        return new(cops, vn)
    end
end

function PCC{T}(cops::Array{T, 1}, vn::Vine)
    return PCC{T}(cops, vn)
end

function PCC(vn::Vine)
    # fill with independence copulas
    nVars = size(vn.trees, 2)
    nCops = sum((nVars-1):-1:1)
    cops = ParamPC_Cpp[IndepPC_Cpp([0.0]) for ii=1:nCops]
    return PCC(cops, vn)
end

## get information
##----------------

function dim(pcc::PCC)
    return size(pcc.vine.trees, 1)
end

#############
## display ##
#############

import Base.Multimedia.display
@doc doc"""
For each pair copula, get its name in terms of conditioning sets
(1,2|3,6 for example) and the respective copula family that is used.
Display pair copulas sorted with increasing conditioning set.
"""->
function display(pcc::PCC)
    # get all conditioning sets
    condSets = Copulas.getAllCopCondSets(pcc.vine)
    nVars = dim(pcc)
    typ = eltype(pcc.pccCops)
    println("$nVars dimensional pair copula construction of type $typ:")

    # get maximum copula name length
    maxNameLen = []
    if nVars < 10
        maxNameLen = nVars + (nVars-1)
    else
        maxNameLen = 9 + (nVars-9)*2 + (nVars-1)
    end

    maxCondSetLen = nVars - 2
    for setLen=0:maxCondSetLen # for each layer
        # get conditioning sets of this layer
        inds = find([length(cs) == setLen for cs in condSets])
        # print name of each copula
        for ii=1:length(inds)
            copName = getCopNam(inds[ii], nVars, condSets)
            copName = rpad(copName, maxNameLen, " ")
            cop = pcc.pccCops[inds[ii]]
            copFam = displayShortString(cop)
            println("$copName : $copFam")
        end
    end
end

@doc doc"""
For a PCC-copula given by linear triangular indexing return its
respective name as string (for example: 1,3|4,5).
"""->
function getCopNam(ind::Int, nVars::Int, condSets::IntArrays)
    # get involved variables
    ind1, ind2 = triang2sub(nVars, ind)
    indFirst = min(ind1, ind2)
    indLast = max(ind1, ind2)

    # get conditioning set
    condSet = condSets[ind]

    if isempty(condSet)
        copName = "$indFirst,$indLast"
    else
        nCondVars = length(condSet)
        condSetSorted = sort(condSet)
        condSetString = "$(condSetSorted[1])"
        for ii=2:nCondVars
            condSetString = join([condSetString, ",", condSetSorted[ii]])
        end
        copName = "$indFirst,$indLast|$condSetString"
    end
    return copName
end

function getCopNam(ind::Int, pcc::PCC)
    condSets = Copulas.getAllCopCondSets(pcc.vine)
    nVars = dim(pcc)

    # get involved variables
    ind1, ind2 = triang2sub(nVars, ind)
    indFirst = min(ind1, ind2)
    indLast = max(ind1, ind2)

    # get conditioning set
    condSet = condSets[ind]

    if isempty(condSet)
        copName = "$indFirst,$indLast"
    else
        nCondVars = length(condSet)
        condSetSorted = sort(condSet)
        condSetString = "$(condSetSorted[1])"
        for ii=2:nCondVars
            condSetString = join([condSetString, ",", condSetSorted[ii]])
        end
        copName = "$indFirst,$indLast|$condSetString"
    end
    return copName
end

##############
## getindex ##
##############

import Base.getindex
@doc doc"""
Getindex functions called with i,j generally return the copula of i
and j. Therefore, subscript indices i and j need to be converted to
linear indexing of upper triangular matrices.
"""->

function getindex(pcc::PCC, ind1::Int64, ind2::Int64)
    # return pair copula for pair of variables
    n = dim(pcc)
    singleInd = []
    if ind1 < ind2
        singleInd = sub2triang(n, ind1, ind2)
    else
        singleInd = sub2triang(n, ind2, ind1)
    end
    return pcc.pccCops[singleInd]
end

import Base.setindex!
function setindex!(pcc::PCC, cop::PairCop,
                   ind1::Int64, ind2::Int64)
    # return pair copula for pair of variables
    n = dim(pcc)
    singleInd = []
    if ind1 < ind2
        singleInd = sub2triang(n, ind1, ind2)
    else
        singleInd = sub2triang(n, ind2, ind1)
    end
    pcc.pccCops[singleInd] = cop
end

function getParNot(pcc::PCC)
    return convert(CTreeParRef, pcc.vine)
end

## conditional probability transforms
##-----------------------------------

# evaluate single d-dimensional point
function pdf(pcc::PCC, vals::Array{Int, 1})
    # check dimensions
    n = dim(pcc)
    if length(vals) != n
        error("Point must coincide with pcc dimension.")
    end

    # iterate over all copulas with increasing layer
    condSets1, condSets2 = getAllVarCondSets(pcc.vine)
    pdfVals = NaN*ones(size(condSets1))

    maxCondSetLen = nVars - 2
    for setLen=0:maxCondSetLen # for each layer
        # get all conditioning sets of this layer
        inds = find([length(cs) == setLen for cs in condSets1])

        for copInd in inds
            # get associated variables
            var1, var2 = triang2sub(n, copInd)

            # get PITs
            if !isempty(condSets1[copInd])
                pit1 = getPit(pcc, var1, condSets1[copInd], vals)
                pit2 = getPit(pcc, var2, condSets2[copInd], vals)
            else
                pit1 = vals[var1]
                pit2 = vals[var2]
            end

            # calculate copula density
            cop = pcc.pccCops[copInd]
            pdfVals[copInd] = pdf(cop, pit1, pit2)
        end
    end
    return prod(pdfVals)
end

function getPit(pcc::PCC, condVar::Int, condSet::Array{Int, 1},
                vals::Array{Float64, 1})
    if isempty(condSet)
        # do nothing
        return vals[condVar]
    else
        # assure correct sequence of condSet
        trPar = CTreeParRef(pcc.vine.trees[:, condVar])
        sortedSet = findAndSortCondSet(trPar, condSet)

        # get last variable in conditioning set
        lastVar = sortedSet[end]

        # get PITs
        pit = []
        if length(sortedSet) > 1
            # get pits first
            pit1 = getPit(pcc, condVar, sortedSet[1:(end-1)], vals)
            pit2 = getPit(pcc, lastVar, sortedSet[1:(end-1)], vals)
        else
            pit1 = vals[condVar]
            pit2 = vals[lastVar]
        end
        if condVar < lastVar
            cop = pcc[condVar, lastVar]
            pit = hfun(cop, pit1, pit2)
        else
            cop = pcc[lastVar, condVar]
            pit = vfun(cop, pit1, pit2)
        end
    end
    return pit
end

## function getPit(pcc::PCC, condVar::Int, condSet::Array{Int, 1})
##     # get correct sequence of condSet
##     vnStruct = getParNot(pcc)

##     # get last variable in conditioning set
##     sortedSet = findAndSortCondSet(vnStruct[:, condVar], condSet)
##     lastVar = sortedSet[end]

##     ## create string for output
##     condSetString = "$(sortedSet[1])"
##     for ii=2:length(sortedSet)
##         condSetString = join([condSetString, ",", sortedSet[ii]])
##     end

##     println("$condVar|$condSetString")
##     if length(sortedSet) > 1
##         # get conditional PITs of lower layer
##         condVar_pit = getPit(pcc, condVar, sortedSet[1:end-1])
##         lastVar_pit = getPit(pcc, lastVar, sortedSet[1:end-1])
##     else
##         condVar_pit = condVar
##         lastVar_pit = lastVar
##     end

##     ## create string for output
##     condSetString = []
##     if length(sortedSet) == 1
##         condSetString = ""
##     else
##         condSetString = "$(sortedSet[1])"
##         for ii=2:(length(sortedSet)-1)
##             condSetString = join([condSetString, ",", sortedSet[ii]])
##         end
##     end

##     println("h($condVar,$lastVar|$condSetString)")
##     # get copula
##     ## vals = []
##     ## if condVar < lastVar
##     ##     cop = pcc.pairCops[condVar, lastVar]
##     ##     vals = hfun(cop, condVar_pit, lastVar_pit)
##     ## else
##     ##     cop = pcc.pairCops[lastVar, condVar]
##     ##     vals = vfun(cop, condVar_pit, lastVar_pit)
##     ## end
##     return condVar_pit
## end