//////////////////////////////////////////////////////////////////////
// Computation of the Stratification  of the Symmetry Group of 
// Nickel-Titan alloys (22.10.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

<"strata.sing";		// load package

///////////////////////////////////////////////////////
//define basering and matrices to compute the invariant ring

print("compute invariants");

ring R = 0, (t(1..6)), dp(6);

matrix M[6][6]=0,1,0,0,0,0, 1,0,0,0,0,0, 0,0,1,0,0,0, 0,0,0,-1,0,0, 0,0,0,0,0
,-1, 0,0,0,0,1,0;
matrix N[6][6]=0,0,1,0,0,0, 0,1,0,0,0,0, 1,0,0,0,0,0, 0,0,0,0,-1,0, 0,0,0,-1,
0,0, 0,0,0,0,0,1;
matrix HH(1..2);

int ti; ti = timer; 
//matrix PI,SI,FI = invariant_ring(M,N);  // benoetig ca. 1h

ideal I = t(1)+t(2)+t(3),t(4)^2+t(5)^2+t(6)^2,t(1)^2+1/2*t(1)*t(2)+t(2)^2+1/2*t(1)*t(3)+1/2*t(2)*t(3)+t(3)^2,t(4)*t(5)*t(6),t(1)*t(4)^2+t(2)*t(4)^2-1/2*t(3)*t(4)^2-1/2*t(1)*t(5)^2+t(2)*t(5)^2+t(3)*t(5)^2+t(1)*t(6)^2-1/2*t(2)*t(6)^2+t(3)*t(6)^2,t(1)^3+t(2)^3+t(3)^3,t(4)^4+1/2*t(4)^2*t(5)^2+t(5)^4+1/2*t(4)^2*t(6)^2+1/2*t(5)^2*t(6)^2+t(6)^4,t(1)^2*t(4)^2+t(2)^2*t(4)^2+t(1)*t(3)*t(4)^2+t(2)*t(3)*t(4)^2+t(1)*t(2)*t(5)^2+t(2)^2*t(5)^2+t(1)*t(3)*t(5)^2+t(3)^2*t(5)^2+t(1)^2*t(6)^2+t(1)*t(2)*t(6)^2+t(2)*t(3)*t(6)^2+t(3)^2*t(6)^2,t(3)*t(4)^4+t(1)*t(5)^4+t(2)*t(6)^4;

list tt = group_reynolds(M,N);		// compute all group elements
list Grp = tt[2..size(tt)];
// ideal I = PI,FI;

ideal J = mstd(I)[2];			// construct a minimal set of fundamental invariants
ideal phi = evaluate_reynolds(tt[1],J);
ti = timer - ti; ti;

print(" invariants after " + string(ti));

///////////////////////////////////////////////////////
// compute the singular locus of the Hilbert map phi
// together with an irredundant decomposition

print("compute decomposition of phi");

ti = 0;
ring R1 = 0,(t(1..6)), lp;

ideal phi = imap(R, phi);
ti = timer;

ideal I = std(minor(jacob(phi),6));
ideal J = radical(I);
list v2 = primdecGTZ(J);

ti = timer - ti; ti;

print(" decomp after " + string(ti));

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

print("compute all intersections");

ti = 0;
ti = timer; list allSubspaces = ReduceIdealList(AllIntersections(v2));ti = timer - ti; ti;
print(" all subspaces after " + string(ti));

///////////////////////////////////////////////////////
// compute the generic orbit length of the subspaces 
// obtained in the previous step.

print("compute generic orbit length");

ring S1 = 0,(s(1..36),t(1..6)), dp;

ideal action;
action[1] = s(1)*t(1) + s(2)*t(2) + s(3)*t(3) + s(4)*t(4) + s(5)*t(5) + s(6)*t(6);
action[2] = s(7)*t(1) + s(8)*t(2) + s(9)*t(3) + s(10)*t(4) + s(11)*t(5) + s(12)*t(6); 
action[3] = s(13)*t(1) + s(14)*t(2) + s(15)*t(3) + s(16)*t(4) + s(17)*t(5) + s(18)*t(6);
action[4] = s(19)*t(1) + s(20)*t(2) + s(21)*t(3) + s(22)*t(4) + s(23)*t(5) + s(24)*t(6);
action[5] = s(25)*t(1) + s(26)*t(2) + s(27)*t(3) + s(28)*t(4) + s(29)*t(5) + s(30)*t(6);
action[6] = s(31)*t(1) + s(32)*t(2) + s(33)*t(3) + s(34)*t(4) + s(35)*t(5) + s(36)*t(6);

ideal phi = imap(R, phi);
list allSubspaces = imap(R1, allSubspaces);
ideal zero;
for(int i = 1; i <= 36; i++) {zero[i] = var(i);}

int ti3 = timer;
ideal G = GroupEquations(imap(R,Grp));
ti3 = timer - ti3;

///////////////////////////////////////////////////////
// collect the subspaces according to their orbit length
// and compute their orbits to obtain the semistrata

print("collect subspaces according to their orbit length");
print("and compute their orbits");

int ti4 = timer;
list result = CollectSubspaces(G, action, allSubspaces);
ti4 = timer - ti4;

int ti5 = timer;
list orbitLength = result[1];
list collectedSpaces = result[2];
list collectedStrata;
for(int i = 1; i <= size(orbitLength); i++) {
	collectedStrata[i] = ComputeStrata(G, action, collectedSpaces[i]);
}
ti5 = timer - ti5; ti5;
print(" strata after " + string(ti));

///////////////////////////////////////////////////////
// collectedStrata[3][1] and collectedStrata[5][1] are
// the semistrata corresponding to the low/high energy
// phase. Their images in the orbit space are computed

int ti6 = timer;
ideal X1 = collectedStrata[5][1];
ideal X2 = collectedStrata[3][1];
X1 = X1,G;
X2 = X2,G;

def S1a = ImageVariety(X1, phi);
def S2a = ImageVariety(X2, phi);
ti6 = timer - ti6; ti6;
print(" images after " + string(ti6));


//////////////////////////////////////////////////////////////
/*
ring S2 = 0,(t(1..6), Y(1..9)), (dp(6),dp(9));

ideal phi = imap(R,phi);

ideal I;
for(int i = 1; i <= size(phi); i++) { I[i] = phi[i] - Y(i);} 
//option(prot);
intvec w = 1,1,1,1,1,1;
intvec hw = HilbertWeights(phi,w);
intvec hs = HilbertSeries(I,hw);
intvec hs1 = HilbertSeries1(hw);

int ti = timer; ideal J = std(I,hs,hw); ti = timer - ti; ti;
ideal J1 = mstd(nselect(J,1,6))[2];
*/
print(" orbit space after " + string(ti));
print(" groupequations " + string(ti3));
print(" genericorbit " + string(ti4));
print(" strata " + string(ti5));