Code.jl

GG = GAP.Globals
GG.LoadPackage(GAP.julia_to_gap("repsn"))
using Markdown

### An extra macro for the input checks 

macro req(args...)
  @assert length(args) == 2
  quote
    if !($(esc(args[1])))
      throw(ArgumentError($(esc(args[2]))))
    end
  end
end

##############################################################################
#                                                                            
# Helpers
#                                                                            
##############################################################################

### Extra converter

"""
  A conversion to transform elements of cyclotomic fields from GAP to Oscar
"""
function _convert_gap_to_julia(w_gap,F::AnticNumberField,z::nf_elem,n::Int64)
  L = GAP.gap_to_julia(GG.CoeffsCyc(w_gap,n))   
  w = sum([F(L[l]*z^(l-1)) for l=1:length(L) if L[l] !=0], init = F(0))
  return w
end

### Lists operations

"""
  Given a sorted list of integers, return the list of distinct integers in
  the last, and a second list giving the number of occurence of each value.
"""
function content(L::Vector{Int64})
  @assert issorted(L)
  Lu = unique(L)
  cont = Int64[]
  for i in Lu
    numb = count(j -> j == i, L)
    push!(cont, numb)
  end
  return Lu, cont
end

content(a::Int64) = content([a])

"""
  An internal to return at which index appears an element `f` in a basis `B`.
"""
function _index(B::Vector{T}, f::T) where T
  @req f in B "f not in B"
  return findfirst(b -> b == f, B)
end

"""
  Return the list of elements in a list `A` which satisfies a function `f`.
"""
function keep(f, A)
  lis = [i for i=1:length(A) if f(A[i])]
  return A[lis]
end

"""
  An internal that checks is all the elements of a list are distinct.
"""
function _distinct(L)
  L2 = [L[1]]
  for l in L
    !(l in L2) ? (push!(L2,l)) : continue 
  end
  return L2 == L
end

##############################################################################
#
# Combinatoric
#
##############################################################################

### Elevators

"""
  An internal recursive function for the elev method
"""

function _rec(A::Vector{Vector{Int64}},L::Vector{Int64},k::Int64)
  Ind_inter = typeof(A[1])[]
  bool = false
  for t in A
    L_inter = Int64[L[t[l]] for l=1:length(t)]
    if sum(L_inter) == k
      push!(Ind_inter,t)
    else
      for l = (t[end]):(length(L))
        if sum(L_inter) + L[l] <= k
          U = deepcopy(t)
          bool = true
          push!(U,l)
          push!(Ind_inter, U)
        end
      end
    end
  end
  bool == false ? (return Ind_inter) : (return _rec(Ind_inter, L, k))
end

"""
    elev(L::Vector{Int64}, k::Int64) -> Vector{Vector{Int64}}

Given a sorted list of integers `L` and an integer `k`, return all the set of
indices for which the corresponding elements in `L` sum up to `k`. Such sets of
indices are ordered increasingly (not strictly) and the output is ordered in
the lexicographic order from the left.
"""
function elev(L::Vector{Int64}, k::Int64)
  @req issorted(L) "L has to be a sorted list of integers"
  Ind = Vector{Int64}[[i] for i = 1:(length(L))]
  return _rec(Ind,L,k)
end

### For exterior powers

"""
  Use to construct all all tuples of `k` distinct among `1, ...,n` for exterior
  powers.
"""
function _rec_perm(n::Int64,k::Int64)
  if n==0 || k>n || k <= 0
    return Int64[]
  elseif k == 1
    return [[i] for i=1:n]
  elseif k == n
    return [[i for i=1:n]]
  else
    L = Vector{Int64}[]
    for j =k:n
      list1 = _rec_perm(j-1,k-1)
      for l in list1
        push!(l,j)
        push!(L,l)
      end
    end
    return L
  end
end

###############################################################
#
#  Tool to construct all possible character of a certain degree
#
###############################################################

### Get list information 

"""
  Given a sorted list of integers `L` and an integer `d`, return the number of
  distinct possible sums of elements in `L` to `d`, the type of this sums, and
  the cumulative numbers of type of sums
"""
function _number_sums(L::Vector{Int64}, d::Int64)
  Lu, numb = content(L)
  El = elev(Lu, d)
  sums = [Lu[el] for el in El]
  len = 0
  cumul = Int64[0]
  sumsum = []
  for sum in sums
    if length(sum) == 1
      for i in [i for i in 1:length(L) if L[i] == sum[1]]
        len += 1
        push!(cumul, len)
        push!(sumsum, (i, sum))
      end
    else
      it = 0
      LL = copy(L)
      while sum[1] in LL
        j = _index(LL, sum[1])
        LLu, numbLL = content(LL)
        Lu2, numb2 = content(sum[2:end])
        len += prod([binomial(numbLL[_index(LLu, Lu2[i])]+numb2[i]-1,numb2[i]) for i=1:length(Lu2)])
        LL = LL[j+1:end]
        it += j
        push!(cumul, len)
        push!(sumsum, (it, sum))
      end
    end
  end
  return len, cumul, sumsum
end

### "Elevations" for big lists

mutable struct ElevCtx
  L::Vector{Int64}
  d::Int64
  length::Int64
  cumul_length::Vector{Int64}
  sums::Vector{Tuple{Int64, Vector{Int64}}}

  function ElevCtx(L::Vector{Int64}, d::Int64)
    z = new(L, d)
    len, cumul, sumsum = _number_sums(L, d)
    z.length = len
    z.cumul_length = cumul
    z.sums = sumsum
    return z
  end
end

len(EC::ElevCtx) = EC.length
cumulated_lengths(EC::ElevCtx) = EC.cumul_length
base_list(EC::ElevCtx) = EC.L
degree(EC::ElevCtx) = EC.d
possible_sums(EC::ElevCtx) = EC.sums

function _rec_get_index(EC2::ElevCtx, i::Int64, sumtype::Vector{Int64})
  @assert 0 < i <= len(EC2)
  cumul2 = cumulated_lengths(EC2)
  sumsum2 = possible_sums(EC2)
  j = findfirst(j -> cumul2[j+1] >= i > cumul2[j] && sumsum2[j][2] == sumtype, 1:length(cumul2)-1)
  indexL2, sumtype2 = sumsum2[j]
  if length(sumtype2) == 1
    return [indexL2]
  end
  d = degree(EC2)
  L2 = base_list(EC2)
  el2 = [indexL2]
  EC3 = ElevCtx(L2[indexL2:end], d-L2[indexL2])
  i2 = -cumul2[j]
  i2 += sum([cumul2[k+1]-cumul2[k] for k =1:j-1 if sumsum2[k][1] == indexL2])
  suite = [indexL2 - 1 + i for i = _rec_get_index(EC3, i+i2, sumtype2[2:end])]
  el2 = vcat(el2, suite)
  return el2
end

function Base.getindex(EC::ElevCtx, i::Int64)
  @assert 0 < i <= len(EC) 
  cumul = cumulated_lengths(EC)
  j = findfirst(j -> cumul[j+1] >= i > cumul[j], 1:length(cumul)-1)
  sumsum = possible_sums(EC)
  indexL, sumtype = sumsum[j]
  if length(sumtype) == 1
    return [indexL]
  end
  d = degree(EC)
  L = base_list(EC)
  el = [indexL]
  EC2 = ElevCtx(L[indexL:end], d-L[indexL])
  i2 = -cumul[j]
  i2 += sum([cumul[k+1]-cumul[k] for k =1:j-1 if sumsum[k][1] == indexL])
  suite = [indexL - 1 + i for i = _rec_get_index(EC2, i+i2, sumtype[2:end])]
  el = vcat(el, suite)
  return el
end

Base.lastindex(EC::ElevCtx) = len(EC)

function Base.iterate(EC::ElevCtx, i::Int64 = 1)
  if i > len(EC)
    return nothing
  end
  return EC[i], i+1
end

##############################################################################
#
# Projective representations
#
##############################################################################


### Associated structures

@attributes mutable struct RepRing{S, T, U, V} 
  field::S
  group::T
  gens::U
  ct::V
  irr

  function RepRing(E::T) where {T <: Oscar.GapObj}
    e = GG.Exponent(E)
    F, _ = cyclotomic_field(e)
    ct = GG.CharacterTable(E)
    Irr = GG.Irr(ct)
    H = GG.GeneratorsOfGroup(E)
    RR = new{typeof(F), T, typeof(H), typeof(ct)}(F, E, H, ct, Irr)
    return RR
  end
end

splitting_field(RR::RepRing) = RR.field

underlying_group(RR::RepRing) = RR.group

generators_of_group(RR::RepRing) = RR.gens

character_table(RR::RepRing) = RR.ct

irreducible_characters(RR::RepRing) = RR.irr

polynomial_algebra(RR::RepRing) = get_attribute(RR, :poly_ring)

@attributes mutable struct LinRep{S, T, U, V, W} 
  rep_ring::RepRing{S, T, U, V}
  mats::Vector{MatElem}
  char_gap::W

  function LinRep(RR::RepRing{S,T,U, V}, MR::Vector{ <: MatElem}, CG::W) where {S, T, U, V, W}
    z = new{S, T, U, V, W}()
    z.rep_ring = RR
    z.mats = MR
    z.char_gap = CG
    return z
  end
end

representation_ring(LR::LinRep) = LR.rep_ring

matrix_representation(LR::LinRep) = LR.mats

char_gap(LR::LinRep) = LR.char_gap

dim(LR::LinRep) = LR.char_gap[1]

@attributes mutable struct ProjRep{S, T, U, V, W} 
  LR::LinRep{S, T, U, V, W}
  G::T

  function ProjRep(LR::LinRep{S, T, U, V, W}, G::T) where {S, T, U, V, W}
    z = new{S, T, U, V, W}()
    z.LR = LR
    z.G = G
    return z
  end
end

lift(PR::ProjRep) = PR.LR

underlying_group(PR::ProjRep) = PR.G

representation_ring(PR::ProjRep) = representation_ring(PR.LR)

matrix_representation(PR::ProjRep) = matrix_representation(PR.LR)

char_gap(PR::ProjRep) = char_gap(PR.LR)

dim(PR::ProjRep) = dim(PR.LR)

function Base.show(io::IO, ::MIME"text/plain", RR::RepRing)
  println(io, "Representation ring of")
  println(io, "$(RR.group)")
  println(io, "over")
  print(io, "$(RR.field)")
end

function Base.show(io::IO, RR::RepRing)
  print(io, "Representation ring of finite group over a field of characteristic 0")
end

function Base.show(io::IO, ::MIME"text/plain", LR::LinRep)
  println(io, "Linear $(dim(LR))-dimensional representation of")
  println(io, "$(LR.rep_ring.group)")
  println(io, "over")
  println(io, "$(LR.rep_ring.field)")
  println(io, "with matrix representation")
  print(io, "$(LR.mats)")
end 

function Base.show(io::IO, LR::LinRep)
  print(io, "Linear representation of finite group of dimension $(dim(LR))")
end

function Base.show(io::IO, ::MIME"text/plain", PR::ProjRep)
  println(io, "Linear lift of a $(dim(PR))-dimensional projective representation of")
  println(io, "$(PR.G)")
  println(io, "over")
  println(io, "$(PR.LR.rep_ring.field)")
  println(io, "with matrix representation")
  print(io, "$(PR.LR.mats)")
end

function Base.show(io::IO, PR::ProjRep)
    print(io, "Linear lift of a projective representation of finite group of dimension $(dim(PR.LR))")
end


### Identification of matrix representations

"""
  An internal function computing the quotient of a finite matrix group `G` by 
  its subgroup of complex multiple of the identity.
"""
function _quo_proj(G::GAP.GapObj)
  e = GG.Exponent(G)
  Elem = GG.Elements(GG.Centre(G))
  mult_id = GAP.julia_to_gap([], recursive=true)
  for elem in Elem
    if GG.IsDiagonalMat(elem) && GG.Size(GG.Eigenvalues(GG.CyclotomicField(e),elem)) == 1
      GG.Add(mult_id,elem)
    end
  end
  return GG.FactorGroup(G,GG.Group(mult_id))
end

"""
  An internal function computing the GAP matrix group associated to the matrix
  form of a representation
"""
function _rep_to_gap_group(rep::Vector{T}) where T <: MatElem
  rep_gap = GAP.julia_to_gap(rep, recursive = true)
  return GG.Group(rep_gap)
end

"""
  An internal function which returns, if it exists, the GAP id of the matrix group
  generated by `rep` modulo the multiple of the identity
"""
function _id_group_rep(rep::Vector{T}) where T <: MatElem
  G_gap = _rep_to_gap_group(rep)
  if !GG.IsFinite(G_gap)
    return "Possibly infinite group"
  end
  G_gap = _quo_proj(G_gap)
  if !GG.IdGroupsAvailable(GG.Size(G_gap))
    return "Group of order $(GG.Size(G_gap))"
  else
    return GG.IdGroup(G_gap)
  end
end

"""
  An internal functions checking whether the matrix form `rep` of a representation
  defines a projectively faitful representation of `GG.SmallGroup(o,n)`.
"""
function _is_good_pfr(rep::Vector{T}, (o,n)::Tuple{Integer, Integer}) where T <: MatElem
  G_gap = _rep_to_gap_group(rep)
  if !GG.IsFinite(G_gap)
    return false
  end
  G_gap = _quo_proj(G_gap)
  GG.Size(G_gap) != o ? (return false) : (return GAP.gap_to_julia(GG.IdGroup(G_gap)) == [o,n])
end

### Projectively faithful representations

"""
  An internal function which checks whether a small group `G_gap` admits
  faitful projective representations of dimension `dim`. If yes, it returns the
  necessary tools to create the projectively faithful representations of a Schur
  cover associated to those faithful projective representations
"""
function _has_pfr(G_gap::GAP.GapObj, dim::Int64; setup = nothing)
  f = GG.EpimorphismSchurCover(G_gap)
  H = GG.Source(f)
  n, p = ispower(GG.Size(H))
  if isprime(p)
    G_cover = GG.Image(GG.EpimorphismPGroup(H, p))
  else
    G_cover = GG.Image(GG.IsomorphismPermGroup(H))
  end
  RR = RepRing(G_cover)
  ct = character_table(RR)
  Irr = irreducible_characters(RR)
  L = [irr[1] for irr in Irr]
  ps = sortperm(L)
  L, Irr = L[ps], Irr[ps]
  RR.irr = Irr
  L = [i for i in L if i<= dim]
  lin_gap = [Irr[i] for i=1:length(L) if L[i] == 1]
  set_attribute!(RR, :lin_gap, lin_gap)
  set_attribute!(RR, :poly_ring, PolynomialRing(RR.field, "x" => 1:dim, cached = false)[1])
  keep_index = Int64[]
  sum_index = Vector{Int64}[]
  char_gap = GAP.GapObj[]
  if setup !== nothing
    @assert setup isa Tuple{Int64, Int64}
    d, t = setup
    cds = Vector{Tuple{Int64, typeof(Irr[1])}}[]
  end
  local El::ElevCtx
  El = ElevCtx(L, dim)
  @info "Classify $(len(El)) many projective representations."
  for l in El
    chi = sum(Irr[l])
    if any(lin -> lin*chi in char_gap, lin_gap)
      continue
    end
    Facto = GG.FactorGroup(G_cover, GG.CenterOfCharacter(chi))
    if !((GG.Size(Facto) == GG.Size(G_gap)) && (GG.IdGroup(Facto) == GG.IdGroup(G_gap)))
      continue
    end
    if setup != nothing
      chi_hom = GG.SymmetricParts(ct, GAP.julia_to_gap([chi], recursive=true), d)[1]
      cd = _character_decomposition(ct, chi_hom)
      filter!(pair -> pair[2][1] <= t, cd)
      cd in cds ? continue : push!(cds, cd)
    end    
    push!(char_gap, chi)
    keep_index = union(keep_index, l)
    push!(sum_index, l)
  end
  bool = length(keep_index) > 0
  return bool, RR, keep_index, sum_index
end

"""
    projectively_faithful_representations(o::Int64, n::Int64, dim::Int64)
                      -> Vector{Vector{MatElem{nf_elem}}, Vector{GAP.GapObj}, 
                         Vector{Vector{MatElem{nf_elem}}, Vector{GAP.GapObj}

Given a small group `G` of GAP index `[o,n]`, return representatives of all classes of
complex projectively faithful representations of a Schur cover `E` of `G` of dimension `dim`.
"""
function projectively_faithful_representations(o::Int64, n::Int64, dim::Int64; setup = nothing)
  G_gap = GG.SmallGroup(o,n)
  bool, RR, keep_index, sum_index = _has_pfr(G_gap, dim, setup = setup)
  if bool == false
    return ProjRep[]
  end
  @info "Construct and compare $(length(sum_index)) many projectively faithful representations."
  F = splitting_field(RR)
  e = get_attribute(F, :cyclo)
  z = gen(F)
  H_cover = generators_of_group(RR)
  ct = character_table(RR)
  Irr = irreducible_characters(RR)
  G_cover = underlying_group(RR)
  Rep = Vector{MatElem{nf_elem}}[]
  lin_gap = get_attribute(RR, :lin_gap)
  lin_oscar = Vector{MatElem{nf_elem}}[]
  for lin in lin_gap
    rep = GG.IrreducibleAffordingRepresentation(lin)
    Mat_cover = [GG.Image(rep,h) for h in H_cover]
    Mat = MatElem{nf_elem}[]
    for u =1:(length(Mat_cover))
      M_cover = Mat_cover[u]
      k = length(M_cover)
      M = [_convert_gap_to_julia(M_cover[i,j],F,z,Int(e)) for i=1:k, j =1:k]
      push!(Mat, matrix(F,M))
    end
    Mat = [m for m in Mat]
    push!(lin_oscar, Mat)
  end
  set_attribute!(RR, :lin_oscar, lin_oscar) 
  for index in keep_index
    irr = Irr[index]
    j = findfirst(j -> irr == lin_gap[j], 1:length(lin_gap))
    if j === nothing
      rep = GG.IrreducibleAffordingRepresentation(irr)
      Mat_cover = [GG.Image(rep,h) for h in H_cover]
      Mat = MatElem{nf_elem}[]
      for u =1:(length(Mat_cover))
        M_cover = Mat_cover[u]
        k = length(M_cover)
        M = [_convert_gap_to_julia(M_cover[i,j],F,z,Int(e)) for i=1:k, j =1:k]
        push!(Mat, matrix(F,M))
      end
      Mat = [m for m in Mat]
    else
      Mat = lin_oscar[j]
    end
    push!(Rep, Mat)
  end
  pfr = ProjRep[]
  for l in sum_index
    chara = sum(Irr[l])
    @assert chara[1] == dim
    Mat = Rep[_index( keep_index, l[end])]
    for i=length(l)-1:-1:1
      Mat = [block_diagonal_matrix([Mat[j], Rep[_index(keep_index, l[i])][j]]) for j=1:length(Mat)]
    end
    #MG = matrix_group(Mat)
    keep = true
    for pr in pfr
      Mat_inter = matrix_representation(pr)
      G_inter = matrix_group(Mat_inter)
      all(mm -> mm in G_inter, Mat) ? keep = false : nothing
    end
    if keep == true
      lr = LinRep(RR, Mat, chara)
      pr = ProjRep(lr, G_gap)
      push!(pfr, pr)
    end
  end
  return pfr
end

### Projective representations 

"""
  An internal function which checks whether a small group `G_gap` admits projective
  representations of dimension `dim` with kernel of dimension at most `index_max`.
  If yes, it returns the tools necessary to compute the linear representations
  of a Schur cover corresponding to those projective representations
"""
function _has_pr(G_gap::GAP.GapObj, dim::Int64; index_max::IntExt = inf)
  f = GG.EpimorphismSchurCover(G_gap)
  H = GG.Source(f)
  n, p = ispower(GG.Size(H))
  if isprime(p)
    G_cover = GG.Image(GG.EpimorphismPGroup(H, p))
  else
    G_cover = GG.Image(GG.IsomorphismPermGroup(H))
  end
  RR = RepRing(G_cover)
  ct = character_table(RR)
  Irr = irreducible_characters(RR)
  Irr = keep(irr -> irr[1] <= dim, Irr)
  L = [irr[1] for irr in Irr]
  ps = sortperm(L)
  L, Irr = L[ps], Irr[ps]
  RR.irr = Irr
  lin_gap = [Irr[i] for i=1:length(Irr) if Irr[i][1] == 1]
  set_attribute!(RR, :lin_gap, lin_gap)
  set_attribute!(RR, :poly_ring, PolynomialRing(RR.field, "x" => 1:dim, cached = false)[1])
  keep_index = Int64[]
  sum_index = Vector{Int64}[]
  char_gap = GAP.GapObj[]
  local El::ElevCtx
  El = ElevCtx(L, dim)
  @info "Classify $(len(El)) many projective representations."
  for l in El
    chi = GG.Character(G_cover, sum(Irr[l]))
    @assert chi[1] == dim
    keep = true
    Facto = GG.FactorGroup(G_cover, GG.CenterOfCharacter(chi))
    idx, r = divrem(GG.Size(G_gap), GG.Size(Facto))
    if (r == 0 && idx <= index_max)
      if idx == 1
        if GG.IdGroup(Facto) != GG.IdGroup(G_gap)
          keep = false
        end
      end
      if any(lin -> lin*chi in char_gap, lin_gap)
        keep = false
      end
    else
      keep = false
    end
    if keep == true
      push!(char_gap, chi)
      keep_index = union(keep_index, l)
      push!(sum_index, l)
    end
  end
  bool = length(keep_index) > 0
  return bool, RR, keep_index, sum_index
end

"""
  Return projective representations of a group, given as linear representations
  of a Schur cover, for which the kernel as maximal index `index_max` in the 
  group
"""
function projective_representations(o::Int64, n::Int64, dim::Int64; index_max::IntExt = inf)
  @req (o >= index_max || index_max == inf) "index_max must be smaller than the order of the group"
  index_max = min(o, index_max)
  if (index_max == 1)
    return projectively_faithful_representations(o, n, dim)
  end
  G_gap = GG.SmallGroup(o,n)
  bool, RR, keep_index, sum_index = _has_pr(G_gap, dim, index_max = index_max)
  if bool == false
    return ProjRep[]
  end
  @info "Construct $(length(sum_index)) many projective representations."
  F = splitting_field(RR)
  e = get_attribute(F, :cyclo)
  z = gen(F)
  H_cover = generators_of_group(RR)
  ct = character_table(RR)
  Irr = irreducible_characters(RR)
  G_cover = underlying_group(RR)
  Rep = Vector{MatElem{nf_elem}}[]
  lin_gap = get_attribute(RR, :lin_gap)
  lin_oscar = Vector{MatElem{nf_elem}}[]
  for lin in lin_gap
    rep = GG.IrreducibleAffordingRepresentation(lin)
    Mat_cover = [GG.Image(rep,h) for h in H_cover]
    Mat = MatElem{nf_elem}[]
    for u =1:(length(Mat_cover))
      M_cover = Mat_cover[u]
      k = length(M_cover)
      M = [_convert_gap_to_julia(M_cover[i,j],F,z,Int(e)) for i=1:k, j =1:k]
      push!(Mat, matrix(F,M))
    end
    Mat = [m for m in Mat]
    push!(lin_oscar, Mat)
  end
  set_attribute!(RR, :lin_oscar, lin_oscar)
  for index in keep_index
    irr = Irr[index]
    j = findfirst(j -> irr == lin_gap[j], 1:length(lin_gap))
    if j === nothing
      rep = GG.IrreducibleAffordingRepresentation(irr)
      Mat_cover = [GG.Image(rep,h) for h in H_cover]
      Mat = MatElem{nf_elem}[]
      for u =1:(length(Mat_cover))
        M_cover = Mat_cover[u]
        k = length(M_cover)
        M = [_convert_gap_to_julia(M_cover[i,j],F,z,Int(e)) for i=1:k, j =1:k]
        push!(Mat, matrix(F,M))
      end
      Mat = [m for m in Mat]
    else
      Mat = lin_oscar[j]
    end
    push!(Rep, Mat)
  end
  pfr = ProjRep[]
  for l in sum_index
    chara = sum(Irr[l])
    @assert chara[1] == dim
    Mat = Rep[_index( keep_index, l[end])]
    for i=length(l)-1:-1:1
      Mat = [block_diagonal_matrix([Mat[j], Rep[_index(keep_index, l[i])][j]]) for j=1:length(Mat)]
    end
    bool = true
    for pr in pfr
      Mat_inter = matrix_representation(pr)
      G_inter = matrix_group(Mat_inter)
      it2 = count(mat -> mat in G_inter, Mat)       
      it2 == length(Mat) ? bool = false : nothing
    end
    if bool == true
      lr = LinRep(RR, Mat, chara)
      pr = ProjRep(lr, G_gap)
      push!(pfr, pr)
    end
  end
  return pfr
end

##############################################################################
#
# Polynomial algebra tools
#
##############################################################################

### Basis for homogeneous components

"""
    homog_basis(R::MPolyRing{nf_elem}, d::Int64}) -> Vector{MPolyElem{nf_elem}}

Given a multivariate polynomial ring `R` with the standard grading, return a 
basis of the `d`-homogeneous part of `R`, i.e. all the monomials of homogeneous
degree `d`. The output is ordered in the lexicographic order from the left for the indices
of the variables of `R` involved.
"""
function homog_basis(R::T,d::Int64) where T <: MPolyRing
  n = nvars(R)
  L = Int64[1 for i =1:n]
  E = elev(L,d)
  G = gens(R)
  B = typeof(G[1])[prod(G[e]) for e in E]
  return B
end

### Functions for conversion between abstract homogeneous polynomials and their coordinates
### to the chosen bases

"""
  An internal to transform a set of polynomials whose coordinates in a basis
  of the corresponding polynomial algebra `R` are the columns of the given matrix `M`.
"""
function _mat_to_poly(M::MatElem{T}, R::MPolyRing{T}) where T <: RingElem
  @assert nrows(M) == nvars(R) 
  G = gens(R)
  return typeof(G[1])[(change_base_ring(R,transpose(M)[i,:])*matrix(G))[1] for i=1:(size(M)[2])]
end

"""
  An internal to change a homogeneous polynomials into a vector whose entries are the 
  coefficients of `P` in the canonical basis of the corresponding homogeneous component
"""
function _homo_poly_in_coord(P::T) where T 
  R = parent(P)
  F = base_ring(R)
  d = total_degree(P)
  B = homog_basis(R,d)
  v = zeros(F,length(B), 1)
  Coll = collect(terms(P))
  for c in Coll
    v[_index(B, collect(monomials(c))[1])] = collect(coeffs(c))[1]
  end
  return matrix(v)
end

"""
  An internal to change an homogeneous polynomial of `R` of degree `d` given in 
  coordinates `w` into its polynomial form.
"""
function _coord_in_homo_poly(w::MatElem{T}, R::MPolyRing, d::Int64) where T <: RingElem
  B = homog_basis(R,d)
  @assert length(B) == size(w)[1]
  return sum([w[:,1][i]*B[i] for i=1:(size(w)[1])])
end

"""
  An internal to transform an homogeneous ideal into a matrix whose columns are the coordinates
  of homogenoues generators of the same degree of `I` in an appropriate basis
"""
function _homo_ideal_in_coord(I::T) where T <:MPolyIdeal
  g = gens(I)
  R = base_ring(I)
  F = base_ring(R)
  d = total_degree(g[1])
  @req all(gg -> sum(degrees(collect(terms(gg))[1])) == d, g) "The generators of `I` must be of the same homogeneous degree"
  B = homog_basis(R,d)
  v = zeros(F,length(B),length(g))
  for i=1:(length(g))
    Coll = collect(terms(g[i]))
    vf = zeros(F, length(B), 1)
    for c in Coll
      vf[_index(B, collect(monomials(c))[1])] = collect(coeffs(c))[1]
    end
    v[:,i] = vf
  end
  return (v,B)
end

##############################################################################
#
# About tensors of homogeneous polynomials
#
##############################################################################

### Basis of exterior power

"""
   basis`_`exterior`_`power(B::Vector{T}, k::Int64) where T -> Vector{Vector{T}}

Given a basis `B` of a vector space `V`, return a basis of the `t`-exterior power of
`V` as a set of lists. Each list is ordered in strictly increasing order of the 
indices of the corresponding elements in `B`. The output is ordered in the lexicographic
order from the right in the indices of the corresponding elements in `B`
"""
function basis_exterior_power(B::Vector{T}, t::Int64) where T
  l = length(B)
  L = _rec_perm(l,t)
  return Vector{T}[B[lis] for lis in L]
end

### Operations on pure tensors

"""
  An internal to check if two pure tensors define the same element up to
  reordering.
"""
function _same_base(v::Vector{T}, w::Vector{T}) where T
  @assert length(v) == length(w)
  l = count(vv -> vv in w, v)
  return l == length(v)
end

"""
  An internal that, given two elements representing the same pure tensors, return
  the sign of the permutation of components from one to another.
"""
function _div(v::Vector{T}, w::Vector{T}) where T
  @assert _same_base(v,w)
  return sign(perm([findfirst(vv -> vv == ww, v) for ww in w]))
end

"""
  An internal to check whether or not a pure tensor is a basis tensor. If yes, it
  returns the given basis tensor and the sign of the permutation from one to another.
"""
function _in_basis(v::Vector{T}, basis::Vector{Vector{T}}) where T
  @assert length(v) == length(basis[1])
  i = findfirst(w -> _same_base(v,w),basis)
  i == nothing ? (return false) : (return basis[i], _div(v,basis[i]))
end

"""
  An internal that given an element `w` of the `t`-exterior power of a vector space 
  expressed by its coordinates in a basis `basis` of this exterior power, returns it in the
  form of a set of pairs consisting of the coefficients and the corresponding element of 
  the basis (given as a matrix as in `_base_to_columns`) read from the coordinates `w`.
"""
function _change_basis(w, basis::Vector{Vector{T}}) where T 
  @assert length(w) == length(basis) 
  gene = [(w[i],basis[i]) for i =1:length(w) if w[i] != parent(w[i])(0)]
  return gene
end

function _standard_basis(F, n)
  v = zero_matrix(F,n,1)
  std_bas = typeof(v)[]
  for j=1:n
    vv =deepcopy(v)
    vv[j,1] = one(F)
    push!(std_bas, vv)
  end
  return std_bas
end

function _decompose_in_standard_basis(v)
  std_bas = _standard_basis(base_ring(v), nrows(v))
  return std_bas, [[v[i], std_bas[i]] for i=1:nrows(v)]
end

function _wedge_product(V::Vector)
  if length(V) == 1
    v = V[1]
		_, dec = _decompose_in_standard_basis(v)
    w = [(l[1], [l[2]]) for l in dec]
    return w
  end
  v = popfirst!(V)
  wp = _wedge_product(copy(V))
  std_bas, dec = _decompose_in_standard_basis(v)
  bas = basis_exterior_power(std_bas, length(V)+1)
  coeffs = zeros(parent(dec[1][1]), length(bas)) 
  for l in dec
    for w in wp
      if l[2] in w[2]
        continue
      end
      w_co = copy(w[2])
      lw_co = pushfirst!(w_co, l[2])
      ba, mult = _in_basis(lw_co, bas)
      coeffs[_index(bas, ba)] += mult*w[1]*l[1]
    end
  end
  return _change_basis(coeffs, bas)
end

function _shuffle_equation(R::MPolyRing, L::PolyRing{nf_elem}, M::MatElem, t::Int64)
  F = base_ring(M)
  @assert base_ring(L) === F
  y = gen(L)
  T = identity_matrix(F, nrows(M))
  Z = map_entries(L, T) + y*map_entries(L, M)
  Q = zero_matrix(L, binomial(nrows(M), t), binomial(ncols(M), t))
  std_bas = _standard_basis(F, nrows(M))
  bas = basis_exterior_power(std_bas, t)
  r = _rec_perm(nrows(M), t)
  for i in 1:length(r)
    l = r[i]
    V = [Z[:,j] for j in l]
    wp = _wedge_product(V)
    for w in wp
      k = _index(bas, [map_entries(p -> evaluate(p,0), ww) for ww in w[2]])
      Q[k,i] = w[1]
    end
  end

  Q -= identity_matrix(L, nrows(Q))
  I = ideal(R, [R(0)])
  x = gens(R)
  
  for j=1:t
    gene = typeof(x[1])[]
    for k=1:ncols(Q)
      my = Q[:,k]
      m = map_entries(p -> divexact(p, y), my)
      v = map_entries(p -> evaluate(p, 0), m)
      Q[:,k] -= y*map_entries(L, v)
      if !iszero(v)
        push!(gene, sum([v[l]*x[l] for l=1:length(x)]))
      end
    end
    J = radical(ideal(R,gene))
    I = radical(I+J)
  end
  return I
end

### wedge product of polynomials

"""
   wedge_product(F::Vector{T}) where T <: MPolyElem 
                                 -> Vector{Tuple{nf_elem}, Vector{MPolyElem}}
                                 
Given a list of homogeneous polynomials `F` of the same polynomial algebra and of the same
total degree, return their wedge product as an abstract tensor, seen as a list of tuples
of coefficients and the respective elements of the basis of the exterior power.
"""
function wedge_product(F::Vector{T}) where T <: MPolyElem
  if length(F) == 1
    f = F[1]
    R = parent(f)
    d = total_degree(f)
    B = homog_basis(R,d)
    basis = basis_exterior_power(B, 1)
    list_mon, list_coeff = collect(monomials(f)), collect(Oscar.coefficients(f))
    part= [(list_coeff[i], [list_mon[i]]) for i=1:length(list_mon)]
    return part
  end
  f = popfirst!(F)
  wp = wedge_product(copy(F))
  R = parent(f)
  d = total_degree(f)
  B = homog_basis(R,d)
  basis = basis_exterior_power(B, length(F)+1)
  list_mon, list_coeff = collect(monomials(f)), collect(Oscar.coefficients(f))
  part = [(list_coeff[i], list_mon[i]) for i=1:length(list_mon)]
  coeffs2 = zeros(base_ring(R), length(basis))
  for l in part
    for w in wp
      if l[2] in w[2]
        continue
      end
      w_co = copy(w[2])
      lw_co = pushfirst!(w_co, l[2])
      ba, mult = _in_basis(lw_co, basis)
      coeffs2[_index(basis, ba)] += mult*w[1]*l[1]
    end
  end
  return _change_basis(coeffs2, basis)
end  

### Operations on arbitrary tensors

"""
  An internal that computes the sum of two tensors in the same tensor space.
"""
function _sum_tensors(v::Vector{T}, w::Vector{T}) where T
  wbis = copy(w)
  for vv in v
    i = findfirst(ww -> vv[2] == ww[2],w)
    if i == nothing
      push!(wbis, vv)
    else
      wbis[i][1] += vv[1]
    end
  end
  return wbis
end

"""
  An internal that changes the coefficient ring of the coefficients of a tensor, if 
  this change is possible.
"""
function _change_coeff_ring(w::Vector{T}, L) where T
  return [[L(ww[1]), ww[2]] for ww in w]
end

"""
  An internal to rescale the coefficients of a tensor by a same element.
"""
function _rescale(w::Vector{T}, r) where T
  @assert parent(w[1][1]) == parent(r)
  return [[r*ww[1],ww[2]] for ww in w]
end

"""
   is_pure(w::Vector{Vector{SetElem}}, basis::Vector{Vector{T}}, B::Vector{T})
           where T <: MPolyElem
                        -> Tuple{Bool, Vector{T}}

Given an element `w` of an exterior power of an homogeneous part of a polynomial algebra
given as a set of pairs consisting of a coefficient and an element of `basis` in matrix form, 
return whether it defines a pure tensor or not.

If it is, it returns `(true,dec)` where `dec` is a vector of polynomials whose wedge product
is the completely decomposed form of `w`.
If not, it returns `(false, w)`
"""
function is_pure(w, basis::Vector{Vector{T}}) where T <: MPolyElem
  R = parent(basis[1][1])
  d = total_degree(basis[1][1])
  B = homog_basis(R,d)
  basis2 = basis_exterior_power(B, length(basis[1])+1)
  m, n = length(B), length(basis[1])
  F  = parent(w[1][1])
  Mat = matrix(zeros(F,length(basis2),length(B)))
  for i=1:length(B)
    col = zeros(F,length(basis2),1)
    v = B[i]
    for ww in w
      if v in ww[2]
        continue
      end
      vv, mult = _in_basis(vcat(v, ww[2]), basis2)
      col[_index(basis2,vv)] = ww[1]*mult
      Mat[:,i] = col
    end
  end
  if Oscar.rank(Mat) != m-n
    return (false, w)
  else
    K = kernel(Mat)[2]
    kern = [K[:,j] for j=1:(size(K)[2])]
    return (true, [_coord_in_homo_poly(k,R,d) for k in kern])
  end
end

### Manipulation of tensors by their coordinates to the basis of the exterior power

"""
  An internal to transform a vector representing a basis element `base` of an exterior power
  of a vector space into a matrix whose columns are the coordinates of each elements
  of the vector in a basis `B` of the corresponding vector space.
"""
function _base_to_columns(base::Vector{T}, B::Vector{T}) where T 
  columns = zeros(base_ring(base[1]),length(B), length(base))
  for i=1:length(base)
    columns[_index(B,base[i]),i] = base_ring(base[1])(1)
  end
  return matrix(columns)
end

### Helper for intersection with Grassmannian

"""
  An internal to compute a matrix of an isomorphism from the `t`-exterior power of a 
  `n`-dimensional homogeneous component of a polynomial algebra to the `(n-t)`-exterior 
  power of its dual.
"""
function _matrix_ext_powers(B::Vector{T},k::Int64, l::Int64) where T <: MPolyElem{nf_elem}
  @req length(B) == k+l "l should be equal to length(B)-k"
  basis = basis_exterior_power(B,k)
  basis2 = basis_exterior_power(B,l)
  F = base_ring(B[1])
  Mat = matrix(zeros(F, length(basis2), length(basis)))
  for j=1:length(basis)
    bt = basis[j]
    i = findfirst(i -> isempty(intersect(bt, basis2[i])), 1:length(basis2))
    b = vcat(bt,basis2[i])
    perma = Int64[_index(B,poly) for poly in b]
    PP = Perm(perma)
    Mat[i,j] = sign(PP)
  end
  return Mat
end

"""
  An internal to compute the matrix of the linear map from which we read the pure tensors 
  from the degeneracy locus at which the rank doesn't exceed a certain bound
"""
function _matrix_multivec_star(V, B, k::Int64)
  @assert base_ring(B[1]) == parent(V[1][1])
  R = parent(B[1])
  F = base_ring(R)
  l = length(B)-k
  Basis = basis_exterior_power(B,k)
  Basis2 = basis_exterior_power(B,l)
  Basis3 = basis_exterior_power(B,l+1)
  L,y = PolynomialRing(F,"y"=>(1:length(V)))
  w = sum([y[i]*map_entries(L, V[i]) for i=1:length(V)])
  Mat = _matrix_ext_powers(B,k,l)
  wstar = _change_basis(Mat*transpose(w), Basis2)
  Mat = matrix(zeros(L,length(Basis3),length(B)))
  for i=1:length(B)
    col = zeros(L,length(Basis3),1)
    v = B[i]
    for ww in wstar
      if v in ww[2]
        continue
      end
      vv, mult = _in_basis(vcat(ww[2],v), Basis3)
      col[_index(Basis3,vv)] = ww[1]*mult
      Mat[:,i] = col
    end
  end
  return Mat
end

##############################################################################
#
# Induced representations
#
##############################################################################

### Helper to compute induced representations on different vector spaces

function _induced_representation(prep::ProjRep, M::Vector{T}, chara::GAP.GapObj) where T <: MatElem
  lr = LinRep(representation_ring(prep), M, chara)
  return ProjRep(lr, underlying_group(prep))
end

### Induced representations on homogeneous component of polynomial algebra

"""
  An internal for computing dual representations in matrix form from those in `Rep`.
"""
function _dual_representation(prep::ProjRep)
  Rep = matrix_representation(prep)
  char = char_gap(prep)
  Rep_dual = [transpose(inv(rep)) for rep in Rep]
  char_dual = GG.ComplexConjugate(GAP.julia_to_gap([char], recursive = true))[1]
  return _induced_representation(prep, Rep_dual, char_dual)
end

"""
  An internal that compute the action of a group on an `d`-homogeneous part of
  a polynomial algebra `R` from a matrix form `M` of the action on the variables of `R`.
"""
function _action_poly_ring(M::MatElem{T},R::MPolyRing{T},d::Int64) where T <: RingElem
  @assert nrows(M) == nvars(R)
  B = homog_basis(R,d)
  v = _mat_to_poly(M,R)
  MM = zeros(base_ring(M), length(B), length(B))
  for i =1:(length(B))
    W = evaluate(B[i],v)
    C = collect(terms(W))
    for c in C
      MM[_index(B,collect(monomials(c))[1]),i] = collect(coeffs(c))[1]
    end
  end
  return matrix(MM)
end

_homo_poly_action(m, R, d) = _action_poly_ring(transpose(inv(m)), R, d)
"""
   homo`_`poly`_`rep(Rep:Vector{Vector{T}}, d::Int64, (char`_`gap, ct)::Tuple{Vector{S}, S})
                 where T <: MatElem, S <: GAP.GapObj
                    -> Vector{Vector{T}}, Vector{S}

Given a set `Rep` of linear representations of a group `G` on a finite dimensional complex 
vector space `V`, return the induced representations on the `d` homogeneous part of the dual
of `V` seen as a polynomial algebra.

- The first output is the set of induced representations given in Oscar matrix form;
- The second output is the set of the corresponding characters given as GAP objects
  computed from the input `char_gap`, the characters of the original representations,
  and `ct` the character table of `G`.
"""
function homo_poly_rep(prep::ProjRep, d::Int64)
  prep_dual = _dual_representation(prep)
  ct = character_table(representation_ring(prep))
  R = polynomial_algebra(representation_ring(prep))
  if R === nothing
    R, _ = PolynomialRing(splitting_field(representation_ring(prep)), "x" => 1:dim(prep), cached = false)
  end
  Rep_homo = [_action_poly_ring(m,R,d) for m in matrix_representation(prep_dual)]
  char_homo = GG.SymmetricParts(ct, GAP.julia_to_gap([char_gap(prep_dual)], recursive = true),d)[1]
  prep_homo = _induced_representation(prep_dual, Rep_homo, char_homo)
  set_attribute!(prep, :prep_homo, prep_homo)
  set_attribute!(prep_homo, :degree_homo, d)
  return prep_homo
end

### Induced representation on exterior power of homogeneous component of polynomials
### algebra

"""
    ext__homo_poly_rep(Rep::Vector{T}, R::MPolyring, d::Int64, t::Int64) 
                                               where T <: MatElem{nf_elem}
                                                                -> Vector{T}

Given a representation of a group on a the `d` homogeneous part of `R`,
return the induced representation on the `t`-exterior power of `R_d` given in
matrix form.
"""
function ext_homo_poly_rep(prep_homo::ProjRep, t::Int64)
  @req 1 <= t <= dim(prep_homo) "t must be contain between 1 and the dimension $(dim(prep_homo)) of prep"
  d = get_attribute(prep_homo, :degree_homo)
  @assert d !== nothing
  R = polynomial_algebra(representation_ring(prep_homo))
  B = homog_basis(R,d)
  basis = basis_exterior_power(B,t)
  chara = char_gap(prep_homo)
  ct = character_table(representation_ring(prep_homo))
  Rep_homo = matrix_representation(prep_homo)
  prep_ext = typeof(Rep_homo[1])[]
  for M in Rep_homo
    M_ext = zeros(base_ring(M), length(basis), length(basis))
    for i =1:length(basis)
      b = basis[i]
      bb = [_coord_in_homo_poly(M*_homo_poly_in_coord(bbb), R, d) for bbb in b]
      wp = wedge_product(bb)
      for bbb in wp
        M_ext[_index(basis, bbb[2]), i] += bbb[1]
      end
    end
    push!(prep_ext, matrix(M_ext))
  end
  chara_ext = GG.AntiSymmetricParts(ct, GAP.julia_to_gap([chara], recursive=true), t)[1]
  return _induced_representation(prep_homo, prep_ext, chara_ext)
end

function _reconstruct_character(G::Vector, RR::RepRing)
  MG = matrix_group(G)
  gene = [MG(m).X for m in G]
  f = GG.GroupHomomorphismByImages(underlying_group(RR), MG.X, generators_of_group(RR), GAP.julia_to_gap(gene))
  coll = []
  H = GG.ConjugacyClasses(underlying_group(RR))
  for h in H
    b = f(h)
    push!(coll, GG.TraceMat(GG.Matrix(b[1])))
  end
  c = GG.Character(character_table(RR), GAP.julia_to_gap(coll))
  return c
end

function _projective_representation(G::Vector, RR::RepRing; gap_grp = nothing)
  chara = _reconstruct_character(G, RR)
  lr = LinRep(RR, G, chara)
  if gap_grp === nothing
    gap_grp = _quo_proj(matrix_group(G).X)
  end
  return ProjRep(lr, gap_grp)
end

function _get_gap_rep(prep::ProjRep)
  RR = representation_ring(prep)
  gene = GAP.julia_to_gap([GAP.julia_to_gap(m, recursive=true) for m in matrix_representation(prep)])
  G2 = GG.Group(gene)
  f = GG.GroupHomomorphismByImages(underlying_group(RR), G2, generators_of_group(RR), gene)
  return GG.MappingByFunction(underlying_group(RR), G2, f)
end

function _small_sums_invertible(B::Vector)
  k = length(B)
  part = typeof(B[1])[]
  L = Int64[]
  max_comp = Int64[]
  for m in B
    push!(L, 1)
    push!(max_comp, 1)
  end
  for t in 1:k
    el = ElevCtx(L, t)
    for l in el
      countL = [count(i -> i == j, l) for j=1:length(L)]
      if any(i -> countL[i] > max_comp[i], 1:length(L))
        continue
      end
      m = sum([B[i] for i in l])
      if det(m) != 0
        return m
      end
    end
  end
  return nothing
end

function _find_base_change(prep::ProjRep, prep2::ProjRep)
  @assert representation_ring(prep) === representation_ring(prep2)
  @assert char_gap(prep) == char_gap(prep2)
  As = matrix_representation(prep)
  Bs = matrix_representation(prep2)
  linind = transpose(LinearIndices(size(As[1])))
  n = nrows(As[1])
  zz = zero_matrix(base_ring(As[1]), 0, n^2)
  for (A, B) in zip(As, Bs)
    z = zero_matrix(base_ring(A), n^2, n^2)
    for i in 1:n
      for j in 1:n
        for k in 1:n
          z[linind[i, j], linind[i, k]] += A[k, j]
          z[linind[i, j], linind[k, j]] -= B[i, k]
        end
      end
    end
    zz = vcat(zz, z)
  end
  zzs = sparse_matrix(zz)
  if sparsity(zzs) >= Float64(7//10)
    r, K = right_kernel(zzs)
  else
    r, K = right_kernel(zz)
  end
  res = eltype(As)[]
  for k in 1:ncols(K)
    cartind = Hecke.cartesian_product_iterator([1:x for x in size(As[1])],
                                         inplace = true)
    M = zero_matrix(base_ring(As[1]), nrows(As[1]), ncols(As[1]))
    for (l, v) in enumerate(cartind)
      M[v[2], v[1]] = K[l, k]
    end
    push!(res, M)
  end
  
  if length(res) == 1
    @assert det(res[1]) != 0
    return res[1]
  end

  j = findfirst(m -> det(m) != 0, res)
  if j !== nothing
    return res[j]
  end
  
  m = _small_sums_invertible(res)
  if m !== nothing
    return m
  end

  t = length(res)
  B = sum(res)
  while det(B) == 0
    v = rand(-t:t, t)
    B = sum([res[i]*v[i] for i in 1:t])
  end
  GC.gc()
  GC.gc()
  return B
end

function _to_equivalent_block_representation(prep::ProjRep)
  rho = _get_gap_rep(prep)
  rho2 = GG.EquivalentBlockRepresentation(rho)
  @info "Found equivalent block representation"
  RR = representation_ring(prep)
  F = splitting_field(RR)
  z = gen(F)
  e = get_attribute(F, :cyclo)
  H = generators_of_group(RR)
  coll = typeof(matrix_representation(prep)[1])[]
  m = [GG.Image(rho2, h) for h in H]
  for M in m
    MM = [_convert_gap_to_julia(M[i,j],F,z,Int(e)) for i=1:length(M), j=1:length(M)]
    push!(coll, matrix(F, MM))
  end
  prep2 = _projective_representation(coll, RR, gap_grp = underlying_group(prep))
  
  B = _find_base_change(prep, prep2)
  @info "Base change computed!"
  
  @assert all(i -> B*matrix_representation(prep)[i]*inv(B) == matrix_representation(prep2)[i], 1:length(H))
  GC.gc()
  GC.gc()
  return prep2, B
end

function _action_on_simple_submodule(prep::ProjRep, chi)
  @req GG.IsIrreducibleCharacter(chi) "chi is not irreducible"
  RR = representation_ring(prep)
  ct = character_table(RR)
  nu = char_gap(prep)
  m = matrix_representation(prep)
  @req GG.ScalarProduct(ct, nu, chi) != 0 "chi is not afforded by a simple submodule"
  
  ic = [irr for irr in irreducible_characters(RR) if irr[1] <= chi[1]]
  index = findfirst(mu -> mu == chi, ic)
  @assert index !== nothing
  
  vals = [GG.ScalarProduct(ct, irr, nu) for irr in ic]
  id_block = sum([vals[i]*ic[i][1] for i=1:index-1])
  mchi = [M[id_block+1:id_block+ic[index][1], id_block+1:id_block+ic[index][1]] for M in m]
  
  prepchi = _projective_representation(mchi, RR) 
  @assert char_gap(prepchi) == chi
  
  return prepchi
end
  
function _basis_simple_submodule(prep::ProjRep, K, chi)
  @req GG.IsIrreducibleCharacter(chi) "chi is not irreducible"
  RR = representation_ring(prep)
  F = splitting_field(RR)
  ct = character_table(RR)
  nu = char_gap(prep) 
  
  cd = _character_decomposition(prep)
  ic = [c[2] for c in cd]
  j = findfirst(mu -> mu == chi, ic)
  if j === nothing
    return zero_matrix(F, dim(prep), 0)
  end
  alpha = cd[j][1]
  
  vals = [c[1] for c in cd]
  id_block = sum([vals[i]*ic[i][1] for i=1:j-1], init = 0)
  M = zero_matrix(F, dim(prep), alpha*chi[1])
  M[id_block+1:id_block+alpha*chi[1], :] = identity_matrix(F, alpha*chi[1])
  
  return inv(K)*M
end

function _character_decomposition(prep::ProjRep)
  char = char_gap(prep)
  RR = representation_ring(prep)
  ct = character_table(RR)
  return _character_decomposition(ct, char)
end

function _character_decomposition(ct, char)
  lis = GG.ConstituentsOfCharacter(ct, char)
  decomp = Tuple{Int64, typeof(lis[1])}[]
  for chi in lis
    alpha = GG.ScalarProduct(ct, chi, char)
    push!(decomp, (alpha, chi))
  end
  return decomp
end

function _partition(prep::ProjRep, t::Int64)
  RR = representation_ring(prep)
  ct = character_table(RR)
  chi = char_gap(prep)
  return _partition(ct, chi, t)
end

function _partition(ct, chi, t)
  cd = _character_decomposition(ct, chi)
  part = typeof(cd[1][2])[]
  L = Int64[]
  max_comp = Int64[]
  for c in cd
    push!(L, c[2][1])
    push!(max_comp, c[1])
  end
  el = ElevCtx(L, t)
  for l in el
    countL = [count(i -> i == j, l) for j=1:length(L)]
    if any(i -> countL[i] > max_comp[i], 1:length(L))
      continue
    end
    chi = sum([cd[i][2] for i in l])
    push!(part, chi)
  end
  return part
end

############################################################################
#
#  Moduli space
#
############################################################################

struct ParameterMapIsotypical
  f::Function
  U::Vector{MatElem{nf_elem}}
  s::Tuple{Int64, Int64}

  function ParameterMapIsotypical(f,U,s)
    z = new{}(f,U,s)
    return z
  end
end

size_input(P::ParameterMapIsotypical) = P.s

splitting_field(P::ParameterMapIsotypical) = base_ring(P.U[1])

base_module(P::ParameterMapIsotypical) = P.U

function weighted_sum(P::ParameterMapIsotypical, M::MatElem{nf_elem})
  F = splitting_field(P)
  @req base_ring(M) === F "Wrong parent of input: call `splitting_field(P)`"

  @req size(M) == size_input(P) "Wrong size of input: call `size_input(P)`"

  @req rank(M) == nrows(M) "Input must be of full rank"

  return P.f(M)
end

(P::ParameterMapIsotypical)(M::MatElem{nf_elem}) = weighted_sum(P, M)

function Base.show(io::IO, ::MIME"text/plain", P::ParameterMapIsotypical)
  println(io, "Parametrization of weighted sums for the isotypical module with basis")
  println(io, P.U)
  print(io, "by a full rank matrix of size $(P.s)")
end

function Base.show(io::IO, P::ParameterMapIsotypical)
  print(io, "Parameter map for submodules of an isotypical module.")
end

struct ParameterMap
  f::Function
  U::Vector{Vector{MatElem{nf_elem}}}
  s::Vector{Tuple{Int64, Int64}}
  
  function ParameterMap(f, U, s)
    z = new{}(f, U, s)
    return z
  end
end

size_input(P::ParameterMap) = P.s

splitting_field(P::ParameterMap) = base_ring(P.U[1][1])

base_module(P::ParameterMap) = P.U

function weighted_sum(P::ParameterMap, M::Vector{T}) where T <: MatElem{nf_elem}
    F = splitting_field(P)
  @req all(m -> base_ring(m) === F, M) "Wrong base ring for the input: call `parent_inputs(P)`"

  s = size_input(P)
  @req length(M) == length(s) "Wrong number of inputs: should be $(length(s))"
  @req all(i -> size(M[i]) == s[i], 1:length(s)) "Wrong size of inputs: call `size_inputs(P)"
  
  @req all(m -> rank(m) == nrows(m), M) "Inputs must be of full rank"
  return P.f(M)
end

(P::ParameterMap)(M::Vector{T}) where T <: MatElem{nf_elem} = weighted_sum(P, M)

(P::ParameterMap)(x::Vararg{T}) where T <: MatElem{nf_elem} = P(collect(x))

function Base.show(io::IO, ::MIME"text/plain", P::ParameterMap)
  println(io, "Parametrization of weighted sums for the module with basis")
  println(io, P.U)
  print(io, "by full rank matrices of size $(P.s)")
end

function Base.show(io::IO, P::ParameterMapIsotypical)
  print(io, "Parameter map for submodules of a module.")
end

struct ModSpaceSubmodulesWithIsotypicalCharacter
  chi
  prep_mod::ProjRep
  P::ParameterMapIsotypical
  d::Int64

  function ModSpaceSubmodulesWithIsotypicalCharacter(chi, prep, P, dim)
    z = new{}(chi, prep, P, dim)
    return z
  end
end

submodule_character(M::ModSpaceSubmodulesWithIsotypicalCharacter) = M.chi

submodule_dimension(M::ModSpaceSubmodulesWithIsotypicalCharacter) = M.chi[1]

module_representation(M::ModSpaceSubmodulesWithIsotypicalCharacter) = M.prep_mod

module_character(M::ModSpaceSubmodulesWithIsotypicalCharacter) = char_gap(module_representation(M))

parametrizing_map(M::ModSpaceSubmodulesWithIsotypicalCharacter) = M.P

dim(M::ModSpaceSubmodulesWithIsotypicalCharacter) = M.d

is_irreducible(M::ModSpaceSubmodulesWithIsotypicalCharacter) = true

function Base.show(io::IO, ::MIME"text/plain", M::ModSpaceSubmodulesWithIsotypicalCharacter)
  chi = submodule_character(M)
  println(io, "Moduli space of $(chi[1])-dimensional isotypical submodules of")
  println(io, module_representation(M))
  println(io, "with character")
  print(io, chi)
end

function Base.show(io::IO, M::ModSpaceSubmodulesWithIsotypicalCharacter)
  print(io, "Moduli space of isotypical submodules of dimension $(dim(M))")
end

struct ModSpaceSubmodulesWithCharacter
  chi
  dec::Vector
  P::ParameterMap
  constituents::Vector{ModSpaceSubmodulesWithIsotypicalCharacter}

  function ModSpaceSubmodulesWithCharacter(prod::Vector{ModSpaceSubmodulesWithIsotypicalCharacter}, P)
    dec = [M.chi for M in prod]
    chi = sum(dec)

    z = new{}(chi, dec, P, prod)
    return z
  end
end

submodule_character(M::ModSpaceSubmodulesWithCharacter) = M.chi

submodule_dimension(M::ModSpaceSubmodulesWithCharacter) = M.chi[1]

moduli_factors(M::ModSpaceSubmodulesWithCharacter) = M.constituents

module_representation(M::ModSpaceSubmodulesWithCharacter) = module_representation(moduli_factors(M)[1])

module_character(M::ModSpaceSubmodulesWithCharacter) = char_gap(module_representation(M))

parametrizing_map(M::ModSpaceSubmodulesWithCharacter) = M.P

dim(M::ModSpaceSubmodulesWithCharacter) = sum(dim.(moduli_factors(M)))

is_irreducible(ModSpaceSubmodulesWithCharacter) = true

character_decomposition(M::ModSpaceSubmodulesWithCharacter) = M.dec

function Base.show(io::IO, ::MIME"text/plain", M::ModSpaceSubmodulesWithCharacter)
  chi = submodule_character(M)
  println(io, "Moduli space of $(chi[1])-dimensional submodules of")
  println(io, module_representation(M))
  println(io, "with character")
  print(io, chi)
end

function Base.show(io::IO, M::ModSpaceSubmodulesWithCharacter)
  print(io, "Moduli space of submodules of dimension $(dim(M))")
end

@attributes mutable struct ModSpaceSubmodules
  irr_comp::Vector{ModSpaceSubmodulesWithCharacter}
  prep::ProjRep
  t::Int64

  function ModSpaceSubmodules(p::Vector{ModSpaceSubmodulesWithCharacter}, prep::ProjRep, t::Int64)
    z = new{}(p, prep, t)
    return z
  end
end

submodule_dimension(M::ModSpaceSubmodules) = M.t

irreducible_components(M::ModSpaceSubmodules) = M.irr_comp

is_empty(M::ModSpaceSubmodules) = length(irreducible_components(M)) == 0

is_irreducible(M::ModSpaceSubmodules) = length(irreducible_components(M)) == 1

module_representation(M::ModSpaceSubmodules) = M.prep

module_character(M::ModSpaceSubmodules) = char_gap(module_representation(M))

dim(M::ModSpaceSubmodules) = begin; is_empty(M) && return -1; return maximum(dim.(irreducible_components(M))); end

character_partitions(M::ModSpaceSubmodules) = _partition(module_character(M), submodule_dimension(M))

function irreducible_component(M::ModSpaceSubmodules, chi)
  @req chi in character_partitions(M) "chi is not the character of a module in M"

  j = findfirst(m -> submodule_character(m) == chi, irreducible_components(M))

  return irreducible_components(M)[j]
end

function Base.show(io::IO, ::MIME"text/plain", M::ModSpaceSubmodules)
  println(io, "Moduli space of $(submodule_dimension(M))-dimensional submodules of")
  print(io, module_representation(M))
end

function Base.show(io::IO, M::ModSpaceSubmodulesWithCharacter)
  print(io, "Moduli space of submodules of dimension $(dim(M))")
end

function _parametrising_isotypical_submodule(prep::ProjRep,B, chi)
  RR = representation_ring(prep)
  ct = character_table(RR)
  cd = _character_decomposition(ct, chi)
  @assert length(cd) == 1
  alpha, chis = cd[1]
  U = _basis_simple_submodule(prep, B, chis)
  U = lcm(denominator.(U)) * U
  n = chis[1]
  d = divexact(ncols(U), n)
  UU = MatElem{nf_elem}[U[:,1+(i-1)*n:i*n] for i in 1:d]

  function _basis_parameterisation(M::MatElem{nf_elem})
    F = base_ring(M)
    V = zero_matrix(F, nrows(UU[1]), alpha*ncols(UU[1]))
    for i=1:nrows(M)
      VV = zero_matrix(F, nrows(UU[1]), ncols(UU[1]))
      for j=1:ncols(M)
        VV += M[i,j]*UU[j]
      end
      V[:, 1+(i-1)*n:i*n] = VV
    end
    return V::MatElem{nf_elem}
  end
  
  return ParameterMapIsotypical(_basis_parameterisation, UU, (alpha, d))
end

function _parametrising_submodule(prep::ProjRep, B, chi)
  RR = representation_ring(prep)
  ct = character_table(RR)
  cd = _character_decomposition(ct, chi)
  chis = [c[1]*c[2] for c in cd]
  Ub = Vector{MatElem{nf_elem}}[]
  Uf = Function[]
  Us = Tuple{Int64, Int64}[]
  for c in chis
    P = _parametrising_isotypical_submodule(prep, B, c)
    push!(Ub, base_module(P))
    push!(Uf, P.f)
    push!(Us, size_input(P))
  end

  function _basis_parameterisation2(x::Vector{T}) where T <: MatElem{nf_elem}
    V = T[]
    for i in 1:length(Uf)
      f = Uf[i]
      push!(V, f(x[i]))
    end
    return V::Vector{T}
  end

  return ParameterMap(_basis_parameterisation2, Ub, Us)
end


function moduli_space_of_submodules_with_isotypical_character(prep::ProjRep, K::MatElem{nf_elem}, chi)
  ct = character_table(representation_ring(prep))
  nu = char_gap(prep)
  @req length(_character_decomposition(ct, chi)) == 1 "chi is not isotypical"
  @req GG.IsCharacter(ct, nu-chi) "chi is not a constituent of the character of prep"
  f = _parametrising_isotypical_submodule(prep, K, chi)
  t, n = size_input(f)
  d = t*(n-t)
  return ModSpaceSubmodulesWithIsotypicalCharacter(chi, prep, f, d)
end

function moduli_space_of_submodules_with_character(prep::ProjRep, K::MatElem{nf_elem}, chi)
  ct = character_table(representation_ring(prep))
  nu = char_gap(prep)
  @req GG.IsCharacter(ct, nu-chi) "chi is not a constituent of the character of prep"
  f = _parametrising_submodule(prep, K, chi)
  dec = _character_decomposition(ct, chi)
  constituents = ModSpaceSubmodulesWithIsotypicalCharacter[]
  P = _parametrising_submodule(prep, K ,chi)
  for d in  dec
    M = moduli_space_of_submodules_with_isotypical_character(prep, K, d[1]*d[2])
    push!(constituents, M)
  end

  return ModSpaceSubmodulesWithCharacter(constituents, P)
end

function moduli_space_of_submodules(prep::ProjRep, t::Int64)
  @req 1 <= t < dim(prep) "t must be positive and (strictly) smaller than the dimension of prep"
  p2, K = _to_equivalent_block_representation(prep)
  @info "Equivalent block representation and base change computed"
  ct = character_table(representation_ring(prep))
  chis = _partition(p2, t)
  irr_comp = ModSpaceSubmodulesWithCharacter[]
  @info "Construct $(length(chis)) irreducible component(s)"
  for chi in chis
    M = moduli_space_of_submodules_with_character(prep, K, chi)
    push!(irr_comp, M)
  end

  return ModSpaceSubmodules(irr_comp, prep, t)
end

function moduli_spaces_generators_ideals(id::Tuple{Int64, Int64}, n::Int64, d::Int64, t::Int64)
  pfr = projectively_faithful_representations(id[1], id[2], n)
  @info "Compute $(length(pfr)) moduli space(s)"
  res = ModSpaceSubmodules[]
  for prep in pfr
    prep_homo = homo_poly_rep(prep, d)
    M = moduli_space_of_submodules(prep_homo, t)
    push!(res, M)
    @info "done"
  end
  return pfr, res
end

##############################################################################
#
# Fitting ideals (with authorization)
#
##############################################################################

"""
  An internal based on a code of Matthias Zach (TU Kaiserslautern) for computation
  of fitting ideals. Here it is used to look for the ideal of coordinates of pure
  tensors (computed from the denegeracy locus for the rank of a sparse matrix) and
  the smoothness test.
"""
function _ideal_degeneracy(X, Mfact, k)
  isempty(X) ? (return [base_ring(X)]) : nothing
  
  Mfact = _refact(Mfact)
  val, row, col = Mfact

  if k == 1
    if length(val) == 0  
      return [X]
    else
      Y = subscheme(X, ideal(OO(X), val))
      isempty(X) ? (return [base_ring(Y)]) : (return [Y])
    end
  end
  
  if length(val) == 0
    return [X]
  end
  
  d = maximum(total_degree.(lifted_numerator.(val)))
  (i, j) = (1, 1)

  for l=1:length(val)
    dv = total_degree(lifted_numerator(val[l]))
    if dv <= d && !iszero(val[l])
      d = dv
      (i, j) = (row[l], col[l])
    end
  end
  
  v = _value(Mfact, i, j)

  if iszero(val)
    return [X]
  end
  
  val2, row2, col2 = typeof(val[1])[], [], []
  for l in union(1:(i-1), (i+1):maximum(row))
    for m=1:maximum(col)
      vall = _value(Mfact, l, m)*v - _value(Mfact, i, m)*_value(Mfact, l, j)
      if !iszero(vall)
        push!(val2, vall)
	      l < i ?	push!(row2, l) : push!(row2, l-1)
	      m < j ? push!(col2, m) : push!(col2, m-1)
      end
    end
  end
    
  Y = subscheme(X, v)
  Mfact = (OO(Y).(Oscar.lift.(val)), row, col)

  U = hypersurface_complement(X, v)
  Mfact2 = (OO(U).(Oscar.lift.(val2)), row2, col2)

  ID = _ideal_degeneracy(Y, Mfact, k)
  _ID2 = _ideal_degeneracy(U, Mfact2, k-1)
  ID2 = [closure(V, X) for V in _ID2 if V != base_ring(X)]

  ID = [V for V in ID if V != base_ring(X) && !isempty(V)]
  filter!(V -> !all(xi-> radical_membership(xi, defining_ideal(V)), gens(ambient_ring(V))), ID) 
  ID2 = [V for V in ID2 if V != base_ring(X) && !isempty(V)]
  filter!(V -> !all(xi-> radical_membership(xi, defining_ideal(V)), gens(ambient_ring(V))), ID2) 
  return union(ID,ID2)
end

##############################################################################
#
# some checks functions
#
##############################################################################

### Smoothness test

"""
   is`_`smooth(I::T) where T <: MPolyIdeal -> Bool

Given an homogeneous ideal `I`, return `true` if `V(I)` is smooth, `false` else.
"""
function is_smooth(I::T) where T <:MPolyIdeal
  v = [f for f in gens(I)]
  R = base_ring(I)
  L = quo(R,I)[1]
  X = Spec(L)
  x = gens(R)
  D = jacobi_matrix(v)
  D = map_entries(OO(X), D)
  Dfact = _sparse_matrix_fact(D)
  k = length(x) - Oscar.dim(L)
  ID = _ideal_degeneracy(X, Dfact, k)
  length(ID) == 0 ? (return true) : nothing
  for Y in ID
    bool = true
    for y in x
      U = hypersurface_complement(Y, y)
      bool = isempty(U) ? bool : false
    end
    bool == false ? (return bool) : nothing
  end
  return true
end

### Invariance test

function _action_on_basis(MG::MatrixGroup, M::Matrix)
  @assert Oscar.degree(MG) == nrows(M)
  F = base_ring(matrix(M))
  coll = MatElem[]
  for m in gens(MG)
    N = solve(matrix(M), m*matrix(M))
    push!(coll, N)
  end
  return matrix_group(coll)
end

function _action_on_basis(prep_homo::ProjRep, I::MPolyIdeal)
  MG = matrix_group(matrix_representation(prep_homo))
  M, _ = _homo_ideal_in_coord(I)
  return _action_on_basis(MG, M)
end

"""
  An internal to check if an ideal is invariant under a group action given by `rep_homo`.
  The generators of `I` has to be of the same degree
"""
function _is_invariant_by_rep(prep::ProjRep, I::T) where {S,T}
  R = base_ring(I)
  @req R == polynomial_algebra(representation_ring(prep)) "polynomial algebra must coincide"
  d = total_degree(gens(I)[1])
  prep_homo = get_attribute(prep, :prep_homo)
  if prep_homo !== nothing
    d0 = get_attribute(prep_homo, :degree_homo)
  end
  if (prep_homo === nothing || d != d0)
    prep_homo = homo_poly_rep(prep, d)
  end
  F = base_ring(R)
  Coord,_ = _homo_ideal_in_coord(I)
  L,y = PolynomialRing(F,"y" => (1:(ncols(Coord)+1)), cached = false)
  Rep_homo = matrix_representation(prep_homo)
  for M in Rep_homo
    for j=1:ncols(Coord)
      act = change_base_ring(L,M*matrix(Coord[:,j]))
      v = y[end]*act-sum([y[i]*change_base_ring(L,matrix(Coord[:,i])) for i=1:ncols(Coord)])
      II = ideal(L,[v[i,:][1] for i=1:nrows(v)])
      J = radical(II)
      J == ideal(L,y) || J == ideal(L,L(1)) ? (return false) : continue
    end
  end
  return true
end

### About (non)-symplectic actions

# This Lemma (2.1) of Mukai '88 "Finite groups of automorphisms of K3 surfaces and the 
# Mathieu group: only works under the hypotheses of the Lemma!
function _induced_action_symplectic_form(prep::ProjRep, I)
  R = base_ring(I)
  @req R == polynomial_algebra(representation_ring(prep)) "polynomial algebra must coincide"
  d = total_degree(gens(I)[1])
  F = base_ring(R)
  Coord, B = _homo_ideal_in_coord(I)
  L,y = PolynomialRing(F,"y" => (1:(ncols(Coord)+1)), cached = false)
  Rep = matrix_representation(prep)
  MG = matrix_group(Rep)
  coll = elem_type(F)[]
  for m in gens(MG)
    PM = matrix(m)
    M = _homo_poly_action(PM, R, d) 
    MM = matrix(zeros(F,ncols(Coord),ncols(Coord)))
    VV = VectorSpace(F,length(B))
    Var,f = sub(VV, [VV(Coord[:,i]) for i=1:ncols(MM)])
    Pas = hcat([transpose((f\VV(transpose(matrix(Coord[:,j])))).v) for j=1:ncols(MM)])
    for j=1:ncols(MM)
      MM[:,j] = inv(Pas)*transpose((f\VV(transpose(M*matrix(Coord[:,j])))).v)
    end
    l = F(det(PM))*F(inv(det(MM)))
    push!(coll, l)
  end
  GF = matrix_group(matrix(F,1,1, [gen(F)]))
  collF = [GF([a]) for a in coll]
  f = hom(MG, GF, collF)
  return f
end
  
  
"""
  An internal to check if the action `rep` of a group on `V(I)` is symplectic.
"""
function _is_symplectic_action(prep::ProjRep, I) 
  f = _induced_action_symplectic_form(prep, I)
  return (kernel(f)[1] == domain(f))
end

"""
  An internak to check whether the action `rep` of a group on `V(I)` is purely 
  non-symplectic.
"""
function _is_purely_non_symplectic_action(prep::ProjRep, I)
  f = _induced_action_symplectic_form(prep, I)
  return is_trivial(kernel(f)[1])
end

"""
  Given `rep` acting on `V(I)`, return the id of the group generated by the
  elements of `rep` acting symplectically on `V(I)`.
"""
function _id_symplectic_subgroup(prep::ProjRep, I)
  f = _induced_action_symplectic_form(prep, I)
  MG = domain(f)
  MG2, MG2toMG = kernel(f)
  Q, p = quo(MG, MG2)
  @assert is_cyclic(Q)
  elem = collect(Q)
  j = findfirst(zz -> order(zz) == order(Q), elem)
  zMG2 = preimage(p, elem[j])
  return MG2, zMG2
end

######################################################################
#
#  Sparse matrices
#
######################################################################

"""
  An internal to reduce a polynomial sparse matrix by deleting all zero rows and 
  columns, and deleting rows and columns appearing several times.
"""
function _reduction_sparse_mat(M::T) where T 
  @assert !iszero(M)
  M = vcat([M[i,:] for i =1:nrows(M) if !iszero(M[i,:])])
  M = hcat([M[:,j] for j =1:ncols(M) if !iszero(M[:,j])])
  rowM = [M[1,:]]
  for i=2:nrows(M)
    rowint = M[i,:]
    bool = !any(roww -> (rowint == roww || rowint == -roww), rowM)
    bool == true ? push!(rowM, rowint) : continue
  end
  M = vcat(rowM)
  colM = [M[:,1]]
  for j=2:ncols(M)
    colint = M[:,j]
    bool = !any(coll -> (colint == coll || colint == -coll),  colM)
    bool == true ? push!(colM, colint) : continue
  end
  return hcat(colM)
end

"""
  An internal to store sparse matrices into 3 lists: a list of non-zero values, and
  the other lists that keep track of the indices of the rows and columns in which each
  value lies.
"""
function _sparse_matrix_fact(M::T) where T  
  @assert !iszero(M)
  M = _reduction_sparse_mat(M)
  values = [M[i,j] for j=1:ncols(M), i=1:nrows(M) if !iszero(M[i,j])]
  indrow = [i for j=1:ncols(M), i=1:nrows(M) if !iszero(M[i,j])]
  indcol = [j for j=1:ncols(M), i=1:nrows(M) if !iszero(M[i,j])]
  return (values, indrow, indcol)
end

"""
  An internal that given indices `i,j` returns the value in the sparse matrix,
  at the `i`th row and `j`th column.
"""
function _value(M, i,j)
  val,row,col = M
  k = findfirst(k -> row[k] == i && col[k] == j, 1:length(val))
  k == nothing ? (return parent(val[1])(0)) : (return val[k])
end

"""
  An internal to refactorise a sparse matrix (deleting the zero entries
  in its factorised form).
"""
function _refact(Mfact)
  val, row, col = Mfact
  if length(val) == 0
    return [],[],[]
  end
  valbis, rowbis, colbis = typeof(val[1])[], typeof(row[1])[], typeof(col[1])[]
  for i=1:length(val)
    if !iszero(val[i])
      push!(valbis,val[i])
      push!(rowbis, row[i])
      push!(colbis, col[i])
    end
  end
  return (valbis, rowbis, colbis)
end

######################################################################
#
#  Intersections with Grassmannians
#
######################################################################

"""
  An internal to compute intersection with Grassmannian varieties in dimension bigger than 
  2. For now, it returns the ideals parametrising the coordinates of pure tensors. It 
  still need the user to construct the (family.ies of) pure tensors 'by hand'
"""
function _intersection_with_grassmannian(V, basis)
  R = parent(basis[1][1])
  d = total_degree(basis[1][1])
  B = homog_basis(R,d)
  k = length(basis[1])
  M = _matrix_multivec_star(V, B, k)
  L = parent(M[1,1])
  Q = quo(L, ideal(L, [L(0)]))[1]
  X = Spec(Q)
  M = map_entries(OO(X), M)
  Mfact = _sparse_matrix_fact(M)
  ID = _ideal_degeneracy(X, Mfact, k+1)
  return modulus.(OO.(ID))
end

### Computations of isotypic components in exterior power

"""
  An internal function that returns the intersection of an isotypic component of weight 1 of
  the `t`-exterior power `W` of the `d`-homogeneous part `V` of a polynomial algebra with
  the Grassmannian `Gr(t,V). `W` is seen as a group algebra via `rep`. The input 'smooth' is
  to allow only the elements in the intersection generating an ideal defining a smooth variety.
"""
function _defining_ideal_moduli_determinantal_character(prep_homo::ProjRep, chi::Vector, t::Int64)
  R = polynomial_algebra(representation_ring(prep_homo))
  d = get_attribute(prep_homo, :degree_homo)
  basis = basis_exterior_power(homog_basis(R,d), t)
  prep_ext = t ==1 ? prep_homo : ext_homo_poly_rep(prep_homo, t)
  rep = matrix_representation(prep_ext)
  rep = [rep[i]-chi[i][1]*identity_matrix(base_ring(R), ncols(rep[1])) for i=1:length(rep)]
  rep = reduce(vcat, rep)
  srep = sparse_matrix(rep)
  dim, V = right_kernel(srep)
  if dim == 0
    return false
  end
  V = [transpose(V[:,i]) for i=1:ncols(V)]
  if t == 1
    return V
  end
  ID = _intersection_with_grassmannian(V, basis_exterior_power(basis, t))
  gene = gens(ID[1])
  for i=2:length(ID)
    append!(gene, gens(ID[i]))
  end
  
  return radical(ideal(R, gene))
end