#(C) 2009 Graham Ellis

###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfQuotientGroupSeries_Int,
function(arg)
		#S is a chain of normal subgroups with S[1]=G.   
local
	S,deg,Resol,RESLST,HOMS,A,B,
	N,P,Q,RP,RQ,T,Persistence,hom,eqmap,map,i,j;

#####INPUT#################
S:=arg[1]; #A normal series of subgroups
deg:=arg[2]; #The degree in which homology will be calculated.
if Length(arg)>2 then Resol:=arg[3]; else 
Resol:=ResolutionFiniteGroup;
fi;
###########################

HOMS:=NormalSeriesToQuotientHomomorphisms(S);;
if Order(Range(HOMS[Length(HOMS)]))=1 then
HOMS:=HOMS{[1..Length(HOMS)-1]};
fi;
RESLST:=[];

RP:=Resol(Source(HOMS[1]),deg+1);
Add(RESLST,RP);

for i in [1..Length(HOMS)] do
RQ:=Resol(Range(HOMS[i]),deg+1);
Add(RESLST,RQ);
od;

Persistence:=List( [1..Length(HOMS)+1], i->List([1..Length(HOMS)+1],j->0));

for i in [1..Length(HOMS)] do
for j in [i..Length(HOMS)] do
if i=j then hom:=IdentityMapping(Source(HOMS[i]));
else
hom:=Product(HOMS{[i..j-1]});
fi;
eqmap:=EquivariantChainMap(RESLST[i],RESLST[j],hom);
map:=TensorWithIntegers(eqmap);
map:=Homology(map,deg);
#map:=Image(map);
if Size(Image(map))=1 then A:=[];
else A:=AbelianInvariants(Image(map));
fi;
if Size(Range(map)/Image(map))=1 then B:=[];
else B:=AbelianInvariants(Range(map)/Image(map));
fi;

Persistence[i][j]:=[A,B];
od;
od;

return Persistence;

end);
###############################################################
###############################################################

###############################################################
###############################################################
InstallGlobalFunction(PersistentCohomologyOfQuotientGroupSeries,
function(arg)
                #S is a chain of normal subgroups with S[1]=G.
local
        S,deg,Resol,RESLST,HOMS,A,B,
        N,P,Q,RP,RQ,T,prime,Persistence,hom,eqmap,map,LN,i,j;

#####INPUT#################
S:=arg[1]; #A normal series of subgroups
deg:=arg[2]; #The degree in which homology will be calculated.
if Length(arg)>2 then 
prime:=arg[3];
else prime:=PrimePGroup(S[1]); fi;
if Length(arg)>3 then
Resol:=arg[4]; else
  if not prime=0 then
  Resol:=ResolutionPrimePowerGroup;
  else
  Resol:=ResolutionFiniteGroup;
  fi;
fi;
###########################


HOMS:=NormalSeriesToQuotientHomomorphisms(S);;

if Order(Range(HOMS[Length(HOMS)]))=1 then
#HOMS:=HOMS{[1..Length(HOMS)-1]};
fi;
RESLST:=[];

RP:=Resol(Source(HOMS[1]),deg+1);
Add(RESLST,RP);

for i in [1..Length(HOMS)] do
RQ:=Resol(Range(HOMS[i]),deg+1);
Add(RESLST,RQ);
od;

LN:=Length(HOMS)+1;
Persistence:=List( [1..LN], i->List([1..LN],j->0));

for i in [1..LN] do
for j in [i..LN] do
if i=j  then 
        if i=LN then
	hom:=IdentityMapping(Image(HOMS[LN-1]));
        else
	hom:=IdentityMapping(Source(HOMS[i]));
	fi;
else
 	hom:=Product(HOMS{[i..j-1]});

fi;
eqmap:=EquivariantChainMap(RESLST[i],RESLST[j],hom);
if not prime=0 then
map:=HomToIntegersModP(eqmap,prime);
map:=Cohomology(map,deg);
Persistence[i][j]:=Log(Order(Image(map)),prime);

else

map:=HomToIntegers(eqmap);
map:=Cohomology(map,deg);
A:=Image(map);
B:=Range(map)/A;
  if Order(A)>1 then A:=AbelianInvariants(A); else A:=[]; fi;
  if Order(B)>1 then B:=AbelianInvariants(B); else B:=[]; fi;
Persistence[i][LN-j]:=[A,B];

fi;
od;
od;

Persistence:=Persistence{[2..Length(Persistence)]};
Persistence:=TransposedMat(Persistence);
Persistence:=Persistence{[2..Length(Persistence)]};
Persistence:=TransposedMat(Persistence);

return Persistence;

end);
###############################################################
###############################################################


###############################################################
###############################################################
InstallGlobalFunction(NormalSeriesToQuotientHomomorphisms,
function(SS)
local
        S,
        HOMS,N,P,G,iso,T,hom,hom1,i;

if Size(SS[1])<=Size(SS[Length(SS)])
then S:=SS;
else
S:=Reversed(SS);
fi;

T:=[];
for i in [1..Length(S)] do
T[i]:=NaturalHomomorphismByNormalSubgroup(S[Length(S)],S[i]);
od;

HOMS:=[T[1]];
for i in [1..Length(S)-1] do
P:=Image(T[i]);
N:=List(GeneratorsOfGroup(S[i+1]),x->Image(T[i],x));
Add(N,Identity(P));
N:=Subgroup(P,N);
hom:=NaturalHomomorphismByNormalSubgroup(P,N);
iso:=IsomorphismGroups(Image(hom),Image(T[i+1]));
hom1:=GroupHomomorphismByImages(P,Image(T[i+1]),
GeneratorsOfGroup(P),
List(GeneratorsOfGroup(P), x-> Image(iso,Image(hom,x))));
Add(HOMS,hom1);
od;

return HOMS;

end);
###############################################################
###############################################################


###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfQuotientGroupSeries,
function(arg)
                #S is a chain of normal subgroups with S[1]=G.
local
        S,deg,prime,Resol,
        HOMS,RP,RQ,T,Persistence,MAPS,eqmap,map,i;

#####INPUT#################
S:=arg[1]; #A normal series of subgroups
deg:=arg[2]; #The degree in which homology will be calculated.
if Length(arg)>2 then
prime:=arg[3]; else prime:=PrimePGroup(S[1]);
fi;
if Length(arg)>3 then Resol:=arg[4]; else
Resol:=ResolutionPrimePowerGroup;
fi;
###########################

if prime=0 then
   if Length(arg)>3 then
   return PersistentHomologyOfQuotientGroupSeries_Int(S,deg,Resol);
   else
   return PersistentHomologyOfQuotientGroupSeries_Int(S,deg);
   fi;
fi;

HOMS:=NormalSeriesToQuotientHomomorphisms(S){[1..Length(S)-1]};
MAPS:=[];

if Length(HOMS)=1 and Resol=ResolutionPrimePowerGroup then

  RP:=Resol(Source(HOMS[1]),deg);
  MAPS:=[
  IdentityMapping(GF(prime)^RP!.dimension(deg))
  ];

else

  RP:=Resol(Source(HOMS[1]),deg+1);
  MAPS:=[
  IdentityMapping(HomologyVectorSpace(TensorWithIntegersModP(RP,prime),deg))
  ];

fi;

for i in [2..Length(HOMS)] do
RQ:=Resol(Image(HOMS[i]),deg+1);
eqmap:=EquivariantChainMap(RP,RQ,HOMS[i]);
map:=TensorWithIntegersModP(eqmap,prime);
map:=HomologyVectorSpace(map,deg);
Add(MAPS,map);
RP:=RQ;
od;

MAPS:=LinearHomomorphismsPersistenceMat(MAPS);
MAPS:=MAPS{[2..Length(MAPS)]};
MAPS:=TransposedMat(MAPS);
MAPS:=MAPS{[2..Length(MAPS)]};
MAPS:=TransposedMat(MAPS);

return MAPS;
end);
###############################################################
###############################################################


###############################################################
###############################################################
InstallGlobalFunction(LinearHomomorphismsPersistenceMat,
function(HOMS)
local
	mat,x,y,L;

L:=Length(HOMS)+1;
mat:=IdentityMat(L)*0;

for x in [1..L] do
for y in [x..L] do
if y=x then
  if x=L then
  mat[x][y]:=IdentityMapping(Range(HOMS[x-1]));
  else
  mat[x][y]:=IdentityMapping(Source(HOMS[x]));
  fi;
else
mat[x][y]:= Product(List([x..y-1],i->HOMS[i]));
fi;
mat[x][y]:=Dimension(Image(mat[x][y]));
od;
od;

return mat;

end);
###############################################################
###############################################################

###############################################################
###############################################################
InstallGlobalFunction(LinearHomomorphismsZZPersistenceMat,
function(HOMS)
local
        L, WrongComposition;

###################################
WrongComposition:=function(LUV,LWV)  #LUV:U-->V, LWV:W-->V
local x, U,W, LUW, ImLUV, ImLWV, InterUW,  BW,  BU, BUW, dim;

ImLUV:=Image(LUV);
ImLWV:=Image(LWV);
InterUW:=Intersection(ImLUV,ImLWV);
U:=Source(LUV);
W:=Source(LWV);
if Dimension(InterUW)=0 then 
   return ZeroMapping(U,W);  

else

   BUW:=Basis(InterUW);
   BUW:=Elements(BUW);
fi;

dim:=Length(BUW);
BU:=List(BUW,x->PreImagesRepresentative(LUV,x));;
BW:=List(BUW,x->PreImagesRepresentative(LWV,x));;


for x in Elements(Basis(U)) do
  if Rank(Concatenation(BU,[x]))>dim then
    dim:=dim+1;
    Add(BU,x);;
    Add(BW,Zero(W));
  fi;
od;

LUW:=LeftModuleHomomorphismByImagesNC(U,W,BU,BW);;
return LUW;
end;
###################################

L:=List([1..Length(HOMS)/2],n->
WrongComposition(HOMS[2*n-1],HOMS[2*n]));

return LinearHomomorphismsPersistenceMat(L); 

end);
###############################################################
###############################################################


###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfSubGroupSeries,
function(arg)
local
	S,deg,prime,hom,HOMS,MAPS,RP,RQ,Resol,eqmap,map,i;

#####INPUT#################
S:=arg[1];
if Order(S[1])>Order(S[Length(S)]) then
S:=Reversed(S); 
fi; #A normal series of subgroups
if Order(S[1])=1 then S:=S{[2..Length(S)]};fi;
deg:=arg[2]; #The degree in which homology will be calculated.
if Length(arg)>2 then
prime:=arg[3]; else prime:=PrimePGroup(S[1]);
fi;
if Length(arg)>3 then Resol:=arg[4]; 
else Resol:=ResolutionPrimePowerGroup;
fi;
###########################

Add(S,S[Length(S)]);
HOMS:=[];

for i in [1..Length(S)-1] do
hom:=GroupHomomorphismByFunction(S[i],S[i+1],x->x);;
Add(HOMS,hom);
od;


MAPS:=[];

RP:=Resol(Source(HOMS[1]),deg+1);
for i in [1..Length(HOMS)] do
RQ:=Resol(Image(HOMS[i]),deg+1);
eqmap:=EquivariantChainMap(RP,RQ,HOMS[i]);
map:=TensorWithIntegersModP(eqmap,prime);
map:=HomologyVectorSpace(map,deg);
Add(MAPS,map);
RP:=RQ;
od;


MAPS:=LinearHomomorphismsPersistenceMat(MAPS);
MAPS:=MAPS{[2..Length(MAPS)]};
MAPS:=TransposedMat(MAPS);
MAPS:=MAPS{[2..Length(MAPS)]};
MAPS:=TransposedMat(MAPS);

return MAPS;

end);
###############################################################
###############################################################


###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfPureCubicalComplex,
function(arg)
local
	LM,deg,prime,
        M,triv,ChnCmps,vsmaps,A,S,TS,TM,L,
        bool,x,y,i,F,PC;

##Input###################
LM:=arg[1];
deg:=arg[2];
prime:=arg[3];
##########################

F:=HomotopyEquivalentMinimalPureCubicalSubcomplex;

###############################
LM[1]:=ContractedComplex(LM[1]);
for i in [2..Length(LM)] do
LM[i]:=F(LM[i],LM[i-1]);
od;
###############################

##Initialize##############
M:=LM[1];
if true then #deg=0 then 
S:=ContractibleSubcomplexOfPureCubicalComplex(M); #SHOULD DELETE THIS!
else
S:=AcyclicSubcomplexOfPureCubicalComplex(M);
fi;
L:=[[M,S]];
triv:=List([1..Dimension(M)+1],i->0); triv[1]:=1;
bool:=true;
vsmaps:=[];
#########################

for i in [2..Length(LM)] do
  TM:=LM[i];
  TS:=HomotopyEquivalentMaximalPureCubicalSubcomplex(TM,L[Length(L)][2]);
  Add(L,[TM,TS]); 
od;


L[1][1]:=F(L[1][1],L[1][2]);

#################
for x in [1..Length(L)-1] do
  L[x+1][1]:=F(L[x+1][1],PureCubicalComplexUnion(L[x][1],L[x+1][2]));
od;
#################


L[1]:=ChainComplexOfPair(L[1][1],L[1][2]);

#################
for x in [1..Length(L)-1] do
  L[x+1]:=ChainComplexOfPair(L[x+1][1],L[x+1][2]);
  vsmaps[x]:= 
  TensorWithIntegersModP( ChainMapOfCubicalPairs(L[x],L[x+1]) , prime);
  vsmaps[x]:=HomologyVectorSpace(vsmaps[x],deg);
  Unbind(L[x]);
od;
#################

if deg>0 then
  return LinearHomomorphismsPersistenceMat(vsmaps);
fi;

L:=LinearHomomorphismsPersistenceMat(vsmaps);
for i in [1..Length(L)] do
for y in [i..Length(L)] do
  L[i][y]:=L[i][y]+1;
od;od;
return L;

end);
###############################################################
###############################################################

###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfPureCubicalComplex_Alt,
function(arg)
local
        M,Bettis,triv,Persist,deg,prime,TM,L,SL,bool,CollapseMat,
        Collapses_dim2,Fun,x,y;

##Input###################
deg:=arg[2];
prime:=arg[3];
if IsList(arg[1]) then L:=arg[1];
M:=L[1];
else 
M:=arg[1];
L:=[M];
fi;
deg:=Minimum(deg,Dimension(M));
##########################

if not deg in [0,Dimension(M)-1,Dimension(M)] then 
return fail; 
fi;

if IsList(arg[1]) then
	Bettis:=List(L,x->[Bettinumbers(x,0),Bettinumbers(x,deg)]);
else
	triv:=List([1..Dimension(M)+1],i->0);
	triv[1]:=1;
	Bettis:=[Bettinumbers(ContractedComplex(M))];

	bool:=true;
	while bool do
	TM:=ThickenedPureCubicalComplex(L[Length(L)]);
	Add(Bettis,Bettinumbers(ContractedComplex(TM)));
	bool:= not Bettis[Length(Bettis)]=triv;;
	Add(L,TM); 
	od;
	Bettis:=List(Bettis,b->[b[1],b[deg+1]]);
fi;


Persist:=List([1..Length(L)],x->List([1..Length(L)],y->0));;

###############################
if deg=0 then
for x in [1..Length(L)] do
for y in [x..Length(L)] do
Persist[x][y]:=Bettis[y][1];   
od;
od;
return Persist;
fi;
###############################

###############################
if deg=Dimension(M) then
return Persist;
fi;
###############################

###############################
Collapses_dim2:=function(X,Y)
##Number of cycles killed going from position X to position Y
local
	P, Acomp,A,Bcomp,TP, rows, cols, cnt, bettizero, Pathcomps, k, i, j;

Acomp:=PureCubicalComplex(FrameArray(L[Y]!.binaryArray));
A:=Acomp!.binaryArray;
Bcomp:=PureCubicalComplex(FrameArray(L[X]!.binaryArray));
rows:=Length(A);
cols:=Length(A[1]);
cnt:=0;
bettizero:=Bettis[X][1];

for k in [1..bettizero] do
P:=PathComponentOfPureCubicalComplex(Bcomp,k);

TP:=HomotopyEquivalentMaximalPureCubicalSubcomplex(Acomp,P);
TP:=TP!.binaryArray;

for i in [2..rows-1] do
for j in [2..cols-1] do
if A[i][j]=1 and TP[i][j]=0 then
if 4=Sum([TP[i-1][j]+TP[i+1][j]+TP[i][j-1]+TP[i][j+1]]) then cnt:=cnt+1;fi;;
fi;
od;
od;

od;

return cnt;
end;
###############################

#SL:=[ContractedComplex(L[1])];
#for x in [2..Length(L)] do
#Add(SL,HomotopyEquivalentMinimalPureCubicalSubcomplex(L[x],SL[x-1]));
#SL[x-1]:=HomotopyEquivalentMaximalPureCubicalSubcomplex(SL[x],SL[x-1]);
#od;
#L:=SL;


###############################
if deg=Dimension(M)-1 then

###############
Fun:=function(x,y);
#Ths function implements some partial recurrence in CollapseMat.
if CollapseMat[x+1][y]=0 then return CollapseMat[x][y-1]; fi;
if CollapseMat[x]=Bettis[x][2] then return CollapseMat[x][y-1]; fi;
if not y=Length(Bettis) then
   if Bettis[y+1]=0 then return Bettis[x][2]-Sum(CollapseMat[x]); fi;
fi;
return Collapses_dim2(x,y);
end;
###############

CollapseMat:=[];
for x in Reversed([1..Length(L)]) do
CollapseMat[x]:=List([1..Length(L)],a->0);
#CollapseMat[x]:=[];
if not x=Length(L) then 
CollapseMat[x][x+1]:=Collapses_dim2(x,x+1);
fi;
	for y in [x+2..Length(L)] do
	CollapseMat[x][y]:=Fun(x,y);;
	od;
od;

for x in [1..Length(L)] do
for y in [x..Length(L)] do
Persist[x][y]:=StructuralCopy(Bettis[x][2]);
if not x=y then
Persist[x][y]:=Persist[x][y]-CollapseMat[x][y];
fi;
od;
od;
return Persist;
fi;
###############################
end);
###############################################################
###############################################################


#########################################################
#########################################################
InstallGlobalFunction(BarCode,
function(P)
local PT,B,newrow,cols,rows,i,j,k,m;

PT:=MutableCopyMat(TransposedMat(P));
Add(PT,PT[1]*0);
PT:=TransposedMat(PT);
PT:=Concatenation([PT[1]*0],PT);
for i in Reversed([2..Length(PT)]) do
PT[i]:=PT[i]-PT[i-1];
od;

B:=[];
rows:=Length(P);
cols:=Length(P[1]);

for i in [1..rows] do
for j in Reversed([i..cols]) do
for k in [1..PT[i+1][j]-PT[i+1][j+1]] do
newrow:=List([1..cols],a->0);
for m in [i..j] do
newrow[m]:=1;
od;
Add(B,newrow);
od;

od;
od;

return B;
end);
###########################################################
###########################################################

###########################################################
###########################################################
InstallGlobalFunction(BarCodeDisplay,
function(arg)
local P,B,i,j,Disp,tmpDir,barcodedot,barcodegif;

###############
P:=arg[1];
if Length(arg)=2 then
Disp:=arg[2];
else
Disp:=DISPLAY_PATH;
fi;
###############

B:=BarCode(P);

if Length(arg)>2 then
B:=Filtered(B,x->Sum(x)>arg[3]);
fi;

tmpDir:=DirectoryTemporary();
barcodedot:=Filename(tmpDir,"tmpIn.log");
barcodegif:=Filename(tmpDir,"basic.gif");


PrintTo(barcodedot,"digraph finite_state_machine {\n\n");
AppendTo(barcodedot, "rankdir=LR;\n\n");
AppendTo(barcodedot,
"node [style=filled,shape=point]\n\n");

for i in [1..Length(B)] do
for j in [1..Length(B[1])] do
if B[i][j]=1 then
AppendTo(barcodedot,"node [color=black,fontcolor=black];\n",i,".",j,"\n");
else
AppendTo(barcodedot,"node [color=white,fontcolor=white];\n",i,".",j,"\n");
fi;
od;
od;

for i in [1..Length(B)] do
for j in [1..Length(B[1])-1] do
if B[i][j]=1 and B[i][j+1]=1 then
AppendTo(barcodedot,i,".",j,"->",i,".",j+1," [color=black,arrowhead=none];\n");
else
AppendTo(barcodedot,i,".",j,"->",i,".",j+1," [color=white,arrowhead=none];\n");
fi;
od;
od;

AppendTo(barcodedot,"}");

Exec(Concatenation("dot -Tgif ",barcodedot ,">", barcodegif));
Exec(Concatenation(Disp," ",barcodegif));
Sleep(2);
Exec(Concatenation("rm -r ",barcodegif{[1..Length(barcodegif)-9]}));

end);
#################################################################
#################################################################

###########################################################
###########################################################
InstallGlobalFunction(BarCodeCompactDisplay,
function(arg)
local P,B,i,j,Disp,tmpDir,barcodedot,barcodegif, lb;

###############
P:=arg[1];
if Length(arg)=2 then
Disp:=arg[2];
else
Disp:=DISPLAY_PATH;
fi;
###############

B:=BarCode(P);
B:=Collected(B);

tmpDir:=DirectoryTemporary();
barcodedot:=Filename(tmpDir,"tmpIn.log");
barcodegif:=Filename(tmpDir,"basic.gif");


PrintTo(barcodedot,"digraph finite_state_machine {\n\n");
AppendTo(barcodedot, "rankdir=LR;\n\n");
AppendTo(barcodedot,
"node [style=filled,shape=point]\n\n");

lb:=Length(B[1][1]);
for i in [1..Length(B)] do
for j in [1..Length(B[1][1])] do
if B[i][1][j]>0 then
AppendTo(barcodedot,"node [color=black,fontcolor=black];\n",i,".",j,"\n");
else
AppendTo(barcodedot,"node [color=white,fontcolor=white];\n",i,".",j,"\n");
fi;
od;
od;


for i in [1..Length(B)] do
for j in [1..Length(B[1][1])-1] do
if B[i][1][j]>0 and B[i][1][j+1]>0 then
  if j=Position(B[i][1],1) or (Position(B[i][1],1)=lb and j=lb-1) then
  AppendTo(barcodedot,i,".",j,"->",i,".",j+1," [label =\"",B[i][2]," \",color=black,arrowhead=none];\n");
  else
  AppendTo(barcodedot,i,".",j,"->",i,".",j+1," [color=black,arrowhead=none];\n");
  fi;
else
  if j=Position(B[i][1],1) or (Position(B[i][1],1)=lb and j=lb-1) then
  AppendTo(barcodedot,i,".",j,"->",i,".",j+1," [label =\"",B[i][2]," \",color=white,arrowhead=none];\n");
  else
  AppendTo(barcodedot,i,".",j,"->",i,".",j+1," [color=white,arrowhead=none];\n");
  fi;
fi;
od;
od;

AppendTo(barcodedot,"}");

Exec(Concatenation("dot -Tgif ",barcodedot ,">", barcodegif));
Exec(Concatenation(Disp," ",barcodegif));
Sleep(2);
Exec(Concatenation("rm -r ",barcodegif{[1..Length(barcodegif)-9]}));

end);
#################################################################
#################################################################




############################################################
############################################################
InstallGlobalFunction(PersistentHomologyOfFilteredChainComplex,
function(C,n,prime)
local
    SubComplex,
    ComplexInclusion,
    bars;

####################################################
SubComplex:=function(C,r)
local D,DIM,BND;

DIM:=function(n); if n=0 then return C!.dimension(0);
else return C!.filteredDimension(r,n); fi; end;

###############
if EvaluateProperty(C,"initial_inclusion")=false then
BND:=function(n,v)
local vv,uu;
vv:=C!.dimension(n)-C!.filteredDimension(r,n)+v;
if n>1 then
uu:=C!.dimension(n-1)-C!.filteredDimension(r,n-1)+1;
else
uu:=C!.dimension(0);;
fi;
return C!.boundary(n,vv){[uu..C!.dimension(n-1)]};
end;
fi;
##############

##############
if EvaluateProperty(C,"initial_inclusion")=true then
BND:=function(n,v)
local uu;
if n>1 then
uu:=C!.filteredDimension(r,n-1);
else
uu:=C!.dimension(0);
fi;
return C!.boundary(n,v){[1..uu]};
end;
fi;
##############

D:= Objectify(HapChainComplex,
                rec(
                dimension:=DIM ,
                boundary:=BND ,
                properties:=
                [["length",EvaluateProperty(C,"length")],
                ["connected",true],
                ["type", "chainComplex"],
                ["initial_inclusion",EvaluateProperty(C,"initial_inclusion")],
                ["characteristic",
                EvaluateProperty(C,"characteristic")] ]));

return D;
end;
####################################################


####################################################
ComplexInclusion:=function(C,k)
local S,T, M,map,prime;

S:=SubComplex(C,k);
T:=SubComplex(C,k+1);
prime:=EvaluateProperty(C,"characteristic");

##################
if EvaluateProperty(S,"initial_inclusion")=false then
map:=function(v,n)
local w,i,a,b;
w:=[1..T!.dimension(n)]*0;
a:=Length(w)-Length(v);
b:=Length(w);
for i in [a+1..b] do
w[i]:= v[i-a]*1;
od;
return w*One(GF(prime));
end;
fi;
##################

##################
if EvaluateProperty(S,"initial_inclusion")=true then
map:=function(v,n)
local w,i,a,b;
w:=[1..T!.dimension(n)]*0;
a:=Length(v);
for i in [1..a] do
w[i]:= v[i]*1;
od;
return w*One(GF(prime));
end;
fi;
##################



return Objectify(HapChainMap,
        rec(
           source:=S,
           target:=T,
           mapping:=map,
           properties:=[ ["type","chainMap"],
           ["characteristic",prime ]
           ]));

end;
####################################################


####################################################
bars:=function(C,n)
local FL,K,i;

FL:=EvaluateProperty(C,"filtration_length");

K:=[];
for i in [0..FL-1] do
Add(K,ComplexInclusion(C,i));
od;

Apply(K,x->HomologyVectorSpace(x,n));
return LinearHomomorphismsPersistenceMat(K);
end;
####################################################

return bars(TensorWithIntegersModP(C,prime),n);
end);
#############################################################
#############################################################

##########################################
##########################################
InstallGlobalFunction(TruncatedGComplex,
function(arg)
local R,a,b,Dimension, Boundary, Stabilizer, Action;

R:=arg[1];
a:=arg[2];
b:=arg[3];

##################
Dimension:=function(n);
if not n+a in [a..b] then return 0; fi;
return R!.dimension(n+a);
end;
##################

##################
Boundary:=function(n,i);
return R!.boundary(n+a,i);
end;
##################

##################
Stabilizer:=function(n,i);
return R!.stabilizer(n+a,i);
end;
##################

##################
Action:=function(n,k,g);
return R!.action(n+a,k,g);
end;
##################



if not IsHapNonFreeResolution(R) then return fail; fi;


return          Objectify( HapNonFreeResolution,
                rec(
                dimension:=Dimension,
                boundary:=Boundary,
                homotopy:=fail,
                elts:=R!.elts,
                group:=R!.group,
                stabilizer:=Stabilizer,
                action:=Action,
                properties:=
                   [["length",EvaluateProperty(R,"length")],
                    ["reduced",true],
                    ["type","resolution"],
                    ["characteristic",EvaluateProperty(R,"characteristic")]  ])
                );

end);
##########################################
##########################################


###############################################################
###############################################################
InstallGlobalFunction(ZZPersistentHomologyOfPureCubicalComplex,
function(arg)
local
        LM,deg,prime,
        ChnCmps,vsmaps,L,
        U,CU,V,CV,W,CW,
        x,y,i,F,PC,correction;

##Input###################
LM:=arg[1];
deg:=arg[2];
prime:=arg[3];
##########################

correction:=[];
F:=HomotopyEquivalentMinimalPureCubicalSubcomplex;

##Initialize##############
U:=LM[1];
CU:=ContractibleSubcomplexOfPureCubicalComplex(U);
L:=[[U,CU]];
vsmaps:=[];
#########################


#########################
for i in [1..Length(LM)-1] do
# U --> V <-- W

U:=L[Length(L)][1];
CU:=L[Length(L)][2];

W:=LM[i+1];

V:=PureCubicalComplexUnion(U,W);
if Size(CU)>0 then
CV:=HomotopyEquivalentMaximalPureCubicalSubcomplex(V,CU);
else
CV:=ContractibleSubcomplexOfPureCubicalComplex(V);
fi;

CW:=PureCubicalComplexIntersection(CV,W);
CW:=ContractibleSubcomplexOfPureCubicalComplex(CW);

if Size(CW)=0 then Add(correction,i+1); fi;

Add(L,[V,CV]);
Add(L,[W,CW]);
od;
#####################

L[1]:=ChainComplexOfPair(L[1][1],L[1][2]);

#################
for i in [1..Length(LM)-1] do
  L[2*i]:=ChainComplexOfPair(L[2*i][1],L[2*i][2]);
  vsmaps[2*i-1]:=
  TensorWithIntegersModP( ChainMapOfCubicalPairs(L[2*i-1],L[2*i]) , prime);
  vsmaps[2*i-1]:=HomologyVectorSpace(vsmaps[2*i-1],deg);

L[2*i+1]:=ChainComplexOfPair(L[2*i+1][1],L[2*i+1][2]);
  vsmaps[2*i]:=
  TensorWithIntegersModP( ChainMapOfCubicalPairs(L[2*i+1],L[2*i]) , prime);
  vsmaps[2*i]:=HomologyVectorSpace(vsmaps[2*i],deg);
Unbind(L[2*i]);
  Unbind(L[2*i-1]);

od;
#################


if deg>0 then
  return LinearHomomorphismsZZPersistenceMat(vsmaps);
fi;


L:=LinearHomomorphismsZZPersistenceMat(vsmaps);
for i in [1..Length(L)] do
for y in [i..Length(L)] do
if Length(Intersection([i..y],correction))=0 then
  L[i][y]:=L[i][y]+1;
fi;
od;od;
return L;

end);
###############################################################
###############################################################


###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfFilteredPureCubicalComplex,
function(F,n)
local FF, C, D;

############################
if n=0 and IsHapFilteredPureCubicalComplex(F) then
D:=Dendrogram(F);
return DendrogramToPersistenceMat(D);;
fi;
############################

if IsHapFilteredPureCubicalComplex(F) then
#FF:=ZigZagContractedFilteredPureCubicalComplex(F,2);
FF:=ContractedFilteredPureCubicalComplex(F);
FF:=FilteredPureCubicalComplexToCubicalComplex(FF);
fi;

if IsHapFilteredPureCubicalComplex(F) then
FF:=F;
fi;

C:=SparseFilteredChainComplexOfFilteredCubicalComplex(FF);
if Dimension(F)<=3 then
C:=TensorWithIntegersModPSparse(C,2);
fi;

return
PersistentHomologyOfFilteredSparseChainComplex(C,n);
end);
###############################################################
###############################################################

###############################################################
###############################################################
InstallGlobalFunction(PersistentHomologyOfFilteredPureCubicalComplex_alt,
function(FF,deg)
local flen, HEMin, HEMax, FT, PH, D, PHfirst, Pmat, PB, F, r, s,n;

#############################
if IsBound(FF!.persistentHomology) then
if IsBound(FF!.persistentHomology[deg+1]) then
return FF!.persistentHomology[deg+1];
fi;
fi;
############################


############################
if deg=0 then 
D:=Dendrogram(FF);
FF!.persistentHomology:= [DendrogramToPersistenceMat(D)];;
return FF!.persistentHomology[1];
fi;
############################


F:=ZigZagContractedFilteredPureCubicalComplex(FF);


PB:=Dendrogram(F);
PB:=[DendrogramToPersistenceMat(PB)];
for r in [1..deg] do
Add(PB,PB[1]*0);
od;

flen:=F!.filtrationLength;

HEMin:=HomotopyEquivalentMinimalPureCubicalSubcomplex;
HEMax:=HomotopyEquivalentMaximalPureCubicalSubcomplex;
FT:=[];


###############################################################
###############################################################
PHfirst:= function(F,deg,r)
local
        Fr,Fs, S, C, s, n;

Fr:=FiltrationTerm(F,r);
ContractPureCubicalComplex(Fr);   

S:=AcyclicSubcomplexOfPureCubicalComplex(Fr);
FT[r]:=[Fr];

C:=SparseChainComplexOfPair(Fr,S);

for n in [1..deg] do
#if not IsBound(PB[n+1][r][r]) then
Pmat[r][r][n]:=Bettinumbers(C,n);
PB[n+1][r][r]:=Pmat[r][r][n];
#fi;
od;

end;
###############################################################
###############################################################


###############################################################
###############################################################
PH:= function(F,deg,r)
local
        Fr,Fs, S, C, s, n;

Fr:=FT[r][1];

for s in [r+1..flen] do

Fs:=FiltrationTerm(F,s);
Fs:=HEMin(Fs,Fr);
S:=HEMax(Fs,Fr);

C:=SparseChainComplexOfPair(Fs,S);

for n in [1..deg] do
#if not IsBound(PB[n+1][r][s]) then 
Pmat[r][s][n]:=Bettinumbers(C,n);
PB[n+1][r][s]:=Pmat[s][s][n]-Pmat[r][s][n]+PB[n][r][r]-PB[n][r][s];
#fi;
od;

if Sum(List([1..deg],n->PB[n+1][r][s]))=0 then break; fi;
od;

end;
###############################################################
###############################################################

Pmat:=[];
for r in [1..flen] do
Pmat[r]:=[];
for s in [1..flen] do
Pmat[r][s]:=[];
od;
od;

for r in [1..flen] do
PHfirst(F,deg,r);
od;

for r in [1..flen] do
PH(F,deg,r);
od;

FF!.persistentHomology:=PB;

return PB[deg+1];
end);
###############################################################
###############################################################