
//========================================================================================
// Functions for computing covers of finite graphs.
//
// These functions are to be used within Magma.
//
// Copyright (C) Rok Pozar 2014.
//========================================================================================


/* to be able to work with a generalised undirected graph X (allowing multiedges, cycles, semiedges) we ``exploit'' a Magma multigraph object GrphMultUnd: 
 * first note that the undirected edge {u, v} is stored twice in the list, and so there are two positions in the list associated with the edge; 
 * these two positions (or indices) are contiguous and could be seen as darts of the corresponding edge, such that the initial vertex of the odd
 * index (or dart) is Min{Index(u),Index(v)}, while the terminal vertex is Max{Index(u),Index(v)}, and vice versa for the even dart;
 * therefore semiedges are represented as follows:
 * each semiedge at vertex v is given as a cycle at v with additional information that enables to distinguish
 * between cycles and semiedges i.e. this information tells when two contiguous indices (or darts) are "identified";
 * it is enough to encode this information with a map f : S -> [true, false] on a sequence S of odd indices(darts) such that f(i) = true iff 
 * i and i + 1 are identified darts iff the edge e = EdgesSet(X).i is semiedge;
 * further, voltages are stored just for odd indices (or darts) too;
 * NOTICE: all base graphs are assumed to be connected;
 */

//========================================================================================
//                       Construction of a voltage space 
//========================================================================================

VoltageSpace := function(t, G, X, S, semi, Volt)
/* Given a generalised undirected graph (i.e. a triple 'X'::GrphMultUnd, 'S'::SetIndx[RngIntElt], 'semi'::SeqEnum[BoolElt], where S is an indexed set of odd indices (darts) of X  and 
 * semi is a sequence such that S[i] is semiedge iff semi[i] eq true), a voltage group 'G'::GrpFin and a sequence  'Volt'::SeqEnum[GrpElt] of group elements (voltages) such that Volt[i]
 * corresponds to S[i] this function constructs a voltage space on X of a type 't'::MonStgElt.
 * String t can be "Cayley", "Permutation" or "CayleyToPermutation" (tells the type of a voltage assignment) and consequently indirectly
 * defines the set on which a finite group 'G'::GrpFin (a voltage group) acts:
 * type "Cayley": group G acts on G by right translation;
 * type "Permutation": permutation group G acts on a set {1,2,..,n}, where n is the cardinality of the fibre;
 * type "CayleyToPermutation": permutation group G acts regularly on set {1,2,..,n}, where n is the cardinality of the fibre.
 * The function returns a "record''::Rec with the following fields:
 * type::MonStgElt - a string describing the type of the voltage assignment,
 * group::Grp - a voltage group G,
 * graph::GrphMultUnd - a multigraph X,
 * assignment::Map[SetIndx[RngIntElt],Grp] - a map from the sequence S to the group G,
 * semiedge::Map[SetIndx[RngIntElt],SetIndx[BoolElt]] - a map from the sequence S to sequence the semi.
 */

	local  IsSemiedge, voltassign, CRF, record;
					
	if Type(t) ne MonStgElt then
		error "Bad argument types\n Argument types given: ",
		ExtendedType(t),",",ExtendedType(G),",",ExtendedType(X),",",ExtendedType(S),",",ExtendedType(semi),",",ExtendedType(Volt);
	else
		if (t ne "Cayley") and (t ne "Permutation") and (t ne "CayleyToPermutation") then
			error "Argument 1 (" cat t cat ") should be in the range {Cayley, Permutation, CayleyToPermutation}";
		end if;
	end if;
	if not ISA(Type(G),Grp) then
		error "Bad argument types\n Argument types given: ",
		ExtendedType(t),",",ExtendedType(G),",",ExtendedType(X),",",ExtendedType(S),",",ExtendedType(semi),",",ExtendedType(Volt);
	else
		if Order(G) eq 0 then 
			error "Argument 2 should be finite";
		end if;
	end if;
	if Type(X) ne GrphMultUnd then
		error "Bad argument types\n Argument types given: ",
		ExtendedType(t),",",ExtendedType(G),",",ExtendedType(X),",",ExtendedType(S),",",ExtendedType(semi),",",ExtendedType(Volt);
	end if;
	if ExtendedType(S) ne SetIndx[RngIntElt]then
		error "Bad argument types\n Argument types given: ",
		ExtendedType(t),",",ExtendedType(G),",",ExtendedType(X),",",ExtendedType(S),",",ExtendedType(semi),",",ExtendedType(Volt);
	else
		for i in S do
			if not i in SequenceToSet(Indices(X)) then
				error "Argument 4 (" cat IntegerToString(i) cat") should be in the range of odd indices of argument 3";
			end if;
		end for;
						
		if #Indices(X) ne #S then
			error "Argument 4 should be of length", #Indices(X);
		end if; 
	end if;
	if ExtendedType(semi) ne SeqEnum[BoolElt] then
		error "Bad argument types\n Argument types given: ",
		ExtendedType(t),",",ExtendedType(G),",",ExtendedType(X),",",ExtendedType(S),",",ExtendedType(semi),",",ExtendedType(Volt);
	else 
		if #Indices(X) ne #semi then
			error "Argument 5 should be of length", #Indices(X);
		end if;		
	end if;
	if not ISA(ExtendedType(Volt),SeqEnum[GrpElt]) then
		error "Bad argument types\n Argument types given: ",
		ExtendedType(t),",",ExtendedType(G),",",ExtendedType(X),",",ExtendedType(S),",",ExtendedType(semi),",",ExtendedType(Volt);
	else 
		if 	not IsCoercible(G,Volt[1]) then
			error "Arguments 2 and 6 have no covering structure";
		end if;
		if #Indices(X) ne #Volt then
			error "Argument 6 should be of length", #Indices(X);
		end if;
		for i in [1..#semi] do 
			if semi[i] eq true then
				if TerminalVertex(EdgeSet(X).S[i]) ne InitialVertex(EdgeSet(X).S[i]) then
					error "Edge (" cat IntegerToString(S[i]) cat") can not be semiedge";
				else
					if Type(G) eq GrpAb then
						if (Volt[i] + Volt[i]) ne Id(G) then
							error "Voltage of semiedge (" cat IntegerToString(S[i]) cat") should be of order 2";
						end if;
					else
						if (Volt[i]*Volt[i]) ne Id(G) then
							error "Voltage of semiedge (" cat IntegerToString(S[i]) cat") should be of order 2";
						end if;
					end if;
				end if;
			end if;
		
		end for;
	end if;
		 
	voltassign := map<S -> G|{<S[i], Volt[i]> : i in [1..#S]}>;

	IsSemiedge := map<S -> {@true,false@}|{<S[i], semi[i]> : i in [1..#S]}>;

	CRF := recformat< type, group, graph, assignment, semiedge>;
	record := rec<CRF | type:= t, group := G, graph := X, assignment := voltassign, semiedge := IsSemiedge>;

	return record;

end function;

//========================================================================================
//                         Computing the local voltage group 
//========================================================================================

LocalVoltageGroup := function(VS, u)
// Given a voltage space 'VS'::Rec and a vertex 'u'::GrphVert this function returns the local voltage group P::Grp at u.

	local X, voltageassignment, G, isAb, colour, Q, potential, voltage, dart, multi, currentVertex, nextVertex, P, V, num;
	X := VS`graph;
		
	P := {};                                             // generating generators of the local voltage group
	voltageassignment := VS`assignment;
	G := VS`group;
	if Type(G) eq GrpAb then
		isAb := true;
	else
		isAb := false;
	end if;

	V := VertexSet(X);
	n := #V;
	num := 0;                                          // number of vertices that  have already been searched
	colour := [0 : i in [1..n]];                    // all vertices are coloured white(0)
	potential := [Id(G) : i in [1..n]];          // potential[i] is the voltage of the unique path on the tree from the vertex V.i to the vertex u
	Q := [Index(u)];                                 // vertices that has to be searched
			
	colour[Index(u)] := 1;                         // vertex u is coloured grey(1), since it has been searched at the moment
	
	while IsEmpty(Q) eq false do         // BFS 
		currentVertex := V.(Q[1]);
		Remove(~Q,1);
		for nextVertex in Neighbours(currentVertex) do                                                     // search all neighbours of the current vertex
			if colour[Index(nextVertex)] eq 1 then                                                              // a cotree edge
				multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
				for dart in multi do;                                                      // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						if isAb then
							voltage := -(potential[Index(nextVertex)]) + voltageassignment(dart) + potential[Index(currentVertex)];
						else
							voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
						end if;
					else
						if isAb then
							voltage := -(potential[Index(currentVertex)]) + voltageassignment(dart) + potential[Index(nextVertex)];
						else
							voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment(dart)*potential[Index(nextVertex)];
						end if;
					end if;  
					if G!voltage notin sub<G|P> then
						Include(~P,voltage);
					end if;
				end for;	
			elif colour[Index(nextVertex)] eq 0 then                                        // a tree edge
				multi := Indices(currentVertex, nextVertex);                                // all multiple edges
				dart := multi[1];                                                           // edge with the first index is a tree edge
				colour[Index(nextVertex)] := 1; 
				Append(~Q,Index(nextVertex));
				Remove(~multi,1);
					
				// computing voltage
				if Index(nextVertex) lt Index(currentVertex) then
					voltage := voltageassignment(dart);
				else
					if isAb then
						voltage := -voltageassignment(dart);
					else
						voltage := voltageassignment(dart)^(-1);
					end if;
				end if;  
				
				if isAb then
					potential[Index(nextVertex)] := voltage + potential[Index(currentVertex)];
				else
					potential[Index(nextVertex)] := voltage*potential[Index(currentVertex)];
				end if;
				for dart in multi do;                                                         // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						if isAb then
							voltage := -(potential[Index(nextVertex)]) + voltageassignment(dart) + potential[Index(currentVertex)];
						else
							voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
						end if;
					else
						if isAb then
							voltage := -(potential[Index(currentVertex)]) + voltageassignment(dart) + potential[Index(nextVertex)];
						else
							voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment(dart)*potential[Index(nextVertex)];
						end if;
					end if;  
					if G!voltage notin sub<G|P> then
						Include(~P,voltage);
					end if;
				end for;
			end if;
		end for;
		colour[Index(currentVertex)] := 2;                                                    //  searching  is completed
		num := num + 1;
	end while;

	if num lt n then
		error " Graph in first argument is not connected";
	else
		return sub<G|P>;
	end if;
	
end function;


//========================================================================================
//               Testing whether the local voltage group acts transitively
//========================================================================================

IsLocallyTransitive := function(VS)
/* Given a voltage space 'VS'::Rec this function tests whether it is locally transitive i.e. 
 * whether the local voltage group at some vertex u (and hence at all) acts transitively.
 */

	local X, u, loc;
	X := VS`graph;
	
	u := Vertices(X)[1];
	loc := LocalVoltageGroup(VS,u);
	
	if VS`type eq "Cayley" then
		return (VS`group eq loc);
	else
		return IsTransitive(loc);
	end if;
	
end function;

//========================================================================================
//                   Construction of a (T,u)-reduced voltage assignment 
//========================================================================================

TReduced := function(VS,T,u)
/* given a voltage space 'VS'::Rec, sequence 'T'::SeqEnum[GrphEdge] of spanning tree edges and vertex 'u'::GrphVert, this function constructs a (T,u)-reduced 
 * voltage assignment
 * we assume that T is a spanning tree
 */

	local X, neighbours, initial, terminal, indices, colour, index, G, isAb, voltageassignment, n, potential, Q, tree, currentVertex, nextVertex, dart, voltage, D, newvoltages, newassignment;
	
	X := VS`graph;
	n := #VertexSet(X);
	neighbours := [[ ] : i in [1..n]];
	indices := [[ ] : i in [1..n]];
	G := VS`group;
	if Type(G) eq GrpAb then
		isAb := true;
	else
		isAb := false;
	end if;
	
	colour := [0 : i in [1..n]];             // all vertices are coloured white(0)
	voltageassignment := VS`assignment;
	potential := [Id(G) : i in [1..n]];      // potential[i] is the voltage of the unique path on the tree from vertex i to vertex u
	Q := [Index(u)];                         // vertices that has to be searched
	tree := {Index(i) : i in T};
	
	for e in T do
		initial := Index(InitialVertex(e));
		terminal := Index(TerminalVertex(e));
		index := Index(e);
		neighbours[initial] := Append(neighbours[initial],terminal);
		neighbours[terminal] := Append(neighbours[terminal],initial);
		indices[initial] := Append(indices[initial],index);
		indices[terminal] := Append(indices[terminal],index);
	end for;
	
	while IsEmpty(Q) eq false do                                                        // BFS of tree T
		currentVertex := Q[1];
		Remove(~Q,1);
		for i in [1..#neighbours[currentVertex]] do                                     // search all neighbours of current vertex
			nextVertex := neighbours[currentVertex][i];
			if colour[nextVertex] eq 0 then
				colour[nextVertex] := 1;
				Append(~Q,nextVertex);
				dart := indices[currentVertex][i];
				
				// computing voltage
				if nextVertex lt currentVertex then
					voltage := voltageassignment(dart);
				else
					if isAb then
						voltage := -voltageassignment(dart);
					else
						voltage := voltageassignment(dart)^(-1);
					end if;
				end if;
				if isAb then
					potential[nextVertex] := voltage + potential[currentVertex];
				else
					potential[nextVertex] := voltage*potential[currentVertex];
				end if;
			end if;
		end for;
	
	end while;

	D := Domain(VS`assignment);	
	newvoltages := [];
	for dart in D do
		if dart in tree then
			Append(~newvoltages,Id(G));
		else
			initial := Index(InitialVertex(EdgeSet(X).dart));
			terminal := Index(TerminalVertex(EdgeSet(X).dart));
			if isAb then
				voltage := -(potential[initial]) + voltageassignment(dart) + potential[terminal];
			else 
				voltage := (potential[initial])^(-1)*voltageassignment(dart)*potential[terminal];
			end if;
			Append(~newvoltages,voltage);
		end if;
		
	end for;
		
	newassignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
	return newassignment;
	
end function;


Reduced := function(VS, u)

/* given a voltage space 'VS'::Rec and a vertex 'u'::GrphVert, this function constructs a BFS tree T rooted at u together with (T,u)-reduced 
 * voltage assignment, where T is given as a sequence of edges;
 */

	local X, voltageassignment, G, isAb, colour, Q, tree, potential, voltage, dart, multi, currentVertex, nextVertex, V, E, domain, newvoltages, newassignment, num;
	
	X := VS`graph;
	
	voltageassignment := VS`assignment;
	G := VS`group;
	if Type(G) eq GrpAb then
		isAb := true;
	else
		isAb := false;
	end if;
	V := VertexSet(X);
	E := EdgeSet(X);
	//e := #E;
	n := #V;
	num := 0;                            // number of vertices that have already been searched
	newvoltages := [];                   // new voltage at the position i corresponding  to dart at the position i in domain
	domain := [];	
	colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
	tree := [];                          // edges of tree T
	potential := [Id(G) : i in [1..n]];  // potential[i] is the voltage of the unique path on the tree from vertex V.i to vertex u
	Q := [Index(u)];                     // vertices that has to be searched
			
	colour[Index(u)] := 1;               // vertex u is coloured grey(1), since it has been searched at the moment
	
	while IsEmpty(Q) eq false do         // BFS 
		currentVertex := V.(Q[1]);
		Remove(~Q,1);
		for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
			if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
				multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
				for dart in multi do;                                                      // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						if isAb then
							voltage := -(potential[Index(nextVertex)]) + voltageassignment(dart) + potential[Index(currentVertex)];
						else
							voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
						end if;
					else
						if isAb then
							voltage := -(potential[Index(currentVertex)]) + voltageassignment(dart) + potential[Index(nextVertex)];	
						else 
							voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment(dart)*potential[Index(nextVertex)];	
						end if;
					end if; 
					Append(~domain,dart);
					Append(~newvoltages,voltage); 
				end for;	
			elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
				multi := Indices(currentVertex, nextVertex);                                // all multiple edges
				dart := multi[1];                                                           // edge with the first index is tree edge
				Append(~tree,E.dart);
				Append(~domain,dart);
				Append(~newvoltages,Id(G));
				colour[Index(nextVertex)] := 1; 
				Append(~Q,Index(nextVertex));
				Remove(~multi,1);
					
				// computing voltage
				if Index(nextVertex) lt Index(currentVertex) then
					voltage := voltageassignment(dart);
				else
					if isAb then
						voltage := -voltageassignment(dart);
					else
						voltage := voltageassignment(dart)^(-1);
					end if;
				end if;  
				
				if isAb then
					potential[Index(nextVertex)] := voltage + potential[Index(currentVertex)];
				else
					potential[Index(nextVertex)] := voltage*potential[Index(currentVertex)];
				end if;
				for dart in multi do;                                                         // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						if isAb then
							voltage := -(potential[Index(nextVertex)]) + voltageassignment(dart) + potential[Index(currentVertex)];
						else
							voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
						end if;
					else
						if isAb then
							voltage := -(potential[Index(currentVertex)]) + voltageassignment(dart) + potential[Index(nextVertex)];
						else 
							voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment(dart)*potential[Index(nextVertex)];
						end if;
					end if; 
					Append(~domain,dart);
					Append(~newvoltages,voltage); 
				end for;
			end if;
		end for;
		colour[Index(currentVertex)] := 2;                                                    //  searching of vertex is completed
		num := num + 1;
	end while;

	if num lt n then
		error "Graph in first argument is not connected";
	else
		newassignment := map<domain -> G|{<domain[i], newvoltages[i]> : i in [1..#domain]}>;
		return newassignment, tree;
	end if;
	
end function;

//========================================================================================
//                Converting the Cayley voltage space to the permutation voltage space 
//========================================================================================

CayleyToPermutation:= function(CVS)
/* function CayleyToPermutation converts a given Cayley voltage space 'CVS': Rec to voltage space of the type "CayleyToPermutation" via its regular
 * representation;
 */ 

	local f, P, F, D, assign, newassignment, newvoltages, im;

	F := CVS`group;
	f, P := CosetAction(F,sub<F|>);
	assign := CVS`assignment;
	D := Domain(CVS`assignment);

	newvoltages := [];
	for i in D do
		Append(~newvoltages,f(assign(i)));
	end for;

	newassignment := map<D -> P|{<D[i], newvoltages[i]> : i in [1..#D]}>;

	CRF := recformat< type, group, graph, assignment, semiedge>;
	record := rec<CRF | type:= "CayleyToPermutation", group := Sym(Degree(P)), graph := CVS`graph, assignment := newassignment, semiedge := CVS`semiedge>;

	return record;

end function;

//========================================================================================
//               Construction of the derived cover
//========================================================================================

DerivedCover := function(VS)
/* given the permutation voltage space 'VS'::Rec (if the Cayley voltage space is given, it is first converted to the permutation voltage space) the derived cover is
 * constructed together with the projective map which maps from the dart set of derived cover to the dart set of the base graph: if an edge e in the derived cover 
 * is mapped to a semiedge s, then both respective darts of e are mapped to even dart of s;
 */
 
	local IsSemiedge, X, G, voltageassignment, D, vertices, edges, n, f, F, semi, num, M, V, newsemiedge, projective, projectiveMap, I;
	
	if VS`type eq "Cayley" then
		VS := CayleyToPermutation(VS);
	end if;
	
	
	X := VS`graph;
	IsSemiedge := VS`semiedge;
	voltageassignment := VS`assignment;
	G := VS`group;
	vertices := {@ @};
	D := [i : i in [1..(#EdgeSet(X)*2)]];
	edges := [];
	semi := [];
	projective := [];

	n := Degree(G);
	for v in VertexSet(X) do                                      // construction of the vertex set of the derived graph 
		for i in [1..n] do
			Include(~vertices,<v,i>);
		end for;
	end for;
	
	M := MultiGraph<vertices|>;
	V := VertexSet(M);
	
	for e in EdgeSet(X) do                                                     // construction of the edge set of the derived graph with the projective map
		if TerminalVertex(e) eq InitialVertex(e) then
			if IsSemiedge(Index(e)) then                                       // semiedge in the base graph
				num := [1 : i in [1..n]];
				for i in [1..n] do
					if num[i] eq 1 then
						if i^voltageassignment(Index(e)) eq i then             // semiedge in the derived cover
							Append(~semi, true);
							Append(~edges, <<InitialVertex(e),i>,{<InitialVertex(e),i>}>);
							Append(~projective, Index(e));
							Append(~projective, Index(e));
						else                                                   // link in the derived cover within the fibre
							Append(~semi, false);
							Append(~edges, <<InitialVertex(e),i>,{<InitialVertex(e),i^voltageassignment(Index(e))>}>);
							Append(~projective, Index(e));
							Append(~projective, Index(e));
							num[i^voltageassignment(Index(e))] := 0;
						end if;
					end if;
				end for;	
			else                                                               // loop in the base graph
				for i in [1..n] do
					Append(~semi, false);
					Append(~edges, <<InitialVertex(e),i>,{<InitialVertex(e),i^voltageassignment(Index(e))>}>);
					if i^voltageassignment(Index(e)) eq i then                 // loop in the derived cover
						Append(~projective, Index(e));
						Append(~projective, Index(e) + 1);
					else                                                       // link in the derived cover within the fibre
						if Index(V!<InitialVertex(e),i>) lt Index(V!<InitialVertex(e),i^voltageassignment(Index(e))>) then
							Append(~projective, Index(e));
							Append(~projective, Index(e) + 1);
						else 
							Append(~projective, Index(e) + 1);
							Append(~projective, Index(e));
						end if;
					end if;
				end for;
			end if;
		else                                                                    // link in the base graph
			for i in [1..n] do
				Append(~semi, false);
				Append(~edges, <<InitialVertex(e),i>,{<TerminalVertex(e),i^voltageassignment(Index(e))>}>);
				if Index(V!<InitialVertex(e),i>) lt Index(V!<TerminalVertex(e),i^voltageassignment(Index(e))>) then
					Append(~projective, Index(e));
					Append(~projective, Index(e) + 1);
				else
					Append(~projective, Index(e) + 1);
					Append(~projective, Index(e));
				end if;
			end for;
		end if;
	end for;
	
	M := MultiGraph<vertices|edges>;
	I := [i : i in [1..2*#EdgeSet(M)]];
	projectiveMap := map<I -> D|{<I[i], projective[i]> : i in [1..2*#EdgeSet(M)]}>;
	newsemiedge := map<Indices(M) -> [true,false]|{<Indices(M)[i], semi[i]> : i in [1..#Indices(M)]}>;

	return  M, newsemiedge, projectiveMap;
	
end function;

//========================================================================================
//                     Testing whether an automorphism lifts 
//========================================================================================


AutOnVertex := function(X, f, v) 

//for an automorphism 'f'::GrpPermElt (given on darts) of the connected graph 'X'::GrphMultUnd this function returns the f-image
// of a vertex  'v'::GrphVert 

	local currentVertex, neighbourVertex , dart;
	
	if not IsConnected(X) then
			error " First argument is not connected";
	end if;
	
	currentVertex := VertexSet(X)!v;
	neighbourVertex := SetToSequence(Neighbours(currentVertex))[1];
	dart := Indices(neighbourVertex,currentVertex)[1];
	
	if currentVertex lt neighbourVertex then
		if dart^f mod 2 eq 0 then
			return TerminalVertex(EdgeSet(X).(dart^f));
		else
			return InitialVertex(EdgeSet(X).(dart^f));
		end if; 
	else                                                   // if currentVertex gt neighbourVertex or currentVertex eq neighbourVertex
		if (dart+1)^f mod 2 eq 0 then
			return TerminalVertex(EdgeSet(X).((dart+1)^f));
		else
			return InitialVertex(EdgeSet(X).((dart+1)^f));
		end if; 
	end if;

end function;


CovProjOnVertex := function(X,Y,p,v)

/* for a covering projection 'p'::Map (given on darts) mapping from covering graph 'X'::GrphMultUnd to base graph 'Y'::GrphMultUnd this function 
 * returns the projection of the vertex 'v'::GrphVert of X; 
 * notice: graphs X and Y can be generalized graphs, the information about the semiedges is encoded in projection p;
 */

	local currentVertex, neighbourVertex , dart;
	
	currentVertex := VertexSet(X)!v;
	neighbourVertex := SetToSequence(Neighbours(currentVertex))[1];
	dart := Indices(neighbourVertex,currentVertex)[1];
	

	if currentVertex lt neighbourVertex then
		if p(dart) mod 2 eq 0 then
			return TerminalVertex(EdgeSet(Y).(p(dart)));
		else
			return InitialVertex(EdgeSet(Y).(p(dart)));
		end if; 
	else                                                    // if currentVertex gt neighbourVertex or currentVertex eq neighbourVertex
		if p(dart) mod 2 eq 0 then
			return InitialVertex(EdgeSet(Y).(p(dart)));
		else
			return TerminalVertex(EdgeSet(Y).(p(dart)));
		end if; 
	end if;
	
end function;


SolveSystem := function(A, B, i) //A,B must define transitive permutation groups

	local n, r, S, TreeA, TreeB, permutation, conflict, cp, tp;
	
	n := Degree(Parent(A[1]));
	r := #A;
	S := Sym(n);
	
	TreeA := [0 : i in [1..n]];
	TreeB := [0 : i in [1..n]];
	permutation := [0 : i in [1..n]];
	TreeA[1] := 1;
	TreeB[1] := i;
	permutation[1] := i;
	conflict := 0;
	cp := 1;
	tp := 2;
	while (cp le n) and (conflict eq 0) do
		for j := 1 to r do
			if permutation[TreeA[cp]^A[j]] eq 0 then
				permutation[TreeA[cp]^A[j]] := TreeB[cp]^B[j];
				TreeA[tp] := TreeA[cp]^A[j];
				TreeB[tp] := TreeB[cp]^B[j];
				tp := tp + 1;
			else
				if permutation[TreeA[cp]^A[j]] ne TreeB[cp]^B[j] then
					conflict := 1;
					break;
				end if;
			end if;
		end for;
		cp := cp + 1;
	end while;
	
	if conflict eq 0 then
		return true;
	else 
		return false;
	end if;
	
end function;	
HasLift := function(VS, f)

/* given a voltage space 'VS'::Rec and an automorphism 'f'::GrpPermElt (given on darts) of the base graph this function tests whether f lifts via solving a system of permutation equations; 
 * in case f has a lift function returns true together with a sequence S of integers;
 * notice that each lift (if exists) is determined by an image of one dart;
 * consequently an integer S[i] represents a lift with a property that the vertex <VertexSet(X).1,1> is mapped to the vertex <f(VertexSet(X).1),S[i]>;
 */

	local X, G, voltageassignment, V, n, A, B, colour, potential, fpotential, Q, currentVertex, currentfVertex, nextVertex,
		  nextfVertex, voltage, fvoltage, num, solution, S;
	
	if VS`type eq "Cayley" then
		VS := CayleyToPermutation(VS);
	end if;
	
	X := VS`graph;
	voltageassignment := VS`assignment;
	G := VS`group;
	V := VertexSet(X);
	n := #V;
	
	// construction of a system of permutation equations
	A := [];
	B := [];
	colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
	potential := [Id(G) : i in [1..n]];  // potential[i] is the voltage of the unique path on the tree from vertex V.i to vertex u
	fpotential := [Id(G) : i in [1..n]]; // fpotential[i] is the voltage of the unique path on the tree from vertex f(V.i) to vertex f(u)
	Q := [1];                            // vertices that has to be searched
	colour[1] := 1;                      // vertex u is coloured grey(1), since it has been searched at the moment
	num := 0;                            // number of searched vertices
			
	while IsEmpty(Q) eq false do         // BFS 
		currentVertex := V.(Q[1]);
		currentfVertex := AutOnVertex(X,f,currentVertex);
		Remove(~Q,1);
		for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
			nextfVertex := AutOnVertex(X,f,nextVertex);
			if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
				multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
				for dart in multi do;                                                      // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
						if dart^f mod 2 eq 0 then
							fvoltage := (fpotential[Index(nextfVertex)])^(-1)*voltageassignment(dart^f-1)^(-1)*fpotential[Index(currentfVertex)];
						else
							fvoltage := (fpotential[Index(nextfVertex)])^(-1)*voltageassignment(dart^f)*fpotential[Index(currentfVertex)];
						end if;
					else
						voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment(dart)*potential[Index(nextVertex)];
						if dart^f mod 2 eq 0 then
							fvoltage := (fpotential[Index(currentfVertex)])^(-1)*voltageassignment(dart^f-1)^(-1)*fpotential[Index(nextfVertex)];
						else
							fvoltage := (fpotential[Index(currentfVertex)])^(-1)*voltageassignment(dart^f)*fpotential[Index(nextfVertex)];
						end if;
					end if; 
					Append(~A,voltage);
					Append(~B,fvoltage);
				end for;	
			elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
				multi := Indices(currentVertex, nextVertex);                                // all multiple edges
				dart := multi[1];                                                           // edge with the first index is tree edge
				colour[Index(nextVertex)] := 1; 
				Append(~Q,Index(nextVertex));
				Remove(~multi,1);
							
				// computing voltage and fvoltage
				if Index(nextVertex) lt Index(currentVertex) then
					voltage := voltageassignment(dart);
					if dart^f mod 2 eq 0 then
						fvoltage := voltageassignment(dart^f-1)^(-1);
					else
						fvoltage := voltageassignment(dart^f);
					end if;
				else
					voltage := voltageassignment(dart)^(-1);
					if dart^f mod 2 eq 0 then
						fvoltage := voltageassignment(dart^f-1);
					else
						fvoltage := voltageassignment(dart^f)^(-1);
					end if;
				end if;  
							
				potential[Index(nextVertex)] := voltage*potential[Index(currentVertex)];
				fpotential[Index(nextfVertex)] := fvoltage*fpotential[Index(currentfVertex)];
				for dart in multi do;                                                      // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
						if dart^f mod 2 eq 0 then
							fvoltage := (fpotential[Index(nextfVertex)])^(-1)*voltageassignment(dart^f-1)^(-1)*fpotential[Index(currentfVertex)];
						else
							fvoltage := (fpotential[Index(nextfVertex)])^(-1)*voltageassignment(dart^f)*fpotential[Index(currentfVertex)];
						end if;
					else
						voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment(dart)*potential[Index(nextVertex)];
						if dart^f mod 2 eq 0 then
							fvoltage := (fpotential[Index(currentfVertex)])^(-1)*voltageassignment(dart^f-1)^(-1)*fpotential[Index(nextfVertex)];
						else
							fvoltage := (fpotential[Index(currentfVertex)])^(-1)*voltageassignment(dart^f)*fpotential[Index(nextfVertex)];
						end if;
					end if; 
					Append(~A,voltage);
					Append(~B,fvoltage);
				end for;
			end if;
		end for;
		colour[Index(currentVertex)] := 2;                                                    //  searching of vertex is completed
		num := num + 1;
	end while;
		
	if num lt n then 
		error "Graph in first argument is not connected";
	end if;
	if not IsTransitive(PermutationGroup<Degree(G)|A>) then
		error "First argument is not locally transitive";
	end if;
	
	if VS`type eq "CayleyToPermutation" then
		solution := SolveSystem(A,B,1);
		if solution then
			S := [i : i in [1..Degree(Parent(A[1]))]];
			return solution, S;
		else
			return solution;
		end if;
	else
		S :=[];
		for i := 1 to Degree(Parent(A[1])) do
			if SolveSystem(A,B,i) then
				Append(~S,i);
			end if;
		end for;
		if IsEmpty(S) then
			return false, _;
		else
			return true, S;
		end if;
	end if;
	
end function;



//========================================================================================
//                 Testing whether two assignments are equivalent 
//========================================================================================

IsEquivalent := function(VS1, VS2)

// Given voltage spaces 'VS1'::Rec and 'VS2'::Rec of the same base graph this function tests whether projections are equivalent via solving a system of permutation equations.


	local X1, X2, G1, G2, voltageassignment1, voltageassignment2, V, n, A, B, colour, potential, fpotential, Q, currentVertex, nextVertex,
		   voltage, fvoltage, num, solution, S;
	
	X := VS1`graph;
	X2 := VS2`graph;
	
	if X ne X2 then
		error "Arguments 1 and 2 should should have common base graph";
	end if;
	
	if VS1`type eq "Cayley" then
		VS1 := CayleyToPermutation(VS1);
	end if;
	
	if VS2`type eq "Cayley" then
		VS2 := CayleyToPermutation(VS2);
	end if;

	G1 := VS1`group;
	G2 := VS2`group;
	
	if Degree(G1) ne Degree(G2) then
		return false;
	end if; 
	
	G := Sym(Degree(G1));
	voltageassignment1 := VS1`assignment;
	voltageassignment2 := VS2`assignment;
	
	V := VertexSet(X);
	n := #V;
	
	// construction of a system of permutation equations
	A := [];
	B := [];
	colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
	potential := [Id(G) : i in [1..n]];  // potential[i] is the voltage of the unique path on the tree from vertex V.i to vertex u
	fpotential := [Id(G) : i in [1..n]]; // fpotential[i] is the voltage of the unique path on the tree from vertex f(V.i) to vertex f(u)
	Q := [1];                            // vertices that has to be searched
	colour[1] := 1;                      // vertex u is coloured grey(1), since it has been searched at the moment
	num := 0;                            // number of searched vertices
			
	while IsEmpty(Q) eq false do         // BFS 
		currentVertex := V.(Q[1]);
		Remove(~Q,1);
		for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
			if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
				multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
				for dart in multi do;                                                      // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment1(dart)*potential[Index(currentVertex)];
						fvoltage := (fpotential[Index(nextVertex)])^(-1)*voltageassignment2(dart)*fpotential[Index(currentVertex)];
					else
						voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment1(dart)*potential[Index(nextVertex)];
						fvoltage := (fpotential[Index(currentVertex)])^(-1)*voltageassignment2(dart)*fpotential[Index(nextVertex)];
					end if; 
					Append(~A,voltage);
					Append(~B,fvoltage);
				end for;	
			elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
				multi := Indices(currentVertex, nextVertex);                                // all multiple edges
				dart := multi[1];                                                           // edge with the first index is tree edge
				colour[Index(nextVertex)] := 1; 
				Append(~Q,Index(nextVertex));
				Remove(~multi,1);
							
				// computing voltage and fvoltage
				if Index(nextVertex) lt Index(currentVertex) then
					voltage := voltageassignment1(dart);
					fvoltage := voltageassignment2(dart);
				else
					voltage := voltageassignment1(dart)^(-1);
					fvoltage := voltageassignment2(dart)^(-1);
				end if;  
							
				potential[Index(nextVertex)] := voltage*potential[Index(currentVertex)];
				fpotential[Index(nextVertex)] := fvoltage*fpotential[Index(currentVertex)];
				for dart in multi do;                                                      // cotree edges
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment1(dart)*potential[Index(currentVertex)];
						fvoltage := (fpotential[Index(nextVertex)])^(-1)*voltageassignment2(dart)*fpotential[Index(currentVertex)];
					else
						voltage := (potential[Index(currentVertex)])^(-1)*voltageassignment1(dart)*potential[Index(nextVertex)];
						fvoltage := (fpotential[Index(currentVertex)])^(-1)*voltageassignment2(dart)*fpotential[Index(nextVertex)];
					end if; 
					Append(~A,voltage);
					Append(~B,fvoltage);
				end for;
			end if;
		end for;
		colour[Index(currentVertex)] := 2;                                                    //  searching of vertex is completed
		num := num + 1;
	end while;
		
	if num lt n then 
		error "Base graph is not connected";
	end if;
	if not IsTransitive(PermutationGroup<Degree(G1)|A>) then
		error "First argument is not locally transitive";
	end if;
	
	solution := false;
	for i := 1 to Degree(Parent(A[1])) do
		if SolveSystem(A,B,i) then
			solution := true;
			break;
		end if;
	end for;
	return solution;
end function;






//========================================================================================
//                      Computing the lift 
//========================================================================================

Lift := function(X,Y,p,f,u,v)

/* for a given covering projection 'p'::Map mapping from covering graph 'X'::GrphMultUnd to base graph 'Y'::GrphMultUnd and an automorphism 
 * 'f'::GrpPermElt  of the base graph this function constructs a lift of f such that a vertex 'u'::GrphVert of X is mapped to 
 * a vertex 'v'::GrphVert of X;
 * notice: graphs X and Y can be generalized graphs - the information about the semiedges is encoded in projection p given on darts; 
 */
	local V, n, E, colour,liftOnDarts, liftOnVertex, Q, currentVertex, nextVertex, multi, dart, find, fcurrentVertex, fnextVertex, s, canlift,
	      num;
	 
	V := VertexSet(X);
	n := #V;	
	E := EdgeSet(X);
	colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
	liftOnDarts := [0 : i in [1..2*#E]];  // lift on darts
	liftOnVertex := [0 : i in [1..n]];   // lift on vertex
	ispermutation := [0 : i in [1..n]];  // tells which vertices are already images
	ispermutation[Index(v)] := 1;
	liftOnVertex[Index(u)] := Index(v);
	Q := [Index(u)];                     // vertices that has to be searched
	canlift := true;
	num := 0;                            // number of searched vertices
	
	colour[Index(u)] := 1;               // vertex u is coloured grey(1), since it has been searched at the moment
		
	while (IsEmpty(Q) eq false) and canlift do         // BFS 
		currentVertex := V.(Q[1]);
		fcurrentVertex := V.(liftOnVertex[Index(currentVertex)]);
		Remove(~Q,1);
		for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
			if colour[Index(nextVertex)] eq 1 then                                         // cotree edge, the image of nextVertex is known
				fnextVertex := V.(liftOnVertex[Index(nextVertex)]);
				multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices
				fmulti := SetToSequence(SequenceToSet(Indices(fcurrentVertex, fnextVertex)));
				for dart in multi do
					find := false;
					for i in fmulti do
						if (p(i) eq p(dart)^f) and (liftOnVertex[Index(InitialVertex(E.dart))] eq Index(InitialVertex(E.i))) then
							find := true;
							liftOnDarts[dart] := i;
							liftOnDarts[dart + 1] := i + 1;
							s := i;
							break;
						end if;
						if (p(i + 1) eq p(dart)^f) and (liftOnVertex[Index(InitialVertex(E.dart))] eq Index(TerminalVertex(E.i))) then
							find := true;
							liftOnDarts[dart] := i + 1;
							liftOnDarts[dart + 1] := i;
							s := i;
							break;
						end if;
					end for;
					
					if find eq false then
						canlift := false;
						break;
					else
						Exclude(~fmulti, s);
					end if;
					 
				end for;
			elif colour[Index(nextVertex)] eq 0 then                                        // tree edge, define new image
				colour[Index(nextVertex)] := 1; 
				Append(~Q,Index(nextVertex));
				multi := Indices(currentVertex, nextVertex);                                // all multiple edge	
				dart := multi[1];                                                           // edge with the first index is tree edge
				// find fnextvertex
				find := false;
				for j in Neighbours(fcurrentVertex) do
					fmulti := SetToSequence(SequenceToSet(Indices(fcurrentVertex, j)));
					for i in fmulti do
						if p(dart)^f eq p(i)  then
							if (InitialVertex(E.dart) eq currentVertex) and (InitialVertex(E.i) eq fcurrentVertex) then
								fnext := j;
								find := true;
								break;
							end if;
							if (InitialVertex(E.dart) eq nextVertex) and (InitialVertex(E.i) eq j) then
								fnext := j;
								find := true;
								break;
							end if;
						end if;
						if p(dart)^f eq p(i+1) then
							if (InitialVertex(E.dart) eq currentVertex) and (TerminalVertex(E.i) eq fcurrentVertex) then
								fnext := j;
								find := true;
								break;	
							end if;
							if (InitialVertex(E.dart) eq nextVertex) and (TerminalVertex(E.i) eq j) then
								fnext := j;
								find := true;
								break;
							end if;
						end if;
					end for;
					
					if find eq true then
						break;
					end if;	
				end for;
				
				if find eq true then
					liftOnVertex[Index(nextVertex)] := Index(fnext);
					if ispermutation[liftOnVertex[Index(nextVertex)]] eq 1 then
						canlift := false;
					else                                                          // good on mapped vertex, check potential multiple edges
						ispermutation[liftOnVertex[Index(nextVertex)]] := 1;
						fnextVertex := V.(liftOnVertex[Index(nextVertex)]);
						multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); 
						fmulti := SetToSequence(SequenceToSet(Indices(fcurrentVertex, fnextVertex)));
						for dart in multi do
							find := false;
							for i in fmulti do
								if (p(i) eq p(dart)^f) and (liftOnVertex[Index(InitialVertex(E.dart))] eq Index(InitialVertex(E.i))) then
									find := true;
									liftOnDarts[dart] := i;
									liftOnDarts[dart + 1] := i + 1;
									s := i;
									break;
								end if;
								if (p(i + 1) eq p(dart)^f) and (liftOnVertex[Index(InitialVertex(E.dart))] eq Index(TerminalVertex(E.i))) then
									find := true;
									liftOnDarts[dart] := i + 1;
									liftOnDarts[dart + 1] := i;
									s := i;
									break;
								end if;
							end for;
										
							if find eq false then
								canlift := false;
								break;
							else
								Exclude(~fmulti, s);
							end if;
											 
						end for;
					end if;
				else
					error "First argument is not a cover of second argument or fourth argument is not an automorphism of second argument";
				end if;
			end if;
			
			if canlift eq false then
				break;
			end if;
			
		end for;
		colour[Index(currentVertex)] := 2;                                                    //  searching of vertex is completed
		num := num + 1;
	end while;
	
	if canlift then
		if num lt n then
			error "First argument is not connected graph";
		else 
			return Sym(2*#E)!liftOnDarts, Sym(#V)!liftOnVertex;
		end if;
	else
		error "Fourth argument has no lift that maps fifth argument to sixth argument";
	end if;		
			
end function;

//========================================================================================
//             Computing the automorphism group of the generalised graph 
//========================================================================================

AutGroup := function(X,semi)

/* given multigraph 'X'::GrphMultUnd together with additional information 'semi'::Map[SetIndx[RngIntElt],SetIndx[BoolElt]] about semi edges this
 * function computes automorphism group of X;
 * automorphism group is given as permutation group of degree 2*e, where e is the number of the edges in X(it is assumed
 * that indices of the edges are in the range [1..2*e] );
 * each automorphism(element of the group) is given as a permutation on darts, when two contiguous darts are identified(representing semiedge)
 * then an odd dart is always mapped to an odd dart and vice versa;
 */

	local BipartGr, V, s, E, newE, S, dartToBlock, blockToDart, dart1, dart2, block1, block2, L, G, stab, autonDarts, P;

	
	V := VertexSet(X);
	E := EdgeSet(X);
	newE := [];
	stab := {};
	dartToBlock := [];
	blockToDart := [];
	s := 1;
	
	for i in Indices(X) do
		if semi(i) eq false then           // not semiedge
			Append(~newE, <Index(InitialVertex(E.i)),{Integers()!(#Vertices(X) + (i+1)/2)}>);
			Append(~newE, <Index(TerminalVertex(E.i)),{Integers()!(#Vertices(X) + (i+1)/2)}>);
			Append(~dartToBlock,s);
			Append(~dartToBlock,s+1);
			Append(~blockToDart,i);
			Append(~blockToDart,i+1);
			s := s + 2;
		else
			Append(~newE, <Index(InitialVertex(E.i)),{Integers()!(#Vertices(X) + (i+1)/2)}>);
			Append(~dartToBlock,s);
			Append(~dartToBlock,0);
			Append(~blockToDart,i);
			s := s + 1;
		end if;
	end for;
	
	BipartGr := MultiGraph<(#Vertices(X) + #Edges(X))| newE>;
	
	L := [];
	
	for i in Indices(BipartGr) do
		Append(~L,{Index(InitialVertex(EdgeSet(BipartGr).i)),Index(TerminalVertex(EdgeSet(BipartGr).i))});
	end for;
	
	S := IncidenceStructure<(#Vertices(X) + #Edges(X))|L>;
	
	for i in Vertices(X) do
		Include(~stab,Index(i));
	end for;
	
	if MinimumDegree(X) gt 2 then
		G := AutomorphismGroup(S);
	else
		G := Stabilizer(AutomorphismGroup(S),stab);
	end if;	
	
	P := [];
	
	if IsSimple(S) then
		for gen in Generators(G) do
				autOnDarts := [];               // automorphism on darts
				for i in Indices(X) do
					if semi(i) eq false then  
						block1 := Blocks(S)[dartToBlock[i]];
						block2 := Blocks(S)[dartToBlock[i+1]];
						dart1 := blockToDart[Index(Blocks(S),block1^gen)];
						dart2 := blockToDart[Index(Blocks(S),block2^gen)];	
						Append(~autOnDarts,dart1);
						Append(~autOnDarts,dart2);
					else
						block1 := Blocks(S)[dartToBlock[i]];
						dart1 := blockToDart[Index(Blocks(S),block1^gen)];
						Append(~autOnDarts,dart1);
						Append(~autOnDarts,dart1+1);
					end if;
				end for;
				Append(~P,Sym(2*#E)!autOnDarts);
			end for;
	
	else	
		for gen in Generators(G) do
			autOnDarts := [];               // automorphism on darts
			for i in Indices(X) do
				if semi(i) eq false then  
					block1 := dartToBlock[i];
					block2 := dartToBlock[i+1];
					dart1 := blockToDart[(#Points(S) + block1)^gen - #Points(S)];
					dart2 := blockToDart[(#Points(S) + block2)^gen - #Points(S)];
					Append(~autOnDarts,dart1);
					Append(~autOnDarts,dart2);
				else
					block1 := dartToBlock[i];
					dart1 := blockToDart[(#Points(S) + block1)^gen - #Points(S)];
					Append(~autOnDarts,dart1);
					Append(~autOnDarts,dart1+1);
				end if;
			end for;
			Append(~P,Sym(2*#E)!autOnDarts);
		end for;
	end if;
	
	return PermutationGroup<2*#E|P>;
	
end function;

//========================================================================================
//             Construction of the group of covering transformations 
//========================================================================================

GroupOfCoveringTransformations := function(VS)

/* given a voltage space 'VS'::Rec this function constructs a group of covering transformation CT; 
 * group CT is given as a permutation group on darts of the  derived cover;
 */

	local Central, X, V, E, u, loc, DerGraph, semi, p, CT, u1, v1;
	
	if VS`type eq "Cayley" then
		VS := CayleyToPermutation(VS);
	end if;
		
	X := VS`graph;
	V := VertexSet(X);
	E := EdgeSet(X);
	u := V.1;
	n := Degree(VS`group);
	
	
	// compute centralizer of the local voltage group
	loc := LocalVoltageGroup(VS, u);
	Central := Centralizer(Sym(n),loc);
	
	
	DerGraph,semi,p:= DerivedCover(VS);
	if Order(Central) eq 1 then
		return PermutationGroup<2*#EdgeSet(DerGraph)|>;
	else
		
		id := Id(Sym(2*#E));
		CT := [];
		u1 := VertexSet(DerGraph)!<u,1>;
		for gen in Generators(Central) do
			v1 := VertexSet(DerGraph)!<u,1^gen>;
			lift := Lift(DerGraph,X,p,id,u1,v1);
			Append(~CT, lift);
		end for;
		return PermutationGroup<2*#EdgeSet(DerGraph)|CT>;
	end if;

end function;

GroupOfCT := function(VS, DerGraph, semi, p)

/* given a voltage space 'VS'::Rec this function constructs a group of covering transformation CT; 
 * group CT is given as a permutation group on darts of the derived cover (DerGraph, semi);
 */

	local Central, X, V, E, u, loc, CT, u1, v1;
	
	if VS`type eq "Cayley" then
		VS := CayleyToPermutation(VS);
	end if;
	
	X := VS`graph;
	V := VertexSet(X);
	E := EdgeSet(X);
	u := V.1;
	n := Degree(VS`group);
	
	// compute centralizer of the local voltage group
	loc := LocalVoltageGroup(VS, u);
	Central := Centralizer(Sym(n),loc);
	
	if Order(Central) eq 1 then
		return PermutationGroup<2*#EdgeSet(DerGraph)|>;
	else
		
		id := Id(Sym(2*#E));
		CT := [];
		u1 := VertexSet(DerGraph)!<u,1>;
		for gen in Generators(Central) do
			v1 := VertexSet(DerGraph)!<u,1^gen>;
			lift := Lift(DerGraph,X,p,id,u1,v1);
			Append(~CT, lift);
		end for;
		return PermutationGroup<2*#EdgeSet(DerGraph)|CT>;
	end if;

end function;

//========================================================================================
//               The structure of the lifted group 
//========================================================================================


/////////////// Testing whether the lifted group is a split extension ////////////////////////////

Evaluation_Non := function(potential, Gr_potential, MatRep, word, FGr, G, iso, X)

	ver := VertexSet(X).1;
	
	if IsElementaryAbelian(G) then
		p := Integers()!Order(G.1);
		V := RSpace(GF(p),#Generators(G));
		M := RMatrixSpace(GF(p),#Generators(G),#Generators(G)*#Gr_potential);
	else
		V := RSpace(Integers(),#Generators(G));
		M := RMatrixSpace(Integers(),#Generators(G),#Generators(G)*#Gr_potential);
	end if;

	A := M!0;
	b := V!Eltseq(Id(G));

	for k in [1..#Eltseq(word)] do
		if Eltseq(word)[k] lt 0 then
			o := Order(iso(FGr.(Abs(Eltseq(word)[k]))));
			mat := MatRep[Abs(Eltseq(word)[k])];
			for j in [1..o-1] do
				zetaQ := potential[Index(ver)];
				mapped := AutOnVertex(X,iso(FGr.(Abs(Eltseq(word)[k]))),ver);
				zetagQ := Gr_potential[Abs(Eltseq(word)[k])][Index(mapped)];

				A := Transpose(mat)*A;

				for l in [1..#Generators(G)] do
					A[l,(Abs(Eltseq(word)[k]) - 1)*(#Generators(G)) + l] := A[l,(Abs(Eltseq(word)[k]) - 1)*(#Generators(G)) + l] + 1;
				end for;
				// update vector b
				b := V!(b*mat) + V!(V!Eltseq(zetaQ)*mat) - V!Eltseq(zetagQ);
				ver := mapped;
			end for;
		else
			mat := MatRep[Abs(Eltseq(word)[k])];
			zetaQ := potential[Index(ver)];
			mapped := AutOnVertex(X,iso(FGr.(Abs(Eltseq(word)[k]))),ver);
			zetagQ := Gr_potential[Abs(Eltseq(word)[k])][Index(mapped)];
			
			A := Transpose(mat)*A;

			for l in [1..#Generators(G)] do
				A[l,(Abs(Eltseq(word)[k]) - 1)*(#Generators(G)) + l] := A[l,(Abs(Eltseq(word)[k]) - 1)*(#Generators(G)) + l] + 1;
			end for;
			// update vector b
			b := V!(b*mat) + V!(V!Eltseq(zetaQ)*mat) - V!Eltseq(zetagQ);
			ver := mapped;
		end if;
	end for;

	return A,b;
end function;

IsSplitAdmissibleEA := function(VS, Gr, FGr, iso, com)

/* given an elementary abelian regular Cayley voltage space 'VS'::Rec and a group of automorphisms 'Gr'::GrpPerm of the base graph such that Gr lifts this function
 * tests whether the lifted group is a split extension of the group of covering transformation by the group Gr via solving a system of linear equations; 
 * FGr is a group presentation of Gr, iso is an isomorphism from FGr to Gr,
 * if parameter com is true, then all potential solutions  up to conjugacy are also constructed
 */

	local u, loc, l, hom, D, newvoltages, X, G, voltageassignment, V, n, f, voltOfCycles, Gr_voltOfCycles , f_voltOfCycles , colour, potential, 
	      f_potential, Q, currentVertex, currentfVertex, nextVertex, nextfVertex, voltage, fvoltage, MatRep, Ve, P, C, M, e, mat, sol, im, Im, Vect, Syst, A, b, x, y, E, block, row, Bmatrix, p, VG, map, perm, Bound, L, Sol, Com, Sub, BJoin;
	
	X := VS`graph;
	V := VertexSet(X);
	n := #V;
	
	G := VS`group;
	voltageassignment := VS`assignment;
	
	// construction of potentials and voltages of fundamental cycles together with their g-images 
	
	Gr_voltOfCycles := [];
	Gr_potential := [];		
				
	for k in [1..#Generators(Gr)] do
		f := Gr.k;
		voltOfCycles := [];
		f_voltOfCycles := [];
		colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
		potential := [Id(G) : i in [1..n]];  // potential[i] is the voltage of the unique path on the tree from vertex V.i to vertex u
		fpotential := [Id(G) : i in [1..n]]; // f_potential[i] is the voltage of the unique path on the tree from vertex f(V.i) to vertex f(u)
		Q := [1];                            // vertices that has to be searched
		colour[1] := 1;                      // vertex u is coloured grey(1), since it has been searched at the moment
		
		while IsEmpty(Q) eq false do         // BFS 
			currentVertex := V.(Q[1]);
			currentfVertex := AutOnVertex(X,f,currentVertex);
			Remove(~Q,1);
			for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
				nextfVertex := AutOnVertex(X,f,nextVertex);
				if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
					multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
					for dart in multi do                                                       // cotree edges
						if Index(nextVertex) lt Index(currentVertex) then
							voltage := -(potential[Index(nextVertex)]) + voltageassignment(dart) + potential[Index(currentVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(nextfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(currentfVertex)];
							else
								fvoltage := - fpotential[Index(nextfVertex)] + voltageassignment(dart^f) + fpotential[Index(currentfVertex)];
							end if;
						else
							voltage := - potential[Index(currentVertex)] + voltageassignment(dart) + potential[Index(nextVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(currentfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(nextfVertex)];
							else
								fvoltage :=  - fpotential[Index(currentfVertex)] + voltageassignment(dart^f) + fpotential[Index(nextfVertex)];
							end if;
						end if; 
						if G!voltage notin sub<G|voltOfCycles> then
							Append(~voltOfCycles,voltage);
							Append(~f_voltOfCycles,fvoltage);
						end if;
					end for;	
				elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
					multi := Indices(currentVertex, nextVertex);                                // all multiple edges
					dart := multi[1];                                                           // edge with the first index is tree edge
					colour[Index(nextVertex)] := 1; 
					Append(~Q,Index(nextVertex));
					Remove(~multi,1);
								
					// computing voltage and fvoltage
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := voltageassignment(dart);
						if dart^f mod 2 eq 0 then
							fvoltage := - voltageassignment(dart^f-1);
						else
							fvoltage := voltageassignment(dart^f);
						end if;
					else
						voltage := - voltageassignment(dart);
						if dart^f mod 2 eq 0 then
							fvoltage := voltageassignment(dart^f-1);
						else
							fvoltage := - voltageassignment(dart^f);
						end if;
					end if;  
							
					potential[Index(nextVertex)] := voltage + potential[Index(currentVertex)];
					fpotential[Index(nextfVertex)] := fvoltage + fpotential[Index(currentfVertex)];
					for dart in multi do                                                       // cotree edges
						if Index(nextVertex) lt Index(currentVertex) then
							voltage := - potential[Index(nextVertex)] + voltageassignment(dart) + potential[Index(currentVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(nextfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(currentfVertex)];
							else
								fvoltage := - fpotential[Index(nextfVertex)] + voltageassignment(dart^f) + fpotential[Index(currentfVertex)];
							end if;
						else
							voltage := - potential[Index(currentVertex)] + voltageassignment(dart) + potential[Index(nextVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(currentfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(nextfVertex)];
							else
								fvoltage := - fpotential[Index(currentfVertex)] + voltageassignment(dart^f) + fpotential[Index(nextfVertex)];
							end if;
						end if; 
						if G!voltage notin sub<G|voltOfCycles> then
							Append(~voltOfCycles,voltage);
							Append(~f_voltOfCycles,fvoltage);
						end if;
						end for;
				end if;
			end for;
			colour[Index(currentVertex)] := 2;                                                    //  searching of the vertex is completed	
		end while;
		Append(~Gr_voltOfCycles,f_voltOfCycles);
		Append(~Gr_potential,fpotential);
	end for;

	// construct g# (as a matrix) for each generator in Gr
	MatRep := [];
	
	mat := [];
	for v in voltOfCycles do
		mat := mat cat Eltseq(v);
	end for;
			
	p := Integers()!Order(G.1);
	Ve := RSpace(GF(p),#Generators(G));
	C := RMatrixSpace(GF(p),#voltOfCycles,#Generators(G));
	M := RMatrixSpace(GF(p),#Generators(G),#Generators(G));	
		
	// for each generator G.i find coefficients x_i such that G.i = x_1*voltOfCycles[1]+...+x_r*voltOfCycles[r] e.i. solve linear system  
	sol := [];
	for i in  [1..#Generators(G)] do 
		e := [0 : i in [1..#Generators(G)]];
		e[i] := 1;
		Append(~sol, Solution(C!mat,Ve!e));
	end for;

	for k in [1..#Generators(Gr)] do
		f_voltOfCycles := Gr_voltOfCycles[k];
		Im := [];
		//compute the image of each generator in G of f
		for i in  [1..#Generators(G)] do 
			im := Id(G);
			for j in [1..#f_voltOfCycles] do
				im := im + Integers()!sol[i][j]*f_voltOfCycles[j];
			end for;
			Im := Im cat Eltseq(im);	
		end for;
		Append(~MatRep, M!Im); 
	end for;
			 
	// for each relation in Gr construct equation
	VG := RMatrixSpace(GF(p),#Generators(G),#Generators(G)*#Gr_potential);	
	A,b := Evaluation_Non(potential, Gr_potential, MatRep, LHS(Relations(FGr)[1])*(RHS(Relations(FGr)[1])^-1), FGr, G, iso, X);
	A := VG!A;
	b := Eltseq(b);
			
	for i in [2..#Relations(FGr)] do
		//evaluation of relator
		x,y := Evaluation_Non(potential, Gr_potential, MatRep, LHS(Relations(FGr)[i])*(RHS(Relations(FGr)[i])^-1), FGr, G, iso, X);
		A := VerticalJoin(A,VG!x);
		b := b cat Eltseq(y);
	end for;
	
	// solve the system
	Vect := RSpace(GF(p),#Generators(G)*#Relations(FGr));
	w,e,t := IsConsistent(Transpose(A), -Vect!b);
	
	if com then 
		if w then 
		
			//all complements up to conjugacy
		
			//construct space representing boundary cocycles
		
			Bound := MatRep;
			for j in [1..#Bound] do
				for i in [1..#Eltseq(Bound[j][1])] do
					Bound[j][i,i] := Bound[j][i,i] - 1;
				end for;
			end for;
		
			BJoin := Transpose(Bound[1]);
			for i in [2..#Bound] do
				BJoin := VerticalJoin(BJoin,Transpose(Bound[i]));
			end for; 
	
			L := [];
			for i in [1..#Eltseq(Bound[1][1])] do
				L := Append(L, Eltseq(Transpose(BJoin)[i]));
			end for;	
			
			Quo, mapquo := quo< t | L>;
			
			Sub := sub< t| L>;
			// solutions up to conjugacy
			Sol := [];
			for i in Quo do
				Sol := Append(Sol, Inverse(mapquo)(i) + e);
			end for;
			return w,Sol;
		else
			return w,_;
		end if;
	else
		return w,_;
	end if;
end function;


DecomposeQuotient := function(VS, N)

	local Gamma, GammaN, quotienMap,VSGammaN, D, voltagesGammaN;
	
	//construct new voltage space on base graph with quotient group Gamma/N
	
	Gamma := VS`group;
	GammaN, quotienMap := Gamma/N;
	VSGammaN := VS;
	//pretvori v abelski tip
	cosetMap, perm := CosetAction(GammaN, sub<GammaN|>);
	GAb, hom := AbelianGroup(perm);
	
	D := Domain(VS`assignment);
		
	voltagesGammaN := [];
	for i in D do
		Append(~voltagesGammaN ,hom(cosetMap(quotienMap(VS`assignment(i)))));
	end for;
		
	VSGammaN`group := GAb;
	VSGammaN`assignment := map<D -> GAb|{<D[i], voltagesGammaN [i]> : i in [1..#D]}>;
	
	return VSGammaN;
end function;


DerivedCoverEA := function(VS)

/* given a Cayley voltage space 'VS'::Rec with elementary abelian group a derived cover is
 * constructed together with projective map which maps from dart set of derived cover to dart set of base graph: if an edge e in derived cover 
 * is mapped to a semiedge s, then both respective darts of e are mapped to even dart of s;
 */
 
	local IsSemiedge, ind, X, G, voltageassignment, D, order, vertices, label,edges, n, f, F, semi, num, M, V, newsemiedge, projective, projectiveMap, g, I;	
	
	X := VS`graph;
	IsSemiedge := VS`semiedge;
	voltageassignment := VS`assignment;
	G := VS`group;
	vertices := {@ @};
	D := [i : i in [1..(#EdgeSet(X)*2)]];
	edges := [];
	semi := [];
	projective := [];
	label := [];
	order := [g : g in G];
	ind := map<G -> {i : i in [1..#order]}|{<order[i], i> : i in [1..#order]}>;
	n := #order;
	
	for v in VertexSet(X) do                                      // construction of vertex set of derived graph 
		for g in order do
			Include(~vertices,<v,g>);
			Append(~label,g);
		end for;
	end for;
	
	M := MultiGraph<vertices|>;
	V := VertexSet(M);
	
	for e in EdgeSet(X) do                                                     // construction of edge set of derived graph  with projective map
		if TerminalVertex(e) eq InitialVertex(e) then
			if IsSemiedge(Index(e)) then                                       // semiedge in base graph
				num := [1 : i in [1..n]];
				for i in [1..n] do
					if num[i] eq 1 then
						g := order[i];
						if voltageassignment(Index(e)) eq G.0 then             // semiedge in derived cover
							Append(~semi, true);
							Append(~edges, <<InitialVertex(e),g>,{<InitialVertex(e),g>}>);
							Append(~projective, Index(e));
							Append(~projective, Index(e));
						else                                                   // link in derived cover within the fibre
							Append(~semi, false);
							Append(~edges, <<InitialVertex(e),g>,{<InitialVertex(e),g + voltageassignment(Index(e))>}>);
							Append(~projective, Index(e));
							Append(~projective, Index(e));
							num[ind(g + voltageassignment(Index(e)))] := 0;
						end if;
					end if;
				end for;	
			else                                                               // loop in base graph
				for i in [1..n] do
					g := order[i];
					Append(~semi, false);
					Append(~edges, <<InitialVertex(e),g>,{<InitialVertex(e),g + voltageassignment(Index(e))>}>);
					if voltageassignment(Index(e)) eq G.0 then                 // loop in in derived cover
						Append(~projective, Index(e));
						Append(~projective, Index(e) + 1);
					else                                                       // link in derived cover within the fibre
						if Index(V!<InitialVertex(e),g>) lt Index(V!<InitialVertex(e), g + voltageassignment(Index(e))>) then
							Append(~projective, Index(e));
							Append(~projective, Index(e) + 1);
						else 
							Append(~projective, Index(e) + 1);
							Append(~projective, Index(e));
						end if;
					end if;
				end for;
			end if;
		else                                                                    // link in base graph
			for i in [1..n] do
				g := order[i];
				Append(~semi, false);
				Append(~edges, <<InitialVertex(e),g>,{<TerminalVertex(e), g + voltageassignment(Index(e))>}>);
				if Index(V!<InitialVertex(e),g>) lt Index(V!<TerminalVertex(e), g + voltageassignment(Index(e))>) then
					Append(~projective, Index(e));
					Append(~projective, Index(e) + 1);
				else
					Append(~projective, Index(e) + 1);
					Append(~projective, Index(e));
				end if;
			end for;
		end if;
	end for;
	
	M := MultiGraph<vertices|edges>;
	I := [i : i in [1..2*#EdgeSet(M)]];
	projectiveMap := map<I -> D|{<I[i], projective[i]> : i in [1..2*#EdgeSet(M)]}>;
	newsemiedge := map<Indices(M) -> [true,false]|{<Indices(M)[i], semi[i]> : i in [1..#Indices(M)]}>;

	return  M, newsemiedge, projectiveMap, label, G;
	
end function;


DecomposeNormal := function(VS, VSQ, N)

	local X, Gamma ,GammaN, quotientMap,VSN, D, voltagesN, DCGraph, newsemiedge, projectiveMap, transversal, GAb, hom, dom,
	       cosetMap, mP, P, perm, transversalMap, darts, pvol, newvol, inivol, tervol, semi;

	//construct new voltage space on derived graph with normal group N
	
	X := VS`graph;
	Gamma := VS`group;
	GammaN, quotientMap := Gamma/N;
	cosetMap, perm := CosetAction(GammaN, sub<GammaN|>);
	GAb, hom := AbelianGroup(perm);
	
	DCGraph, newsemiedge, projectiveMap, label, G:= DerivedCoverEA(VSQ);
	transversal, transversalMap := Transversal(Gamma,N);
	
	darts := {@i : i in Indices(DCGraph)@};
	semi := [newsemiedge(i) : i in Indices(DCGraph)];
	// compute voltage assignment
	voltagesN := [];
	for i in Indices(DCGraph) do
		// voltage pvol of the projected dart i
		if projectiveMap(i) mod 2 eq 0 then
			pvol := VS`assignment(projectiveMap(i))^(-1);
		else
			pvol := VS`assignment(projectiveMap(i));
		end if;
		inivol := transversalMap(Inverse(quotientMap)(Inverse(cosetMap)(Inverse(hom)(GAb!Eltseq(label[Index(InitialVertex(EdgeSet(DCGraph).i))])))));
		tervol := transversalMap(Inverse(quotientMap)(Inverse(cosetMap)(Inverse(hom)(GAb!Eltseq(label[Index(TerminalVertex(EdgeSet(DCGraph).i))])))));
		newvol := N!(inivol*pvol*(tervol^(-1)));
		Append(~voltagesN , newvol);
	end for;
	
	VSN := VoltageSpace("Cayley", N, DCGraph, darts, semi, voltagesN);
	
	return VSN, G, DCGraph, projectiveMap;
end function;
	

ComplementConstruction := function(A, B, VSA, compsdown, CTdown, FGr)
 
	VSAmodB := DecomposeQuotient(VSA, B);
	allcompsup := [];
	constructed := false;
	
	for com in compsdown do
		C := PermutationGroup< #Eltseq(com[1])|com>;
		iso := hom<FGr -> C| com>;
		w,Sol := IsSplitAdmissibleEA(VSAmodB, C, FGr, iso, true);
		
		if w then
			if  not constructed then
				VSB, GAb, DCGraph, projectiveMap := DecomposeNormal(VSA, VSAmodB, B);
				Gen := [];
				id := Id(Sym(2*#Edges(VSAmodB`graph)));
				for gen in Generators(GAb) do
					idtilde,_:=Lift(DCGraph,VSAmodB`graph,projectiveMap,id,VertexSet(DCGraph)!<VertexSet(VSAmodB`graph).1,GAb.0>,VertexSet(DCGraph)!<VertexSet(VSAmodB`graph).1,gen>);
					Append(~Gen, idtilde);
				end for;	
				CTAmodB := PermutationGroup< #Eltseq(Gen[1])| Gen>;
				constructed := true;
			end if;
			// construct all solutions upstairs
			compsup := [];
			for e in Sol do
				gen := [];
				for j  in [1..#com] do
					g := com[j]; // j-th generator;
					dimens := #Eltseq(GAb.0); // dimension of abelian group
					k := GAb!([e[i] :i in [((j-1)*dimens + 1)..(j*dimens)]]);
					gtilde,_:=Lift(DCGraph,VSAmodB`graph,projectiveMap,g,VertexSet(DCGraph)!<VertexSet(VSAmodB`graph).1,GAb.0>,VertexSet(DCGraph)!<AutOnVertex(VSAmodB`graph, g, VertexSet(VSAmodB`graph).1),k>);
					Append(~gen, gtilde);
				end for;
				Append(~compsup, gen);
			end for;

			//reducing modulo normalizer
			
			Gdown := PermutationGroup<Degree(C)|CTdown,C>;
			Nor := Normaliser(Gdown,sub<Gdown|C>);
			NorinCTdown := Nor meet CTdown;
			lifttran := [];
			if #Generators(NorinCTdown) ne 0 and #compsup gt 1 then
				for j  in [1..#Generators(NorinCTdown)] do
					g := NorinCTdown.j; // j-th generator;
					gtilde,_:=Lift(DCGraph,VSAmodB`graph,projectiveMap,g,VertexSet(DCGraph)!<VertexSet(VSAmodB`graph).1,GAb.0>,VertexSet(DCGraph)!								<AutOnVertex(VSAmodB`graph, g, VertexSet(VSAmodB`graph).1),GAb.0>);
					Append(~lifttran, gtilde);
				end for;
				liftedcom := [];
				for gen in compsup[1] do
					Append(~liftedcom,Eltseq(gen));
				end for;	
				for gen in Generators(CTAmodB) do
					Append(~liftedcom,Eltseq(gen));
				end for;
				// simulate action of normaliser on H(Gamma,G)
				genAction := [];
				for gen in lifttran do
					gn := [];
					for j in [1..#compsup] do
						con := PermutationGroup< #Eltseq(compsup[j][1])|compsup[j]>^gen;
						for k in [1..#compsup] do
							if IsConjugate(PermutationGroup< #Eltseq(liftedcom[1])|liftedcom>, con, PermutationGroup< #Eltseq(compsup[k][1])|										compsup[k]>) then
								gn[j] := k;
								break;
							end if;
						end for;
					end for;
					Append(~genAction, gn);
				end for;
				ngorbs := OrbitRepresentatives( sub<Sym(#compsup)|genAction> );
				compreduced := [ compsup[ngorbs[k][2]] : k in [1..#ngorbs] ];
				allcompsup := allcompsup cat compreduced;
			else
				allcompsup := allcompsup cat compsup;
			end if;
		end if;
	end for;
	
	if #allcompsup gt 0 then
		// construct also lift of CTdown along VSAmodB
		gen := [];
		if #Generators(CTdown) ne 0 then
			for j  in [1..#Generators(CTdown)] do
				g := CTdown.j; // j-th generator;
				gtilde,_:=Lift(DCGraph,VSAmodB`graph,projectiveMap,g,VertexSet(DCGraph)!<VertexSet(VSAmodB`graph).1,GAb.0>,VertexSet(DCGraph)!<AutOnVertex(VSAmodB`graph, g, VertexSet(VSAmodB`graph).1),GAb.0>);
				Append(~gen, gtilde);
			end for;
		end if;
		for g in Generators(CTAmodB) do
			Append(~gen, g);
		end for;
		CTup := PermutationGroup< #Eltseq(gen[1])|gen>;
		
		return allcompsup, VSB,  CTup;
	else
		return allcompsup, _,_;
	end if;
end function;


LastStep := function(A, B, VSA, compsdown, CTdown, FGr)

	VSAmodB := DecomposeQuotient(VSA, B);
	allcompsup := [];
	
	issplit := false;
	for com in compsdown do
		iso := hom<FGr -> PermutationGroup< #Eltseq(com[1])|com>| com>;
		w,_:= IsSplitAdmissibleEA(VSAmodB, PermutationGroup< #Eltseq(com[1])|com>, FGr, iso,false);
		if w then
			issplit := true;
			break;
		end if;
	end for;
	
	return issplit;
end function;


IsSplitAdmissible := function(VS, Gr)

/* given a solvable regular covering projection 'VS'::Rec of connected graphs along which a given group of automorphisms 'Gr'::GrpPerm of the base graph
 * lifts this function tests whether the lifted group is split extension of the group of covering transformation by the group Gr;
 * if voltage space 'VS'::Rec is of type 'Cayley', it is assumed that VS is locally transitive with solvable voltage group; 
 * if voltage space 'VS'::Rec is of type 'permutation', it is assumed that the corresponding local group is regular and solvable;
 */

	local u, loc, l, hom, map, perm, D, newvoltages, G, VSGPC, series, VSnew, VSQ, w, issplit, e, gen, g, gtilde, k, DCGraph,
	      dimens, projectiveMap, f, F;

	// convert to Cayley voltage space with voltage group of type GrpGPC or GrAb if needed
	if VS`type ne "Cayley" then
		u := Vertices(VS`graph)[1];
		loc := LocalVoltageGroup(VS,u);
		if not IsRegular(loc) then
			error "First argument should be connected regular voltage space";
		end if;
		if not IsSolvable(loc) then
			error "Group of covering transformations should be solvable";
		else
			//convert to Cayley 
			VSGPC := VS;
			VSGPC`assignment, _ := Reduced(VS, u);
			VSGPC`type := "Cayley";
	
			if IsElementaryAbelian(loc) then
				G, hom := AbelianGroup(loc);
			else
				G, hom := GPCGroup(loc);
			end if;
		
			VSGPC`group := G;
			D := Domain(VS`assignment);	
			newvoltages := [];
			for i in D do
				Append(~newvoltages,hom(VSGPC`assignment(i)));
			end for;
			VSGPC`assignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
		end if;
	else 
		if IsLocallyTransitive(VS) then
			if Type(VS`group) eq GrpFP then
				map, loc := CosetAction(VS`group,sub<VS`group|>);
				D := Domain(VS`assignment);
				newvoltages := [];
				for i in D do
					Append(~newvoltages,map(VS`assignment(i)));
				end for;
				VS`group := loc;
				VS`assignment := map<D -> loc|{<D[i], newvoltages[i]> : i in [1..#D]}>;	
			end if;
			if not IsSolvable(VS`group) then
					error "Group of covering transformations should be solvable";
			else
				//Cayley voltage
				VSGPC := VS;
				loc := VS`group;
				if IsElementaryAbelian(loc) then
					if Type(VS`group) ne GrpAb then
						G, hom := AbelianGroup(loc);
						D := Domain(VS`assignment);	
						newvoltages := [];
						for i in D do
							Append(~newvoltages,hom(VS`assignment(i)));
						end for;
						VSGPC`group := G;
						VSGPC`assignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
					end if;
				else
					if Type(loc) ne GrpGPC then
						G, hom := GPCGroup(loc);
						D := Domain(VS`assignment);	
						newvoltages := [];
						for i in D do
							Append(~newvoltages,hom(VS`assignment(i)));
						end for;
						VSGPC`group := G;
						VSGPC`assignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
					end if;
				end if;
				
			end if; 
		else
			error "First argument should be connected regular voltage space";
		end if;
	end if;
	

	/*
	// check whether group Gr lifts
	for gen in Generators(Gr) do
		l, _ := HasLift(VS,gen);
		if not l then
			error "Second argument should lift";
		end if;
	end for;
	*/
	
	timeInit := Cputime();
	FGr, iso := FPGroup(Gr);
	
	if IsElementaryAbelian(loc) then
		issplit,_:= IsSplitAdmissibleEA(VSGPC, Gr, FGr, iso,false);
	else	
		VSA := VSGPC;
		issplit := true;
		series := EFASeries(VSGPC`group);
		CT := PermutationGroup<2*#Edges(VS`graph)| >;
		c := [];
		for i in [1..#Generators(Gr)] do
			Append(~c, iso(FGr.i));
		end for;
		com := [c];
		for i in [1..#series - 1] do
			// if i is last level then do not construct complements upstairs
			if i eq (#series - 1 ) then
				issplit := LastStep(series[i],series[i+1] , VSA, com, CT, FGr);
			else
				com, VSA, CT := ComplementConstruction(series[i], series[i+1] , VSA, com, CT, FGr);
				if #com eq 0 then
					issplit := false;
					break;
				end if;
			end if;
		end for;	
	end if;
	ti := Cputime(timeInit);
	return issplit, ti;
end function;



IsSplitAdmissibleNaive := function(VS,G)
	timeInit := Cputime();
	Xtilde,semi,proj := DerivedCover(VS);
	CT := GroupOfCT(VS, Xtilde,semi,proj); 
	S := Sym(#Eltseq(CT.1));
	genSym := {};
	for gen in Generators(CT) do
		genSym := Include(genSym,S!gen);
	end for;
	
	X := VS`graph;
	
	generatorsGtilde := genSym;
	
	for gen in Generators(G) do
		lift,_ := Lift(Xtilde,X,proj,gen,VertexSet(Xtilde)!<VertexSet(X).1,1>,VertexSet(Xtilde)!<AutOnVertex(X, gen, VertexSet(X).1),1>);
		generatorsGtilde := Include(generatorsGtilde,S!lift);
	end for;
	 Gtilde := PermutationGroup<#Eltseq(CT.1)|generatorsGtilde>;
	 //timecom := Cputime();
	 a,b:=HasComplement(Gtilde,CT);
	 //tt:= Cputime(timecom);
	 ti := Cputime(timeInit);
	return a, ti; 
end function;




Cocycle := function(VS, Gr, words, FGr, iso)
	
	/* computes partial 2-cocylce on words along with an "action" (in fact it is an action on Out(Gamma)) Gr -> Aut(Gammma);
	 * for each generator g in Gr g^# is given as a map: Gamma to Gamma
	 * word in words is given as a sequence of positive integers representing generators of Gr 
	 */	
	
	local X, V, n, G, voltageassignment, Gr_voltOfCycles, Gr_potential, f, voltOfCycles, f_voltOfCycles, colour, fpotential, Q, currentVertex,
		  currentfVertex,  nextfVertex, multi, voltage, fvoltage, dart, cocycles; 
	
	// construction of potentials and voltages of fundamental cycles together with their g-images
	X := VS`graph;
	V := VertexSet(X);
	n := #V;
	
	G := VS`group;
	voltageassignment := VS`assignment;
	
	Gr_voltOfCycles := [];
	Gr_potential := [];		
				
	for k in [1..#Generators(Gr)] do
		f := Gr.k;
		voltOfCycles := [];
		f_voltOfCycles := [];
		colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
		potential := [Id(G) : i in [1..n]];  // potential[i] is the voltage of the unique path on the tree from vertex V.i to vertex u
		fpotential := [Id(G) : i in [1..n]]; // f_potential[i] is the voltage of the unique path on the tree from vertex f(V.i) to vertex f(u)
		Q := [1];                            // vertices that has to be searched
		colour[1] := 1;                      // vertex u is coloured grey(1), since it has been searched at the moment
		
		while IsEmpty(Q) eq false do         // BFS 
			currentVertex := V.(Q[1]);
			currentfVertex := AutOnVertex(X,f,currentVertex);
			Remove(~Q,1);
			for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
				nextfVertex := AutOnVertex(X,f,nextVertex);
				if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
					multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
					for dart in multi do                                                       // cotree edges
						if Index(nextVertex) lt Index(currentVertex) then
							voltage := (potential[Index(nextVertex)])^(-1)*voltageassignment(dart)*potential[Index(currentVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := (fpotential[Index(nextfVertex)])^(-1)*(voltageassignment(dart^f-1))^(-1)* 
											fpotential[Index(currentfVertex)];
							else
								fvoltage := (fpotential[Index(nextfVertex)])^(-1)*voltageassignment(dart^f)*fpotential[Index(currentfVertex)];
							end if;
						else
							voltage := (potential[Index(currentVertex)])^(-1) * voltageassignment(dart) * potential[Index(nextVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := (fpotential[Index(currentfVertex)])^(-1)*(voltageassignment(dart^f-1))^(-1) *
										    fpotential[Index(nextfVertex)];
							else
								fvoltage := (fpotential[Index(currentfVertex)])^(-1) * voltageassignment(dart^f) * fpotential[Index(nextfVertex)];
							end if;
						end if; 
						if G!voltage notin sub<G|voltOfCycles> then
							Append(~voltOfCycles,voltage);
							Append(~f_voltOfCycles,fvoltage);
						end if;
					end for;	
				elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
					multi := Indices(currentVertex, nextVertex);                                // all multiple edges
					dart := multi[1];                                                           // edge with the first index is tree edge
					colour[Index(nextVertex)] := 1; 
					Append(~Q,Index(nextVertex));
					Remove(~multi,1);
								
					// computing voltage and fvoltage
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := voltageassignment(dart);
						if dart^f mod 2 eq 0 then
							fvoltage :=  (voltageassignment(dart^f-1))^(-1);
						else
							fvoltage := voltageassignment(dart^f);
						end if;
					else
						voltage := (voltageassignment(dart))^(-1);
						if dart^f mod 2 eq 0 then
							fvoltage := voltageassignment(dart^f-1);
						else
							fvoltage := (voltageassignment(dart^f))^(-1);
						end if;
					end if;  
							
					potential[Index(nextVertex)] := voltage * potential[Index(currentVertex)];
					fpotential[Index(nextfVertex)] := fvoltage * fpotential[Index(currentfVertex)];
					for dart in multi do                                                       // cotree edges
						if Index(nextVertex) lt Index(currentVertex) then
							voltage :=  (potential[Index(nextVertex)])^-1 * voltageassignment(dart) * potential[Index(currentVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := (fpotential[Index(nextfVertex)])^(-1) * (voltageassignment(dart^f-1))^(-1) *
											 fpotential[Index(currentfVertex)];
							else
								fvoltage := (fpotential[Index(nextfVertex)])^(-1) * voltageassignment(dart^f) * fpotential[Index(currentfVertex)];
							end if;
						else
							voltage := (potential[Index(currentVertex)])^(-1) * voltageassignment(dart) * potential[Index(nextVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := (fpotential[Index(currentfVertex)])^(-1) * (voltageassignment(dart^f-1))^(-1) *
											fpotential[Index(nextfVertex)];
							else
								fvoltage := (fpotential[Index(currentfVertex)])^(-1) * voltageassignment(dart^f) * fpotential[Index(nextfVertex)];
							end if;
						end if; 
						if G!voltage notin sub<G|voltOfCycles> then
							Append(~voltOfCycles,voltage);
							Append(~f_voltOfCycles,fvoltage);
						end if;
						end for;
				end if;
			end for;
			colour[Index(currentVertex)] := 2;                                                    //  searching of the vertex is completed	
		end while;
		Append(~Gr_voltOfCycles,f_voltOfCycles);
		Append(~Gr_potential,fpotential);
	end for;
	
	// F(g,h) = g^#_u_0(zeta(Q))zeta(g(Q))^-1; Q: h(u_0) -> u_0; if working with T-reduced voltage space F(g,h) = zeta(g(Q))^-1
	//compute homomorphisms g^# : Gamma -> Gamma

	homo := [];	
	for volt in Gr_voltOfCycles do
		homo := Append(homo,hom<G -> G| [<G!voltOfCycles[i],G!volt[i]>: i in [1..#voltOfCycles]]>);
	end for;
	
	ver := VertexSet(X).1;	
	cocycles := [];
	for word in words do
		cocycle := [];
		k := #word;
		for i in [1..(k-1)] do
			h := Gr.0;
			for j in [k-i+1..k] do
				h := iso(FGr.(word[j]))*h;
			end for;
			mapped := Index(AutOnVertex(X,h,ver)); // path Q
			mappedG := Index(AutOnVertex(X,iso(FGr.(word[k-i])),VertexSet(X).mapped));
			cocycle := Append(cocycle,homo[word[k-i]](potential[mapped])*(Gr_potential[word[k-i]][mappedG])^(-1));
		end for;
		cocycles := Append(cocycles,cocycle);
	end for;
	
	return cocycles, homo;
end function;


ConstructSystem := function(homoMatrix, cocyclesVe, words, p,d)

	M := RMatrixSpace(GF(p),d*#homoMatrix,d);	
	V := RSpace(GF(p),d);
	MM := RMatrixSpace(GF(p),d*#homoMatrix,0);	
	AA := MM!0;
	VV := RSpace(GF(p),d*#words);
	bb := [];
	for j in [1..#words] do	
		word := words[j];
		k := #word;		
		cycle := cocyclesVe[j];
		A := M!0;
		for l in [1..d] do
				A[(word[k] - 1)*d + l, l] := A[(word[k] - 1)*d + l, l] + 1;
			end for;
		b := V!0;
		
		for i in [1..(k-1)] do
			ind := k-i;
			A := A*homoMatrix[word[ind]];
			for l in [1..d] do
				A[(word[ind] - 1)*d + l, l] := A[(word[ind] - 1)*d + l, l] + 1;
			end for;
			b := b*homoMatrix[word[ind]] + cycle[i]; 
		end for;
		AA := HorizontalJoin(AA,A);
		bb := bb cat Eltseq(b);
	end for;
	
	return AA,VV!bb;
end function;


ExtendDerivation := function(word, derivation, cocycleVe, homoMatrix, d, p)
	Ve := RSpace(GF(p),d);
	M := RMatrixSpace(GF(p),d,d);	
	k := word[#word];
	mapped := Ve!([Integers()!(derivation[j]): j in [((k-1)*d + 1)..k*d]]);
	for i in [1..(#word -1)] do
		kk := word[#word - i]; 
		mapped := mapped*(homoMatrix[kk]) + Ve!([Integers()!(derivation[j]): j in [((kk-1)*d + 1)..kk*d]]) + cocycleVe[i]; 
	end for;
	return [Integers()!(Eltseq(mapped)[i]): i in [1..#Eltseq(mapped)]];
end function;




ComConstruction := function(COCYCLES, HOMO, words, FGr, iso, Gamma, N)
	
	NEWCOYCLES := <>;
	NEWHOMO := <>;
	for i in [1..#COCYCLES] do
		cocycles := COCYCLES[i]; homo := HOMO[i];
		GammaModN, q := Gamma/N; //elementary abelian group: translate to vector space  
		p := Integers()!Order(GammaModN.1);
		d := #Generators(GammaModN);
		Ve := RSpace(GF(p),d);
		M := RMatrixSpace(GF(p),d,d);	
	
		cocyclesVe := [];
		for cocycle in cocycles do
			newcocycle := [Ve!Eltseq(q(elt)) : elt in cocycle];
			cocyclesVe := Append(cocyclesVe, newcocycle);
		end for;
		homoMatrix := [];
		
		for h in homo do
			im := [];
			for i in [1..#Generators(GammaModN)] do
				im := Append(im,Eltseq(q(h((GammaModN!(Eltseq(GammaModN.i))) @@ q)))); 
			end for;
			homoMatrix := Append(homoMatrix, M!im);
		end for; 
		
		A, b := ConstructSystem(homoMatrix, cocyclesVe, words, p,d);
		w,e,t := IsConsistent(A, -b);	
		
		if w then 
		
			//all complements up to conjugacy
		
			//construct space representing boundary cocycles			
			B := RMatrixSpace(GF(p),d,0)!0;	
			
			for M in homoMatrix do
				MI := M;
				for i in [1..d] do
					MI[i,i] := MI[i,i] - 1;
				end for;
				B := HorizontalJoin(B,MI);
			end for;
			
			L := [];
			for i in [1..d] do
				L:= Append(L,Eltseq(B[i]));
			end for;
			
			Quo, mapquo := quo< t | L>;
			
			Sub := sub< t| L>;
			// solutions up to conjugacy
			Sol := [];
			for i in Quo do
				Sol := Append(Sol, Inverse(mapquo)(i) + e);
			end for;
			
			tran, tranMap := Transversal(Gamma,N);
			for sol in Sol do
				newhomo := [];
				newCYCLES := [];
				for i in [1..#homo] do
					preim := tranMap((GammaModN![Integers()!(sol[j]): j in [((i-1)*d + 1)..i*d]]) @@ q);
					newhomo :=  Append(newhomo, hom<N -> N| [<gen, preim*(homo[i](gen))* preim^(-1)>: gen in Generators(N)]>);
				end for;				
				NEWHOMO := Append(NEWHOMO, newhomo);
				
				for l in [1..#words] do
					word := words[l];
					newcocycle := [];
					k := #word;
					for i in [1..(k-1)] do
					//pozor, kaj pa Extend derivation
						pre := tranMap((GammaModN!ExtendDerivation([word[m]: m in [k-i..k]], sol, cocyclesVe[l],homoMatrix, d, p)) @@ q);
						pre1 := tranMap((GammaModN!ExtendDerivation([word[k-i]], sol, cocyclesVe[l],homoMatrix, d, p)) @@ q);						
						pre2 := tranMap((GammaModN!ExtendDerivation([word[m]: m in [k-i+1..k]], sol,cocyclesVe[l], homoMatrix, d, p)) @@ q);
						newcocycle := Append(newcocycle,pre1*(homo[word[k-i]](pre2))*(cocycles[l][i])*pre^(-1));
					end for;
					newCYCLES := Append(newCYCLES,newcocycle);
				end for;
				NEWCOYCLES := Append(NEWCOYCLES, newCYCLES);
			end for;
		end if;
	end for;	
	return NEWCOYCLES, NEWHOMO;
end function;


LastS := function(COCYCLES, HOMO, words, FGr, iso, Gamma, N)
	
	issplit := false;
	for i in [1..#COCYCLES] do
		cocycles := COCYCLES[i]; homo := HOMO[i];
		GammaModN, q := Gamma/N; //elementary abelian group: translate to vector space  
		p := Integers()!Order(GammaModN.1);
		d := #Generators(GammaModN);
		Ve := RSpace(GF(p),d);
		M := RMatrixSpace(GF(p),d,d);	
	
		cocyclesVe := [];
		for cocycle in cocycles do
			newcocycle := [Ve!Eltseq(q(elt)) : elt in cocycle];
			cocyclesVe := Append(cocyclesVe, newcocycle);
		end for;
		homoMatrix := [];
		
		for h in homo do
			im := [];
			for i in [1..#Generators(GammaModN)] do
				im := Append(im,Eltseq(q(h((GammaModN!(Eltseq(GammaModN.i))) @@ q)))); 
			end for;
			homoMatrix := Append(homoMatrix, M!im);
		end for; 
		
		A, b := ConstructSystem(homoMatrix, cocyclesVe, words, p,d);
		w,e,t := IsConsistent(A, -b);	
		
		if w then 
			issplit := true;
			break;
		end if;
		
	end for;
	return issplit;
end function;



IsSplitAdmissibleI := function(VS, Gr)

/* given a solvable regular covering projection 'VS'::Rec of connected graphs along which a given group of automorphisms 'Gr'::GrpPerm of the
 * base graph lifts this function tests whether the lifted group is split extension of the group of covering transformation by the group Gr;
 * if voltage space 'VS'::Rec is of type 'Cayley', it is assumed that VS is locally transitive with solvable voltage group; 
 * if voltage space 'VS'::Rec is of type 'permutation', it is assumed that the corresponding local group is regular and solvable;
 */

	local u, loc, l, hom, map, perm, D, newvoltages, G, VSGPC, series, VSnew, VSQ, w, issplit, e, gen, g, gtilde, k, DCGraph,
	      dimens, projectiveMap, f, F;

	// convert to Cayley voltage space with voltage group of type GrpGPC or GrAb if needed
	if VS`type ne "Cayley" then
		u := Vertices(VS`graph)[1];
		loc := LocalVoltageGroup(VS,u);
		if not IsRegular(loc) then
			error "First argument should be connected regular voltage space";
		end if;
		if not IsSolvable(loc) then
			error "Group of covering transformations should be solvable";
		else
			//convert to Cayley 
			VSGPC := VS;
			VSGPC`assignment, _ := Reduced(VS, u);
			VSGPC`type := "Cayley";
	
			//if IsElementaryAbelian(loc) then
			//	G, hom := AbelianGroup(loc);
			//else
				G, hom := GPCGroup(loc);
			//end if;
		
			VSGPC`group := G;
			D := Domain(VS`assignment);	
			newvoltages := [];
			for i in D do
				Append(~newvoltages,hom(VSGPC`assignment(i)));
			end for;
			VSGPC`assignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
		end if;
	else 
		if IsLocallyTransitive(VS) then
			if Type(VS`group) eq GrpFP then
				map, loc := CosetAction(VS`group,sub<VS`group|>);
				D := Domain(VS`assignment);
				newvoltages := [];
				for i in D do
					Append(~newvoltages,map(VS`assignment(i)));
				end for;
				VS`group := loc;
				VS`assignment := map<D -> loc|{<D[i], newvoltages[i]> : i in [1..#D]}>;	
			end if;
			if not IsSolvable(VS`group) then
					error "Group of covering transformations should be solvable";
			else
				//Cayley voltage
				VSGPC := VS;
				loc := VS`group;
				if IsElementaryAbelian(loc) then
					if Type(VS`group) ne GrpAb then
						G, hom := AbelianGroup(loc);
						D := Domain(VS`assignment);	
						newvoltages := [];
						for i in D do
							Append(~newvoltages,hom(VS`assignment(i)));
						end for;
						VSGPC`group := G;
						VSGPC`assignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
					end if;
				else
					if Type(loc) ne GrpGPC then
						G, hom := GPCGroup(loc);
						D := Domain(VS`assignment);	
						newvoltages := [];
						for i in D do
							Append(~newvoltages,hom(VS`assignment(i)));
						end for;
						VSGPC`group := G;
						VSGPC`assignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;
					end if;
				end if;
				
			end if; 
		else
			error "First argument should be connected regular voltage space";
		end if;
	end if;
	

	/*
	// check whether group Gr lifts
	for gen in Generators(Gr) do
		l, _ := HasLift(VS,gen);
		if not l then
			error "Second argument should lift";
		end if;
	end for;
	*/
	
	timeInit := Cputime();
	FGr, iso := FPGroup(Gr);

	//compute relators as sequences of positive integers
	
	words := [];
	for rel in Relations(FGr) do
		word := [];
		a := Eltseq(LHS(rel)*(RHS(rel)^-1));
		for int in a do
			if int lt 0 then
				o := Order(iso(FGr.(Abs(int)))); 
				for i in [1..o-1] do
					word := Append(word,Abs(int));
				end for;
			else
				word := Append(word,int);
			end if;
		end for;
		words := Append(words, Reverse(word));
	end for;
	
	cocycles, homo := Cocycle(VSGPC, Gr, words, FGr, iso);
	
	series := EFASeries(VSGPC`group);
	
	COCYCLES := <cocycles>;
	HOMO := <homo>; 
	issplit := true;
	
	for i in [1..#series - 1] do
		if i eq (#series - 1) then
			issplit := LastS(COCYCLES, HOMO, words, FGr, iso, series[i], series[i+1]);
		else
			COCYCLES, HOMO := ComConstruction(COCYCLES, HOMO, words, FGr, iso, series[i], series[i+1]);
			print #COCYCLES;
			print 0;
			if #COCYCLES eq 0 then
				issplit := false;
				break;
			end if;
		end if;
	end for;
	
	ti := Cputime(timeInit);
	return issplit, ti;
	
end function;



/////////////// Testing whether the lifted group is a strong (sectional) split extension ////////////////////////////

IsSectionalSplitAdmissible := function(VS, Gr, Omega)

/* given an abelian regular voltage space 'VS'::Rec, a group of automorphisms 'Gr'::GrpPerm of the base graph such that G lifts, and
 * Gr-invariant subset of vertices 'Omega'::SetIndx[GrphVert] of the base graph this function
 * tests whether the lifted group is split extension of the group of covering transformation by the group Gr with an invariant section
 * over Omega; 
 */

	local u, loc, l, hom, D, newvoltages, X, G, voltageassignment, V, n, f, voltOfCycles, Gr_voltOfCycles , f_voltOfCycles , colour, potential, 
	      f_potential, Q, currentVertex, currentfVertex, nextVertex, nextfVertex, voltage, fvoltage, MatRep, Ve, C, M, e, mat, sol, im, Im,
		  A, b, B, E, block, row, Bmatrix, base, Delta, IsDefined, numv, MS, Tr, zetaQ, mapped, zetagQ, Lab, lab, LabT, column, I, Iv, rowA, Esp, Vsp, ma;
	
	
	// check whether voltage space is regular and abelian
	if VS`type ne "Cayley" then
		u := Vertices(VS`graph)[1];
		loc := LocalVoltageGroup(VS,u);
		if not IsRegular(loc) then
			error "First argument should be regular voltage space";
		end if;
	end if;
	
	if not IsAbelian(CosetImage(VS`group,sub<VS`group|>)) then
		error "Voltage group should be abelian";
	end if;
	
	// check whether group Gr lifts
	for gen in Generators(Gr) do
		l, _ := HasLift(VS,gen);
		if not l then
			error "Second argument should lift";
		end if;
	end for;
	
	// convert voltage group to type GrpAb if needed
	if Type(VS`group) ne GrpAb then
		G, hom := AbelianGroup(CosetImage(VS`group,sub<VS`group|>));
		D := Domain(VS`assignment);

		newvoltages := [];
		for i in D do
			Append(~newvoltages,hom(VS`assignment(i)));
		end for;
			
		voltageassignment := map<D -> G|{<D[i], newvoltages[i]> : i in [1..#D]}>;	
	else
		G := VS`group;
		voltageassignment := VS`assignment;
	end if;
	
	X := VS`graph;
	V := VertexSet(X);
	n := #V;
	base := Omega[1];
	
	// construction of potentials and voltages of fundamental cycles together with their g images 
	
	Gr_voltOfCycles := [];
	Gr_potential := [];		
				
	for k in [1..#Generators(Gr)] do
		f := Gr.k;
		voltOfCycles := [];
		f_voltOfCycles := [];
		colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
		potential := [Id(G) : i in [1..n]];  // potential[i] is the voltage of the unique path on the tree from vertex V.i to vertex u
		fpotential := [Id(G) : i in [1..n]]; // f_potential[i] is the voltage of the unique path on the tree from vertex f(V.i) to vertex f(u)
		Q := [Index(base)];                            // vertices that has to be searched
		colour[Index(base)] := 1;                      // vertex u is coloured grey(1), since it has been searched at the moment
		
		while IsEmpty(Q) eq false do         // BFS 
			currentVertex := V.(Q[1]);
			currentfVertex := AutOnVertex(X,f,currentVertex);
			Remove(~Q,1);
			for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
				nextfVertex := AutOnVertex(X,f,nextVertex);
				if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
					multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex))); // function Indices(u,u) returns doubled indices?
					for dart in multi do                                                       // cotree edges
						if Index(nextVertex) lt Index(currentVertex) then
							voltage := -(potential[Index(nextVertex)]) + voltageassignment(dart) + potential[Index(currentVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(nextfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(currentfVertex)];
							else
								fvoltage := - fpotential[Index(nextfVertex)] + voltageassignment(dart^f) + fpotential[Index(currentfVertex)];
							end if;
						else
							voltage := - potential[Index(currentVertex)] + voltageassignment(dart) + potential[Index(nextVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(currentfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(nextfVertex)];
							else
								fvoltage :=  - fpotential[Index(currentfVertex)] + voltageassignment(dart^f) + fpotential[Index(nextfVertex)];
							end if;
						end if; 
						if G!voltage notin sub<G|voltOfCycles> then
							Append(~voltOfCycles,voltage);
							Append(~f_voltOfCycles,fvoltage);
						end if;
					end for;	
				elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
					multi := Indices(currentVertex, nextVertex);                                // all multiple edges
					dart := multi[1];                                                           // edge with the first index is tree edge
					colour[Index(nextVertex)] := 1; 
					Append(~Q,Index(nextVertex));
					Remove(~multi,1);
								
					// computing voltage and fvoltage
					if Index(nextVertex) lt Index(currentVertex) then
						voltage := voltageassignment(dart);
						if dart^f mod 2 eq 0 then
							fvoltage := - voltageassignment(dart^f-1);
						else
							fvoltage := voltageassignment(dart^f);
						end if;
					else
						voltage := - voltageassignment(dart);
						if dart^f mod 2 eq 0 then
							fvoltage := voltageassignment(dart^f-1);
						else
							fvoltage := - voltageassignment(dart^f);
						end if;
					end if;  
							
					potential[Index(nextVertex)] := voltage + potential[Index(currentVertex)];
					fpotential[Index(nextfVertex)] := fvoltage + fpotential[Index(currentfVertex)];
					for dart in multi do                                                       // cotree edges
						if Index(nextVertex) lt Index(currentVertex) then
							voltage := - potential[Index(nextVertex)] + voltageassignment(dart) + potential[Index(currentVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(nextfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(currentfVertex)];
							else
								fvoltage := - fpotential[Index(nextfVertex)] + voltageassignment(dart^f) + fpotential[Index(currentfVertex)];
							end if;
						else
							voltage := - potential[Index(currentVertex)] + voltageassignment(dart) + potential[Index(nextVertex)];
							if dart^f mod 2 eq 0 then
								fvoltage := - fpotential[Index(currentfVertex)] - voltageassignment(dart^f-1) + fpotential[Index(nextfVertex)];
							else
								fvoltage := - fpotential[Index(currentfVertex)] + voltageassignment(dart^f) + fpotential[Index(nextfVertex)];
							end if;
						end if; 
						if G!voltage notin sub<G|voltOfCycles> then
							Append(~voltOfCycles,voltage);
							Append(~f_voltOfCycles,fvoltage);
						end if;
					end for;
				end if;
			end for;
			colour[Index(currentVertex)] := 2;                                                    //  searching of the vertex is completed	
		end while;
		Append(~Gr_voltOfCycles,f_voltOfCycles);
		Append(~Gr_potential,fpotential);
	end for;

	// construct g# (as a matrix) for each generator in Gr
	MatRep := [];
	Ve := RSpace(Integers(),#Generators(G));
	C := RMatrixSpace(Integers(),#voltOfCycles,#Generators(G));
	M := RMatrixSpace(Integers(),#Generators(G),#Generators(G));

	// for each generator G.i find coefficients x_i such that G.i = x_1*voltOfCycles[1]+...+x_r*voltOfCycles[r] e.i. solve linear system  
	mat := [];
	for v in voltOfCycles do
		mat := mat cat Eltseq(v);
	end for;
	
	// rows of BMatrix generates relations of an abelian voltage group
	block := [];
	for i in [1..#Generators(G)] do
		row := Eltseq(G.i);
		row[i] := Integers()!Order(G.i);
		block := block cat row;
	end for;
	Bmatrix := M!block;
	
	// extend matrix by Bmatrix in order to solve system in abelian group
	sol := [];
	for i in  [1..#Generators(G)] do 
		e := [0 : i in [1..#Generators(G)]];
		e[i] := 1;
		Append(~sol, Solution(Transpose(HorizontalJoin(Transpose(C!mat),Bmatrix)),Ve!e));
	end for;
	
	for k in [1..#Generators(Gr)] do
		f_voltOfCycles := Gr_voltOfCycles[k];
		Im := [];
		//compute the image of each generator in G of f#
		for i in  [1..#Generators(G)] do 
			im := Id(G);
			for j in [1..#f_voltOfCycles] do
				im := im + sol[i][j]*f_voltOfCycles[j];
			end for;
		Im := Im cat Eltseq(im);	
		end for;
		Append(~MatRep, M!Im);
	end for;
	
	//initialization of variables
	Delta := [];
	IsDefined := [0 : i in [1..n]];
	IsDefined[Index(base)] := -1;
	numv := 0;
	
	for i in [1..#Generators(Gr)] do
		if IsDefined[Index(AutOnVertex(X,Gr.i,base))] eq 0 then // image of the base vertex of Gr.i is not defined yet
			numv := numv + 1;
			Append(~Delta,AutOnVertex(X,Gr.i,base));
			IsDefined[Index(AutOnVertex(X,Gr.i,base))] := numv;	
		end if;
	end for;
	
	MS := RMatrixSpace(Integers(),#Generators(G),#Generators(G)*numv); // matrices
	Tr := <<MS!0,Ve!0> : i in [1..n]>;
	
	for ver in Delta do
		for k in [1..#Generators(G)] do
			Tr[Index(ver)][1][k, (IsDefined[Index(ver)] - 1)*#Generators(G) + k] := 1;
		end for;
	end for;

	A := [];
	b := [];
	
	while IsEmpty(Delta) eq false do
		for i in [1..#Generators(Gr)] do
			ma := MatRep[i];
			mapped := AutOnVertex(X,Gr.i,Delta[1]);
			zetaQ := potential[Index(Delta[1])];
			zetagQ := Gr_potential[i][Index(mapped)];
			Lab := Tr[Index(Delta[1])][1];
			lab := Tr[Index(Delta[1])][2];
			I := Tr[Index(AutOnVertex(X,Gr.i,base))][1];
			
			//update matrix Lab!!!! it is actually just multiplication of two matrices
			for l in [1..#Generators(G)*numv] do
				column := Transpose(Lab)[l];
				LabT := Transpose(Lab);
				LabT[l] := (column*ma);
				Lab := Transpose(LabT);
			end for;
			
			Lab := Lab + I;
			// update vector lab
			lab := Ve!(lab*ma) + Ve!(Ve!Eltseq(zetaQ)*ma) - Ve!Eltseq(zetagQ);
			
			if IsDefined[Index(mapped)] eq 0 then 
				Append(~Delta,mapped);
				IsDefined[Index(mapped)] := 1;		
				Tr[Index(mapped)][1] := Lab;
				Tr[Index(mapped)][2] := lab;
			else
				A := A cat Eltseq(Lab - Tr[Index(mapped)][1]);	
				b := b cat Eltseq(Tr[Index(mapped)][2] - lab);
			end if;
		end for;
		Remove(~Delta,1);
	end while;

	for ver in Omega do
		if IsDefined[Index(ver)] eq 0 then
			Delta := [ver];
			while IsEmpty(Delta) eq false do
					for i in [1..#Generators(Gr)] do
						ma := MatRep[i];
						mapped := AutOnVertex(X,Gr.i,Delta[1]);
						zetaQ := potential[Index(Delta[1])];
						zetagQ := Gr_potential[i][Index(mapped)];
						Lab := Tr[Index(Delta[1])][1];
						lab := Tr[Index(Delta[1])][2];
						I := Tr[Index(AutOnVertex(X,Gr.i,base))][1];
						
						//update matrix Lab
						for l in [1..#Generators(G)*numv] do
							column := Transpose(Lab)[l];
							LabT := Transpose(Lab);
							LabT[l] := (column*ma);
							Lab := Transpose(LabT);
						end for;
						
						Lab := Lab + I;
						// update vector lab
						lab := Ve!(lab*ma) + Ve!(Ve!Eltseq(zetaQ)*ma) - Ve!Eltseq(zetagQ);
						
						if IsDefined[Index(mapped)] eq 0 then 
							Append(~Delta,mapped);
							IsDefined[Index(mapped)] := 1;		
							Tr[Index(mapped)][1] := Lab;
							Tr[Index(mapped)][2] := lab;
						else
							A := A cat Eltseq(Lab - Tr[Index(mapped)][1]);	
							b := b cat Eltseq(Tr[Index(mapped)][2] - lab);
						end if;
					end for;
					Remove(~Delta,1);
				end while;
		end if;
	end for;		
	
	rowA := Integers()!(#A/(#Generators(G)*numv));	
	E:= Bmatrix;
	for i in [2..Integers()!(rowA/#Generators(G))] do
		E := DiagonalJoin(E,Bmatrix);
	end for; 
		
	Vsp := RSpace(Integers(),rowA);
	Esp := RMatrixSpace(Integers(),rowA,#Generators(G)*numv);
	// solve the system
	w,e,t := IsConsistent(Transpose(HorizontalJoin(Esp!A,E)),Vsp!b);
	
	if w then
		return w,e,t,Esp!A,Vsp!b;
	else
		return w,_,_,Esp!A,Vsp!b;
	end if;
end function;



//========================================================================================
//                Construction of admissible solvable regular covers 
//========================================================================================


ReconstructCover := function(Xtilde, X, semi, p)
/* Given covering projection 'p'::Map mapping from covering graph 'X'::GrphMultUnd to base graph 'Y'::GrphMultUnd with semiedges 'semi'::Map
 * this function constructs a permutation voltage space on X.
 */

	freeIndex := [1 : i in [1..#VertexSet(X)]];
	labelofvertex := [0 : i in [1..#Vertices(Xtilde)]];
	
	for i in [1..#Vertices(Xtilde)] do
		proj := Index(CovProjOnVertex(Xtilde,X,p,VertexSet(Xtilde).i));
		labelofvertex[i] := freeIndex[proj];
		freeIndex[proj] := freeIndex[proj] + 1;
	end for;

	n := Integers()!(#VertexSet(Xtilde)/#VertexSet(X));
	VoltGroup := Sym(n);
	
	volt := [[0: j in [1..n]] : i in [1..#EdgeSet(X)]];
	
	for i in Indices(Xtilde) do
		projDart := p(i);
		ini := labelofvertex[Index(InitialVertex(EdgeSet(Xtilde).i))];
		ter := labelofvertex[Index(TerminalVertex(EdgeSet(Xtilde).i))];
		
		if projDart mod 2 eq 1 then
			volt[Integers()!((projDart + 1)/2)][ini] := ter;
		else
			volt[Integers()!((projDart + 1)/2)][ter] := ini;
		end if;
	end for;
	
	voltP := [];
	for i in volt do
		Append(~voltP, VoltGroup!i);
	end for;
	
	I := {@2*i-1 : i in [1..#EdgeSet(X)]@};
	darts := {@2*i-1 : i in [1..#EdgeSet(X)]@};
	VS := VoltageSpace("Permutation",VoltGroup,X,darts,[semi(i) : i in I],voltP);
		
	return VS;	
end function;

StandardBasisVector := function(d, i) 
//Returns a standard basis vector of dimension d with 1 at i'th place.
	vec := [];
	
	for j:=1 to d do
		if  not i eq j then
			vec := vec cat [0];
		else
			vec := vec cat [1];
		end if;
	end for;
	
	return vec;
end function;

DiagonalVector := function(d, p)
//Returns a vector of dimension d with all components equal to p.
 	
 	diag := [];
 	
 	for i in {1..d} do
 		diag := diag cat [p];
 	end for;
 	
 	return diag;
end function;

BFSCycles := function(X)
//Given a multigraph X this function returns BFS tree stored in sequences parentVertex, parentDart together with positive cotree darts CAp.
	V := VertexSet(X);
	n := #V;
	colour := [0 : i in [1..n]];         // all vertices are coloured white(0)
	parentVertex := [0 : i in [1..n]];        // BFS tree vertices
	parentDart := [0 : i in [1..n]];   // BFS tree(i-th position represents the vertex V.i and we store dart connecting V.i and its ancestor)
	u := V.1;                            // base vertex
	Q := [Index(u)];                     // vertices that has to be searched
	CAp := [];                           // odd darts of cotree edges e.g. positive cotree arcs

	colour[Index(u)] := 1;               // vertex u is coloured grey(1), since it has been searched at the moment
	
	while IsEmpty(Q) eq false do         // BFS 
		currentVertex := V.(Q[1]);
		Remove(~Q,1);
		for nextVertex in Neighbours(currentVertex) do                                     // search all neighbours of current vertex
			if colour[Index(nextVertex)] eq 1 then                                         // cotree edge
				multi := SetToSequence(SequenceToSet(Indices(currentVertex, nextVertex)));
				for dart in multi do
					Append(~CAp, dart);	
				end for;	
			elif colour[Index(nextVertex)] eq 0 then                                        // tree edge
				multi := Indices(currentVertex, nextVertex);                                // all multiple edges
				dart := multi[1];                                                           // edge with the first index is tree edge
				
				parentVertex[Index(nextVertex)] := Index(currentVertex); 
				if Index(nextVertex) le Index(currentVertex) then
					parentDart[Index(nextVertex)] := dart;
				else
					parentDart[Index(nextVertex)] := dart + 1;
				end if;
				colour[Index(nextVertex)] := 1; 
				Append(~Q,Index(nextVertex));
				Remove(~multi,1);
				
				for dart in multi do;                                                         // cotree edges
					Append(~CAp, dart);	
				end for;
			end if;
		end for;
		colour[Index(currentVertex)] := 2;                                                    //  searching of vertex is completed
	end while;

	return CAp, parentVertex, parentDart;
	
end function;


Hash := function(X, CAp, parentVertex, parentDart, generators, R)
/* Given a multigraph X and a sequence CAp of positive cotree darts together with BFS tree parentVertex, parentDart,
 * a set of generators (permutations g on darts in Aut X) and a ring R this function returns the matrices Transpose(g^#), relative to the basis 
 * of Hom_1(X;R) given by cotree darts in CAp for each g in generators
 */
 
	X := StandardGraph(X); 
	XV := Vertices(X);
	
	beta := #CAp;                         // betti number of a graph
	
	Alg := MatrixAlgebra(R, beta);
	N := RModule(Alg);
	
	// voltage assignment on darts of X
	I := {@2*i-1 : i in [1..#EdgeSet(X)]@};
	S := [N!0 : i in [1..#EdgeSet(X)]];
	for j in [1..#CAp] do
		k := Integers()!((CAp[j] + 1)/2);
		S[k] := N!StandardBasisVector(beta, j);
	end for;
	
	vol := map < I -> N | { < 2*i - 1, N!S[i] > : i in [1..#I]} >;                          //homological cover e.i. voltage assignment
		
	setGenMat := [];
	
	for g in generators do 
		//edges of mapped BFStree with permutation g
		EgX:={{Index(AutOnVertex(X, g, VertexSet(X).i)),Index(AutOnVertex(X, g, VertexSet(X).(parentVertex[i])))}: i in [2..#XV]};
		gX:=Graph<#XV|EgX>;
		gVX:=Vertices(gX);
	
		gparentDart := [0 : i in [1..#gVX]]; 
		for i in [2..#XV] do
			gparentDart[Index(AutOnVertex(X, g, VertexSet(X).i))] := parentDart[i]^g;
		end for;
		gparent:=[0 : i in [1..#gVX]];                               // BFS tree
		voltparent:=[N!0 : i in [1..#gVX]];
		Q:=[Index(AutOnVertex(X, g, VertexSet(X).1))];                                                    // rooted in mapped vertex 1^g 
		gparent[Index(AutOnVertex(X, g, VertexSet(X).1))]:=Index(AutOnVertex(X, g, VertexSet(X).1));
		
		while IsEmpty(Q) eq false do                                 // BFS search on mapped tree
			currentNode:=Q[1];
			Remove(~Q,1);
			for u in Neighbours(gVX!currentNode) do
				nextNode:=Index(u);
				if gparent[nextNode] eq 0 then
					gparent[nextNode]:= currentNode;
					dart := gparentDart[nextNode];
					if dart mod 2 eq 1  then 
						voltparent[nextNode]:= voltparent[currentNode]+ vol(dart); //potential direction towards base vertex
					else
						voltparent[nextNode]:= voltparent[currentNode] - vol(dart - 1);
					end if;
					Append(~Q,nextNode);
				end if;
			end for;
		end while;

		volFCg := [];
		for i in CAp do
			a := EdgeSet(X).i;
			s := N!0;
			beg := InitialVertex(a);
			ter := TerminalVertex(a);
			gbeg := Index(AutOnVertex(X, g, beg));
			gter := Index(AutOnVertex(X, g, ter));
			gdart := i^g;
			if 1 eq gdart mod 2 then
				s:= - voltparent[gbeg] + vol(gdart) + voltparent[gter];
			else
				s:= - voltparent[gbeg] - vol(gdart - 1) + voltparent[gter]; 
			end if;
			volFCg := volFCg cat [s];
		end for;

		MS:=[];
		for x in volFCg do 
			MS := MS cat [x[i] : i in [1..beta]]; 
		end for;
		setGenMat := Append(setGenMat,Transpose(Alg!MS));
	end for;
	
	return setGenMat;
	
end function;


InvariantMinimalSubspaces := function(GL , setGenMat, bound)
/* given a general linear group GL containing a subgroup H generated by matrices 
 * in the set setGenMat  returns the sequence of minimal H-invariant subspaces  
 * up to dimension bound (represented by its basis) 
 */
 
	H := sub< GL | setGenMat >;
	N := GModule(H);
	L := MinimalSubmodules(N);

	invH := <>;
		
	for M in L do 
		if Dimension(M) le bound then
			phi := Morphism(M,N);
			gen :=[ phi(x) : x in Basis(M)];
			invH := Append(invH,gen);
		end if; 
	end for;

	return invH;
end function;

HomCoordinates := function(d, p) 
//Returns homogeneous coordinates of space ZZ_p^d
  coord := <<1>>;
  for i in [1..(d-1)] do
 		coor1 := <>;
  		for co in coord do 
  			coor1 := Append(coor1,Append(co,0));
  		end for;
		C:= CartesianPower({0..(p-1)},i);
		for c in C do
			coor1 := Append(coor1,Append(c,1));
		end for;
		coord:=coor1;
	end for;
	return coord;
end function;

OneDimInvariantSubspaces := function(p ,setGenMat) 
/* given a set of matrices setGenMat over field of prime p that generate group H this function returns the sequence of one-dimensional H-invariant subspaces  
 * (represented by its basis) 
 */
	
	beta := #Eltseq(setGenMat[1][1]);
	Alg := MatrixAlgebra(GF(p), beta);
	
	numM := #setGenMat;
	all := Set(CartesianPower({0..(p-1)},numM));
	EigenVec :=<>;
	
	for comb in all do
		A := AddScaledMatrix(Alg!setGenMat[1],-comb[1],Id(Alg));
		for i in [2..numM] do
			B := AddScaledMatrix(Alg!setGenMat[i],-comb[i],Id(Alg));
			A := HorizontalJoin(A,B);
		end for;
		K := Basis(Kernel(A));
		if not IsEmpty(K) then
			if #K eq 1 then
				//print K;
				EigenVec := Append(EigenVec,[K[1]]);	
			else
				Coor := HomCoordinates(#K, p);
				for co in Coor do
					v := co[1]*K[1];
					for i in [2..#K] do
						v := v + co[i]*K[i];
					end for;
					EigenVec := Append(EigenVec,[v]);	
				end for;
			 end if; 
		end if;
	end for;
	
	return EigenVec;

end function;

SequenceElementsToIntegers := function(Seq)
//given a sequence of integers of some other type, returns the sequence of integer type assigned to elements

	seqInts := [];
	for i in {1..#Seq} do
		seqInts := seqInts cat [Integers()!Seq[i]];	
	end for;
	
	return seqInts;
end function;

Voltages := function(InvS, d) 
//given list of invariant subspaces InvS it returns a list of voltages where d is the betti number of a graph

	voltages := [];
	
	for inv in InvS do
		voltage := [];
		for i in {1..d} do
			vol := [];
			for g in inv do
				vol := vol cat [g[i]];
			end for;
			voltage := voltage cat [vol];
		end for;
		voltages := voltages cat [voltage];
	end for;

	return voltages;
end function;

ConstructCover := function(X, semi, CAp, parentVertex, parenDart, Inv, p, max, generators, basegraph, projection,all);

	graphGenPairs := [];	

	// Betti number
	betti := #CAp;
	
	XV := Vertices(X);
	
	// Make the voltage matrices from generators of invariant spaces
	voltMats := Voltages(Inv, betti);
	
	// Main loop
	for voltMat in voltMats do
		// Dimension of invariant subspace
		size := #voltMat[1];
		
		// We check if the lifted graph has at most max vertices (otherwise, we make no lift)
		if p^size * #XV le max then

			// Create a voltage group Z_p^size
			voltGroup := AbelianGroup(DiagonalVector(size, p));
						
			// voltage assignment on darts of X
			I := {@2*i-1 : i in [1..#EdgeSet(X)]@};
			S := [voltGroup!0 : i in [1..#EdgeSet(X)]];
			for j in [1..#CAp] do
				k := Integers()!((CAp[j] + 1)/2);
				S[k] := voltGroup!SequenceElementsToIntegers(voltMat[j]);
			end for;	
			
			VS := VoltageSpace("Cayley",voltGroup,X,I,[semi(i) : i in I],S);
			
			// Construction of the covering graph
			Xtilde, semi1, proj := DerivedCover(VS);
			
			VS := ReconstructCover(Xtilde,basegraph, semi, proj*projection);
			// check whether VS is equivalent to any of the previous covers
			isequiv := false;
			for rec in all do
				VS1 := rec`space;
				if IsEquivalent(VS,VS1) then
					isequiv := true;
					break;
				end if;
			end for;
			
			if not isequiv then
				// Lifting the generators
				// it is enough to lift only one generator of the group of covering transformation since we are interested in minimal covers
			
				id := Id(Sym(2*#EdgeSet(X)));
				u1 := VertexSet(Xtilde)!<VertexSet(X).1,1>;
				u2 := VertexSet(Xtilde)!<VertexSet(X).1,2>;
				lift := Lift(Xtilde,X,proj,id,u1,u2);		
				genSet1 := {lift};	
			
				genSet2 := {};
				for g in generators do
					u1 := VertexSet(Xtilde)!<VertexSet(X).1,1>;
					u2 := VertexSet(Xtilde)!<AutOnVertex(X,g,VertexSet(X).1),1>;
					lift := Lift(Xtilde,X,proj,g,u1,u2);
					genSet2 := genSet2 join {lift};	
				end for;			
				genSet := genSet1 join genSet2;
				// Write a list of records.
				// Cover Record Format: <graph, generators>
				CRF := recformat< graph, semi, generators, space, project>;
									
				record := rec<CRF | graph := Xtilde, semi:= semi1, generators := genSet, space := VS , project := proj*projection>;
				graphGenPairs := graphGenPairs cat [record];
			end if;
		end if;
	end for;
	
	return graphGenPairs;
end function;


FindCover := function(X, semi, basegraph, generators, setGenMat, prime, max, projection,all);
	
	// Find a set of cotree arcs and corresponding BFS tree
	CAp, parentVertex, parentDart := BFSCycles(X);
	
	// Betti number
	betti := #Edges(X) - (#Vertices(X) - 1);
	
	bound := Floor((Log(max/#Vertices(X))/(Log(prime))));

	GL := GeneralLinearGroup(betti, GF(prime));

	
	// Find minimal invariant subspaces
	if (bound eq 1) and (prime lt 4 ) then 
		Inv := OneDimInvariantSubspaces(prime ,setGenMat);
	elif (bound eq 1) and  (#setGenMat lt 3) and (prime lt 24)  then
		Inv := OneDimInvariantSubspaces(prime ,setGenMat);
	elif (bound eq 1) and  (#setGenMat lt 4) and (prime lt 12) then
		Inv := OneDimInvariantSubspaces(prime ,setGenMat);
	elif (bound eq 1) and  (#setGenMat lt 5) and (prime lt 6) then
		Inv := OneDimInvariantSubspaces(prime ,setGenMat);
	else
		Inv := InvariantMinimalSubspaces(GL, setGenMat, bound);
	end if;
	// Find covering graphs		

	return ConstructCover(X, semi, CAp, parentVertex, parentDart, Inv, prime, max, generators, basegraph, projection, all);
	
end function;


SolvableCovers := function(X, semi, genSet, max);
	/* Given a connected graph 'X'::GrphMultUnd with semiedges 'semi'::Map,  a set 'genSet'::SetEnum of automorphisms of X and 
	 * a natural number 'max'::RngIntElt this function constructs all permutation voltage assignments giving rise to pairwise nonequivalent
	 * <genSet>-admissible solvable regular covers of X up to order max of the respective covering graphs
	 */
	
	CAp, parentVertex, parentDart := BFSCycles(X);
	basegraph := X;
	// The set of all lifts
	allLifts := [];
	liftsOnLevel := [];
	maxPrime := Floor(max/#Vertices(X));
	
	// construct a matrix for each generator
	setGenMat := Hash(X, CAp ,parentVertex, parentDart, genSet, Integers());		
	
	for p in {2..maxPrime} do
		if(IsPrime(p)) then
			//initialize projection to be trivial
			I := [i : i in [1..2*#EdgeSet(basegraph)]];
			id := map<I -> I | i :-> i>;
			liftsOnOriginal := FindCover(X, semi, basegraph, genSet, setGenMat, p, max, id, liftsOnLevel);
			liftsOnLevel := liftsOnLevel cat liftsOnOriginal;
		end if;
	end for;

	allLifts := allLifts cat liftsOnLevel;
	
	// loop through all lifts in current level
	liftsOnNext := [];
	
	while #liftsOnLevel gt 0 do
		liftsNewNext := [];
		
		for rec in liftsOnLevel do
			tmp := #Vertices(rec`graph);
			maxPrime := Floor(max / tmp);
			
			if maxPrime gt 1 then
				CAp, parentVertex, parentDart := BFSCycles(rec`graph);
				setGenMat := Hash(rec`graph, CAp ,parentVertex, parentDart, rec`generators, Integers());	
				for p in {2..maxPrime} do
					if(IsPrime(p)) then
						liftsOnNext := FindCover(rec`graph, rec`semi, basegraph, rec`generators, setGenMat, p, max, rec`project, allLifts);
						allLifts := allLifts cat liftsOnNext;
						liftsNewNext := liftsNewNext cat liftsOnNext;
					end if;
				end for;
			end if;
		end for;		
		liftsOnLevel := liftsNewNext;
	end while;

	return allLifts;

end function;


