################################################################
#
#
# Condensation for Tensor Product 10/2/2015
#
# Updated 8/12/16
#
##############################################################

#packages that need to be loaded
LoadPackage("basic");
LoadPackage("cond");
Read(Concatenation(PackageInfo("cond")[1].InstallationPath, "/gap/matching.g"));
#this next line can be removed if homcond1 gets changed
#Read("/u6/akraft/homcond2");
#Read("/u4/klux/mypkg/pkg/basic/devKlaus/condensationmethods/homcond/homcond1.g");
Read(Concatenation(PackageInfo("basic")[1].InstallationPath,
           "/devKlaus/condensationmethods/homcond/homcond1.g"));

#change the base field of the representations, cgens should be the generators for the represenation (as cmats) and field should be the field you want the generators to be over

ChangeField := function(cgens,field)
local i,mats;
#first check that the cgens lie in a subfield of field
if IsSubspace(field, BaseField(cgens[1])) = false then
	return "Invalid field";
fi;
mats := [];
for i in [1..Size(cgens)] do
	Add(mats,CMat(List(cgens[i],x -> CVec(Unpack(x),field))));
od;
return mats;
end;

LcmList := function(list)
local l, lcm, i,j, counter;
l := Collected(list);
#check if all things in list are 1
if Length(l) = 1 then
	return 1;
fi;
#don't need to take lcm of (1,n)
if l[1][1] = 1 then
	Remove(l,1);
fi;
#keep taking pairwise lcm until you get one number
counter := Length(l);
while counter > 1 do
	lcm := [];
	for i in [1..Length(l)-1] do
		for j in [i+1..Length(l)] do
			Add(lcm, LcmInt(l[i][1],l[j][1]));
		od;
	od;
	l := Collected(lcm);
	counter := Length(l);
od;
return l[1][1];
end;

#find splitting field for representations of the group
MinField := function(G,p)
local list,powers,ct, irr, i,j;
list := [];
powers := [];
ct := CharacterTable(G) mod p;
irr := Irr(ct);
for i in [1..Length(irr)] do
	Add(list, List(irr[i],x->FrobeniusCharacterValue(x,p)));
od;
for j in [1..Length(list)] do
	Add(powers, LogInt(Size(Field(list[j])),p));
od;
return GF(p^LcmList(powers));
end;

#extract the diagonal blocks of a set of upper/lower block diagonal matrices
#sub should be a list of upper/lower block diagonal matrices
#dimlist should be an increasing list of numbers the difference between consecutive terms corresponding to the sizes of each block


GetSubreps := function(sub,dimlist)
local i, d, submods;
#since dimlist is not necessarily strictly increasing
d := Collected(dimlist);
submods := [];
submods[1] := List(sub,x -> ExtractSubMatrix(x,[1..d[1][1]],[1..d[1][1]]));
if Length(d) > 1 then
for i in [2..Length(d)] do
	submods[i] := List(sub, x-> ExtractSubMatrix(x,[d[i-1][1]+1..d[i][1]],[d[i-1][1]+1..d[i][1]]));
od;
fi;
return submods;
end;

#Creates a basis vector of length dim in matrix form with 1 in kth position (over the field GF(p^d))

BasisVec := function(p,d,dim,k)
local v;
v := CVec(CVecClass(p,d,dim));
v[k] := 1;
return v;
end;

#computes the Brauer character of an absolutely irreducible representation of a group g for a given p and generating cmats gens

CharOfRep := function(g,p,gens)
local gp,a,c,d,traces,ct,pos;
gp :=AllAtlasGeneratingSetInfos(g);
#create conjugacy classes
a:= AtlasStraightLineProgram(g,"classes");
	if a = fail then
		a := LocalClasses(g);
			if a = fail then
				return fail;
			fi;
	fi;
c:= ResultOfStraightLineProgram(a.program, gens);
#restrict to certain conjugacy classes for Brauer characters
d:= Filtered(c, x -> Order(x) mod p > 0);
#exception for trivial representation
if Length(c) = Length(d) then
	return "1a";
fi;
traces:= List(d, x -> BrauerCharacterValue(Unpack(x)));
ct := CharacterTable(g) mod p;
pos := Position(Irr(ct), traces);
#in case the representation is not absolutely irreducible
if pos = fail then 
	return fail;
else
	return AtlasCharacterNames(ct)[pos];
fi;
end;

#Klaus's code for creating an empty basis
EmptyBasis := function(conddim,alg)
local vec, vec1, mvec, mvec1;
vec:=ListWithIdenticalEntries(conddim,0*Z(alg.prime));
ConvertToVectorRep(vec);
vec1:=CVec(vec);
vec1[1]:=Z(alg.prime)^0;
# make a copy to construct empty basis
#
mvec:=[ShallowCopy(vec1)];
mvec1:=CMat(mvec);
return EmptySemiEchelonBasis(mvec1);
end;

#The main code. G should be a string, p a prime, m,n the position of the desired representations in AllAtlasGeneratingSetsInfos.

TensorRep := function(alg,m,n)
local G, p, ag, mats1, mats2, res, M, decomp, sbasis, basis, i, k, vec,decM,dimlist, sub, diff, bmats1, bmats2, submods, chars, mult, gens, val, lastirred;

#check if condensation subgroup is trivial
if alg.condinfo.horder = 1 then
	return "Condensation subgroup is trivial.";
fi;

G := alg.group;
p := alg.prime;
ag :=AllAtlasGeneratingSetInfos(G);
#get the generators as cmats
mats1 := ChangeField(List(AtlasGenerators(ag[m].identifier).generators, x-> CMat(x)),MinField(G,p));
mats2 := ChangeField(List(AtlasGenerators(ag[n].identifier).generators, x-> CMat(x)),MinField(G,p));
#alg := InitializeRecordAtlasGroup(G,p,1);
#AutoCalcBasic(alg);

res := homcond(alg,mats1,mats2);

#check if the the condensed tensor is the zero tensor
if res.info.Dim = 0 then
	return "The condensed tensor is 0 tensor. Use a different block";
fi;

#need to add information to the record
res.info.Mdim:=ag[m].dim;
res.info.Ndim:=ag[n].dim;
#rewrite the generating matrices for the modules with respect to the new bases. The first one needs transpose inverse
bmats1 := List(mats1, x-> res.rtbasis*TransposedMat(x)^-1*res.rtbasis^-1);
bmats2 := List(mats2, x-> res.rbasis*x*res.rbasis^-1);
M := Module(res.condgens);
decM := Chop(M,rec(compbasis := true));
sbasis := [];
#uncondense
for i in [1..Size(decM.basis)] do
	Add(sbasis, HomUncond(decM.basis[i],res.info));
od;
#create empty basis
vec := CH0P_Unfurl(sbasis[1]);
basis := EmptySemiEchelonBasis(CMat([vec],vec));
#spin the vectors, keeping track of dimensions
dimlist := [];
for i in [1..Size(sbasis)] do
	CH0P_TensorSpin(CH0P_Unfurl(sbasis[i]),List(bmats1,x->x^-1),bmats2,basis,ag[m].dim*ag[n].dim);
	Add(dimlist,Size(basis));
od;
#need to make sure the vectors span the entire space, and if not, add pivots to get basis
if Size(basis) < ag[m].dim*ag[n].dim then
	lastirred := false;
	Add(dimlist,ag[m].dim*ag[n].dim);
	#Add extra pivots to get full basis
	diff := Difference([1..ag[n].dim*ag[m].dim],Pivots(basis));
	k := Log(Size(Field(basis[1])),p);
	for i in diff do
	CH0P_TensorSpin(BasisVec(p,k,ag[n].dim*ag[m].dim,i),
	List(bmats1, x-> x^-1),bmats2,basis,ag[n].dim*ag[m].dim);
	od;
	else lastirred := true;
fi;
#get generating matrices for the tensor product in almost block diagonal form
sub:=CH0P_TensorActionOnSubspace(List(bmats1,x->x^-1),
        bmats2,basis);

#extract all the subreps:
submods := GetSubreps(sub,dimlist);

#decompose the last block if its not already irreducible
if lastirred = false then
	decomp := Chop(Module(submods[Length(submods)]));
	Remove(submods);
fi;

#Find characters, multiplicities of each rep
chars := [];
mult := [];
gens := [];

for i in [1..Length(submods)] do
	val := CharOfRep(G,p,submods[i]);
	if Position(chars,val) = fail then
		Add(chars,val);
		Add(gens,submods[i]);
		Add(mult,1);
	else mult[Position(chars,val)] := mult[Position(chars,val)] +1;
	fi;
od;

if lastirred = false then
	for i in [1..Length(decomp.db)] do
		val := CharOfRep(G,p,
			RepresentingMatrices(decomp.db[i]));
		if Position(chars,val) = fail then
			Add(chars,val);
			Add(gens,decomp.db[i]);
			Add(mult,1);
		else mult[Position(chars,val)] := mult[Position(chars,val)] +1;
		fi;
	od;
fi;

#output is a record containing
#1. A list of generators for the irreducible subreps of M tensor N
#2. The characters of each of these representations
#3. The multiplicities of each of these representations
return rec(characters := chars, multiplicities := mult, generators := gens);
end;