//////////////////////////////////////////////////////////////////////
// Computation of a Stratification  of finite Group actions
// Preliminary version (2.11.2001)
// Implementation by : Thomas Bayer, Institut fuer Informatik, 
//                     Technische Universitaet Muenchen
// www:    http://wwwmayr.informatik.tu-muenchen.de/personen/bayert/
// email : bayert@in.tum.de
//

LIB "rinvar.lib";
LIB "finvar.lib";
LIB "primdec.lib";

//////////////////////////////////////////////////////////////////////
/*

PROCEDURES:
DecomposeHilbertMap(ideal phi)
FindStrata(list matrices)
GetGroupEquations(list groupElements)
HilbertMap(list matrices)
Stratify(ideal G, ideal action, ideal phi)

SUBPROCEDURES:
AllIntersections(list subspaces)
CollectSubspaces(ideal G, ideal action, list vs)
ComputeStrata(ideal G, ideal action, list vs)
GenericOrbitLength(ideal G, ideal action, ideal V)
GroupEquations(list G)
IdealContainedQ(ideal I, ideal J)
Orbit(ideal G, ideal action, ideal I)
ReduceIdealList(list data)


set printlevel=1 for a trace of 'FindStrata' and printlevel=2 for 
additional output in 'CollectSubspaces'
HINT: 	For the first application of 'FindStrata' one should set
	printlevel=2 (invariant_ring might take long)

*/


//////////////////////////////////////////////////////////////////////


proc DecomposeHilbertMap(ideal phi)
"USAGE:  DecomposeHilbertMap(phi); ideal phi
PUROPSE: compute all intersections of the irreducible subspaces of the
         singular locus of phi
ASSUME:  each entry of 'subspaces' is of the form list(ideal, ideal)
RETURN:  list of ideals, each ideal defines an intersection	       
EXAMPLE: example
"
{
	def DHMR = basering;
	string ringSTR = "ring DHMS = 0, (" + varstr(basering) + "), lp;";
	execute(ringSTR);
	ideal I, J, phi;
	list primaryDec;

	phi = imap(DHMR, phi);
	I = std(minor(jacob(phi), nvars(basering)));
	J = radical(I);
	primaryDec = primdecGTZ(J);
	list allSubspaces = ReduceIdealList(AllIntersections(primaryDec));
	dbprint(dbPrt,"Singular locus generated by " + string(size(J)) + " elements");
	dbprint(dbPrt, string(size(primaryDec)) + " irreducible components - compute all intersections");

	setring DHMR;
	return(imap(DHMS, allSubspaces)));
}

proc FindStrata(list matrices)
"USAGE:   FindStrata(matrices); list matrices
PUROPSE: compute all intersections of the ideals of the list 'subspaces'
ASSUME:  each entry of 'subspaces' is of the form list(ideal, ideal)
RETURN:  list: 	_[1] = list of strata
		_[2] = ring, containing G, action, eqns of strata, phi
		       
EXAMPLE: example
"
{
	def FSR1 = basering;
	ideal phi;
	int i, j, n, nrStrata;
	int ti1, ti2, ti3, ti4, ti5, ti6;
	list grpList, allSubspaces;
	string ringSTR1, ringSTR2;

	int dbPrt = printlevel-voice+2;

	// compute the Hilbert map and all group elements

	ti1 = timer;
	grpList = HilbertMap(matrices);
	phi = grpList[1];
	ti1 = timer - ti1;

	dbprint(dbPrt,"Hilbert map has " + string(size(phi)) + " coordinates");

	// compute the singular locus of phi
	// compute all intersections of the irreducible components
	// obtained in the previous step 

	ti2 = timer;
	allSubspaces = DecomposeHilbertMap(phi);	
	ti2 = timer - ti2;
	dbprint(dbPrt, string(size(allSubspaces)) + " different subspaces");

	// define the action, compute equations for the group
	
	setring FSR1;
	
	def FSR2 = GetGroupEquations(grpList[2]);
	setring FSR2;
	ideal phi = imap(FSR1, phi);
	list grpList, allSubspaces, collectedSpaces, collectedStrata, orbitLength, result;

	// collect subspaces according to their generic orbit length
	
	ti4 = timer;
	allSubspaces = imap(FSR1, allSubspaces);
	result = CollectSubspaces(G, action, allSubspaces);
	ti4 = timer - ti4;
	dbprint(dbPrt, "List of all orbit lengths: [" + string(result[1]) + "]");

	// compute equations of the Semistrata

	ti5 = timer;
	orbitLength = result[1];
	collectedSpaces = result[2];
	nrStrata = 0;
	for(i = 1; i <= size(orbitLength); i++) {
		collectedStrata[i] = ComputeStrata(G, action, collectedSpaces[i]);
		nrStrata = nrStrata + size(collectedStrata[i]);
	}
	ti5 = timer - ti5;
	dbprint(dbPrt, string(nrStrata) + " different strata ");
	dbprint(dbPrt, "-- Semistrata ------- ");
	for(i = 1; i <= size(orbitLength); i++) {
		dbprint(dbPrt, " - Generic orbit length = " + string(orbitLength[i]));
		for(j = 1; j <= size(collectedStrata[i]); j++){
			print(string(j) + " : " + string(collectedStrata[i][j]));
		}
	}

	// prepare result

	dbprint(dbPrt, "-- Timing -------------");
	dbprint(dbPrt, " hilbert map                     : " + string(ti1));
	dbprint(dbPrt, " singular locus and decomposition: " + string(ti2));
	dbprint(dbPrt, " compute all subspaces           : " + string(ti3));
	dbprint(dbPrt, " collect subspaces               : " + string(ti4));
	dbprint(dbPrt, " collect strata                  : " + string(ti5));

	export(phi);
	export(collectedStrata);

	setring FSR1;
	return(list(imap(FSR2, collectedStrata), orbitLength, FSR2));
}

proc GetGroupEquations(list groupElements)
"USAGE:   GetGroupEquations(groupElements)
PUROPSE: computes the action and the defining ideal of the group
         given by the list 'groupElements'
ASSUME:  
RETURN:  ring:	ideal G = defining ideal of the group
		ideal action = ideal defining the action
EXAMPLE: example
"
{
	int i, j, n;
	string ringSTR;

	def GGER = basering;
	n = ncols(groupElements[1]);
	ringSTR = "ring GGES = 0, (s(1.." + string(n^2) + "), ";
	ringSTR = ringSTR + varstr(basering) + "), dp;";

	execute(ringSTR);
	ideal action, G;
	list grpList;
	poly p;

	for(i = 1; i <= n; i++){
		p = 0;
		for(j = 1; j <= n; j++){
			p = p + s(j + (i - 1)*n ) * var(n^2 + j);
		}
		action[i] = p;
	}
	//grpList = imap(GGER, groupElements);
	G = GroupEquations(imap(GGER, groupElements));
	export(G);
	export(action);
	return(GGES);
}

proc HilbertMap(list matrices)
"USAGE:   HilbertMap(matrices); list matrices
PUROPSE: compute all intersections of the ideals of the list 'subspaces'
ASSUME:  each entry of 'subspaces' is of the form list(ideal, ideal)
RETURN:  list: 	_[1] = ideal defining phi
		_[2] = list of all elements of the group generated by
		       'matrices'
EXAMPLE: example
"
{
	ideal phi, inv, invBase;
	list grp_reynolds, elementList;
 
	int dbPrt = printlevel-voice+2;

	dbprint(dbPrt, " computing the invariant ring");
	matrix PI, SI, FI = invariant_ring(matrices);
	inv = PI, FI;
	invBase = mstd(inv)[2];		// minimal generating set

	grp_reynolds = group_reynolds(matrices);// compute all group elements
	elementList = grp_reynolds[2..size(grp_reynolds)];
	phi = evaluate_reynolds(grp_reynolds[1], invBase);
	return(list(phi, elementList));
}

proc Stratify(ideal G, ideal action, ideal phi)
"USAGE:   Stratify(G,action, phi); ideal G, action, phi
PUROPSE: compute the stratification of the representation space of G 
ASSUME:  phi = Hilbert map (projection of K^n onto the quotient K^n/G
         phi = generated by fundamental invariants
RETURN:  basering = K[s,t]
NOTE:    
EXAMPLE: example
"
{
	int i, n;
	ideal sing;
	list maxSubspaces, allSubspaces, orbitLength;
	list collectedSpaces, collectedStrata, result; 

	int dbPrt = printlevel-voice+2;

	n = ncols(action);
	sing = radical(std(minor(jacob(phi), n)));
	print(" sing = " + string(sing));
	maxSubspaces = primdecGTZ(sing);
	allSubspaces = ReduceIdealList(AllIntersections(maxSubspaces));

	print(" ----------------------------- ");
	print(" Number of subspaces = " + string(size(allSubspaces)));

	// collect subspaces having the same orbit length

	result = CollectSubspaces(G, action, allSubspaces);
	orbitLength = result[1];
	collectedSpaces = result[2];

print("--------------------------------");
	print(" orbitlengths = " + string(orbitLength));
	print(" collectedSpaces"); 
	for(i = 1; i <= size(collectedSpaces); i++) {
		print(string(i) + " : " + string(collectedSpaces[i]));
	}

	// compute orbits of subspaces, each orbit corresponds 
	// to an orbit type

	for(i = 1; i <= size(orbitLength); i++) {
		collectedStrata[i] = ComputeStrata(G, action, collectedSpaces[i]);
	}
	return(list(orbitLength, collectedStrata));	
}

//////////////////////////////////////////////////////////////////////


proc AllIntersections(list subspaces)
"USAGE:   AllIntersections(subspaces); list subspaces
PUROPSE: compute all intersections of the ideals of the list 'subspaces'
ASSUME:  each entry of 'subspaces' is of the form list(ideal, ideal)
RETURN:  
EXAMPLE: example
"
{
        int i, j, k, l, m, n, maxLen;
        ideal I, J;  
        list elem, all, allId, currentElems, newElems, currentIdeals, newIdeals, index;

        n = size(subspaces);
        for(i = 1; i <= n; i++) { 
                all[i] = list(i); 
                allId[i] = subspaces[i][2];
        }

        k = 0;
        currentElems = all;
        currentIdeals = allId;

        for(l = 2; l <= n; l++){
                for(i = 1; i <= size(currentElems); i++) {
                        elem = currentElems[i]; 
                        I = currentIdeals[i];
                        m = size(elem);
                        index = list();
                        for(j = elem[m] + 1; j <= n; j++) { index[j - elem[m]] = j;}
                        
                        for(j = 1; j <= size(index); j++) {
                                k = k + 1;
                                newElems[k] = elem + list(index[j]);
                                J = I, subspaces[index[j]][2];
                                newIdeals[k] = std(J);
                        }
                }
                currentElems = newElems;
                currentIdeals = newIdeals;
                all = all + newElems;
                allId = allId + newIdeals;
                newElems = list();
                newIdeals = list();
                k = 0;
        }
        return(allId);
}   

proc CollectSubspaces(ideal G, ideal action, list vs)
"USAGE:   CollectSubspases(G,action,vs); list vs
PUROPSE: collect invariant subspaces according to their generic orbit lenght 
ASSUME:  basering = K[s,t]
RETURN:  
NOTE:    
EXAMPLE: example
"
{
	int i, j, count, index, glen;
	list collectedSpaces, orbitLength, result;
	int dbPrt = printlevel-voice  + 2;

	count = 0;
	for(i = 1; i <= size(vs); i++) {
		dbprint(dbPrt, "  CollectSubspaces:" + string(i) + " of " + string(size(vs)));
		glen = GenericOrbitLength(G,action, vs[i]);
		index = FirstEntry(orbitLength, glen);
		if(index == 0) {
			count++;
			orbitLength[count] = glen;
			collectedSpaces[count] = list(vs[i]);
		}
		else {
			collectedSpaces[index] = collectedSpaces[index] + list(vs[i]);
		}
	}
	return(list(orbitLength, collectedSpaces));
}

proc ComputeStrata(ideal G, ideal action, list vs)
"USAGE:   ComputeStrata(G,action, I); list vs
PUROPSE: compute the strata out of the vectorspaces contained in vs.  
ASSUME:  basering = K[s,t], all vectorspaces of 'vs' have the same 
         generic orbit length.
RETURN:  ideal in K[t], no change of basering 
"
{
	int i, j, k, len;
	ideal I;
	list strata, subspaces;
	int dbPrt = printlevel-voice;

	i = 0; 
	k = 0;
	subspaces = vs;
	len = size(subspaces);

	while(i < len) {
		dbprint(dbPrt, "  ComputeStrata, subspace #:" + string(i));
		i++;
		if(subspaces[i] != 0) {
			I = std(Orbit(G, action, subspaces[i]));
			k++;
			strata[k] = I;
			for(j = i + 1; j <= len; j++) {
				if(IdealContainedQ(subspaces[j], I)) {	
					subspaces[j] = ideal(0);
				}
			}
		}
	}  
	return(strata);
}

proc GenericOrbitLength(ideal G, ideal action, ideal V)
"USAGE:   GenericOrbitLength(G,action,V); ideal V
PUROPSE: compute the length of the orbits of points from a dense set of V
ASSUME:  basering = K[s,t]
RETURN:  
NOTE:    computes the image of the generic point of the vectorspace V
EXAMPLE: example
"
{
        int i, k, n, N, c;
	ideal zero, tempV;
        string parSTR, varSTR, ringSTR1, ringSTR2;

 	n = size(action);
        N = nvars(basering);
	for(i = 1; i <= N - n; i++) {zero[i] = var(i);}
	tempV = V,zero;
	k = dim(std(tempV));
	// if(k == 0) { error(" variety is zerodimensional ");}
       
        parSTR = "(0, a(1.." + string(k) + "))";
	varSTR = "(" + varstr(basering) + ", Y(1.." + string(n) + "))";
        ringSTR1 = "ring RGO = " + parSTR + ", " + varSTR;
	ringSTR1 = ringSTR1 + ", (dp(" + string(N) + "), dp(" + string(n) + "));";
	ringSTR2 = "ring SGO = " + parSTR + ", (Y(1.." + string(n) + ")), dp;";	

	def BGO = basering;
	execute(ringSTR1);
	ideal action, G, I, image, V;

	G = imap(BGO, G);	
	action = imap(BGO, action);
	V = imap(BGO, V);

	// compute the correct equations !

	I = std(V);
	c = 0;
	for(i = N - n + 1; i <= N; i++) {
		if(deg(reduce(var(i), I)) > 0) {
			c++;
			I = I, var(i) - a(c); I = std(I);
		}
	}
	for(i = 1; i <= n; i++) {
		I = I, action[i] - var(N + i);
	}
	I = I,G; 
	image = nselect(std(I),1,N);
	execute(ringSTR2);
	ideal image = std(imap(RGO, image));
	return(vdim(image)));
}

proc GroupEquations(list G)
"USAGE:   GroupEquations(G); list G
PUROPSE: compute the ideal of the finite group G
ASSUME:  G is a list of matrices
RETURN:  ideal
EXAMPLE: example 
"
{
        int i, j, k, n; 
        ideal I, J, action;

        n = ncols(G[1]);        // number of variables for G

        // M[i][j] -> var((i - 1) * n + j)

        for(i = 1; i <= n; i++) {
                action[i] = 0;
                for(j = 1; j <= n; j++) {
                        I = I, var((i - 1) * n + j) - G[1][i, j];
                        action[i] = action[i] + var((i - 1) * n + j) * var(n^2 + j);
                }
        }
        for(k = 2; k <= size(G); k++) {
                J = 0;
                for(i = 1; i <= n; i++) {
                        for(j = 1; j <= n; j++) {
                                J = J, var((i - 1) * n + j) - G[k][i, j];
                        }
                }
                I = std(intersect(I, J));
        }
        return(I);
}


proc IdealContainedQ(ideal I, ideal J)
// is J contained in I
{
	int i, containedQ;
	ideal I1 = std(I);

	i = 0; 
	containedQ = 1;

	while(i < ncols(J) && containedQ == 1) {
		i++;
		if(reduce(J[i], I1) != 0) { containedQ = 0;}
	}
	return(containedQ);
}

proc Orbit(ideal G, ideal action, ideal I)
"USAGE:   Orbit(G,action, I); list vs
PUROPSE: compute the orbit of V(I) w.r.t. the action of G 
ASSUME:  basering = K[s,t]
RETURN:  ideal in K[t], no change of basering 
NOTE:    
EXAMPLE: example
"
{
	int i, n, N;
	ideal mapid, X;

	n = ncols(action);
	N = nvars(basering);
	X = G,I;
	def ROB = basering;
	def ROR = ImageVariety(X, action);

	for(i = 1; i <= n; i++) {mapid[i] = var(N - n + i);}
	map F = ROR,mapid;
	return(F(imageid));
}

proc ReduceIdealList(list data)
"USAGE:   ReduceIdealList(dc); list dc
PUROPSE: eliminates double entries
ASSUME:  
RETURN:  
EXAMPLE: example
"
{
        int i, j, k, len, loop;
        ideal I, J;
        list idList, allSubspaces;

	allSubspaces = data;
        len = 0;

        for(i = 1; i <= size(data); i++) {
                I = std(allSubspaces[i]);
		if(size(I) > 0) {
			len++;
			idList[len] = I;
                	
			for(j = i + 1; j <= size(data); j++) {
				J = std(allSubspaces[j]);
				if(IdealContainedQ(I, J) && IdealContainedQ(J, I)) {
					allSubspaces[j] = 0;
				}
				else {	allSubspaces[j] = J; }
			}
                }
        }
        return(idList);
}

proc FirstEntry(data, elem)
// Type : int
// Purpose : position of first entry equal to elem in data (from left to right)
{
	int i, pos;
 
	i = 0;
	pos = 0;
	while(!pos && i < size(data)) {
		i = i + 1;
		if(data[i] == elem) { pos = i;}
	}
	return(pos);
}