newPackage(
	"StatePolytope",
    	Version => "1.0", 
    	Date => "July 9, 2008",
    	Authors => {
	     {Name => "Dave Swinarski", Email => "swinarsk@math.columbia.edu"}
	     },
    	Headline => "for computing the state polytope of an ideal",
	-- DebuggingMode should be true while developing a package, 
	--   but false after it is done
    	DebuggingMode => true 
    	)

export { --testGfan --to be added
     --testPolymake  --to be added
     isGfanCompatible, --documented
     gfanRingInput, 
     gfanIdealInput, 
     initialIdeals, --documented
     hilbertPt, 
     statePolytopePoints, 
     fullStatePolytopePoints, 
     printHilbPt, 
     createPolymakeInputFile, 
     statePolytopeInDegreeM, --documented
     maxUGBDegree, 
     fullStatePolytope,  --documented
     testMHilbertStability --documented
     }




isGfanCompatible = method(TypicalValue => Boolean)
isGfanCompatible(Ideal) := (I) -> (R:= ring(I);
       if isPolynomialRing(R) == false then (print "Not an ideal in a polynomial ring.  As of June 2008, gfan only supports ideals in polynomial rings over QQ or ZZ/p with p < 32749.", return false); 
       --I don't think large characteristics are supported in m2 yet, but in case they get added to m2 before gfan, let's check...
       if char(R) >= 32749 then (print "Characteristic is too large.  As of June 2008, gfan only supports ideals in polynomial rings over QQ or ZZ/p with p < 32749.", return false);
     if char(R) > 0 and coefficientRing(R) =!=  ZZ/char(R) then (print "As of June 2008, gfan only supports ideals in polynomial rings over QQ or ZZ/p with p < 32749.", return false);
     if char(R) == 0 and coefficientRing(R) =!=  QQ  then (print "As of June 2008, gfan only supports polynomial rings over QQ or ZZ/p with p < 32749.", return false);
     gensStringList := apply(numgens R, i -> toString (gens(R))_i);
     P := tally flatten flatten apply(numgens(R)-1, i-> apply(numgens(R)-i-1, j-> {match(gensStringList_i,gensStringList_(i+j+1)),match(gensStringList_(i+j+1),gensStringList_i)}));
     if P_true > 0 then (print "Don't use variable names in gfan such that one variable's name is a substring of another's.", return false);
     return true
     )

TEST ///
R = QQ[x,y];
I = ideal(x^2,y);
assert(isGfanCompatible I == true)
R = QQ[a,b]/(ideal(a^2));
I = ideal(b);
assert(isGfanCompatible I == false)
R = ZZ[a,b];
I = ideal(a^2,b);
assert(isGfanCompatible I == false)
R = QQ[a1,a11];
I = ideal(a1^2,a11);
assert(isGfanCompatible I == false)
///




--The next two functions print out the ring and the ideal in the format required by gfan

gfanRingInput = (I) -> (R:= ring(I);
     if isGfanCompatible(I) =!= true then error toString(isGfanCompatible(I));
     str:="";
     if char(R) == 0 then str = "Q[" else str = concatenate(concatenate("Z/",toString(char(R))) ,"Z[");
for i from 0 to (numgens(R)-2) do str = concatenate(str,concatenate(toString R_i, "," ));
return concatenate(str,concatenate(toString R_(numgens(R)-1), "]" ));
)


gfanIdealInput = (I) -> (toString apply((numgens(I)), k -> I_k)
    )



initialIdeals = method(
     TypicalValue => List
     )
initialIdeals(Ideal,String,RR) := (I,gfanCommand,gfanVersion) -> (
if isGfanCompatible(I) =!= true then error toString(isGfanCompatible(I));
gfanoutputfile := openInOut gfanCommand;
if gfanVersion >= 0.3 then gfanoutputfile << gfanRingInput(I) << gfanIdealInput(I);
if gfanVersion < 0.3 then (print "The user is responsible for making sure the ideal and its ring are compatible with this version of gfan", gfanoutputfile << gfanIdealInput(I));
gfanoutputfile << closeOut;
gfanoutputstring := get gfanoutputfile;
markedInitialIdealsString:=replace("[Q,Z].*]","",gfanoutputstring);
use(ring(I));
initialIdealsList = value replace("[+-][^,{}]*}","}",replace("[+-][^,{}]*,",",",markedInitialIdealsString));
return initialIdealsList
)


--How can I write a test for this, since the user needs to put in the gfan commeand and version?

hilbertPt = (n,m,L,monomialsList) -> (
     totalDegree := binomial( n-1 + m, n-1) * m / n ;
     totalDegreeVector := apply(n, k-> totalDegree);
     complementHilbPt := flatten exponents product flatten apply(#monomialsList, i -> {if product apply(#L, j -> monomialsList_i % L_j) != 0 then monomialsList_i else 1});
     idealHilbertPt := totalDegreeVector - complementHilbPt;
     return idealHilbertPt
     )

statePolytopePoints = method(
     TypicalValue => List)
statePolytopePoints(Ideal,List,ZZ) := (I,initialIdealsList,m) -> (
     n := numgens ring(I);
     monomialsList := flatten entries(basis(m,ring(I)));
     return apply(#initialIdealsList, k -> hilbertPt(n,m,initialIdealsList_k,monomialsList))
	  )

fullStatePolytopePoints = method(
     TypicalValue => List)
fullStatePolytopePoints(Ideal,List,ZZ) := (I,initialIdealsList,m) -> (
     sum apply(m+1,k-> statePolytopePoints(I,initialIdealsList,k))
     )

printHilbPt = (L) -> ( str:= "1";
for j from 0 to (numgens ring(I)-1) do str = concatenate(str,concatenate(" ",toString lift(L_j,ZZ)));
return str
)

createPolymakeInputFile = (statePolytopePointsList) -> (
openOut "temporarypolymakefile.txt";
polymakeinput =  concatenate("POINTS",newline);
for i from 0 to (#statePolytopePointsList - 1) do polymakeinput = concatenate(polymakeinput,concatenate(printHilbPt(statePolytopePointsList_i),newline));
"temporarypolymakefile.txt" << polymakeinput << closeOut
)

statePolytopeInDegreeM = method(
     TypicalValue => String)
statePolytopeInDegreeM(Ideal,ZZ,String,RR) := (I,m,gfanCommand,gfanVersion) -> ( initialIdealsList := initialIdeals(I,gfanCommand,gfanVersion);
createPolymakeInputFile(statePolytopePoints(I,initialIdealsList,m));
polymakesession = "!polymake temporarypolymakefile.txt VERTICES";
polymakesession << closeOut;
return get polymakesession
)

maxUGBDegree = (L) -> (
     (max apply(#L, i -> max apply(#(L_i), j -> degree L_i_j)))_0
	  )
     
fullStatePolytope = method(
     TypicalValue => String)
fullStatePolytope(Ideal,String,RR) := (I,gfanCommand,gfanVersion) -> ( initialIdealsList := initialIdeals(I,gfanCommand,gfanVersion);
createPolymakeInputFile(fullStatePolytopePoints(I,initialIdealsList,maxUGBDegree(initialIdealsList)));
polymakesession = "!polymake temporarypolymakefile.txt VERTICES";
polymakesession << closeOut;
return get polymakesession
)


testMHilbertStability = method(
     TypicalValue => String)
testMHilbertStability(Ideal,ZZ,String,RR) := (I,m,gfanCommand,gfanVersion) -> (   
     st := statePolytopeInDegreeM(I,m,gfanCommand,gfanVersion);
     str:= st;
     barycenter := toString ((sum(statePolytopePoints(I,initialIdeals(I,gfanCommand,gfanVersion),m))_0) / (numgens ring(I) ));
     barycenterstring:= concatenate("POINTS",concatenate(newline,"1 "));
     for i from 0 to (numgens(ring(I))-1) do barycenterstring = concatenate(barycenterstring, concatenate(barycenter," "));
     str = replace("VERTICES",barycenterstring,str);
     openOut "augmentedtemporarypolymakefile.txt";
     "augmentedtemporarypolymakefile.txt" << str << closeOut;
     polymakesession2 = "!polymake augmentedtemporarypolymakefile.txt VERTICES";
     polymakesession2 << closeOut;
     augmentedstatepolytope := get polymakesession2;
     if augmentedstatepolytope == st then "This ideal is GIT stable or semistable with respect to the maximal torus action" else "This ideal is GIT unstable with respect to the maximal torus action"
          )
     
     



beginDocumentation()
document { 
     
	Key => StatePolytope,
	Headline => "computes state polytopes of ideals",
	EM "StatePolytope", " computes state polytopes of ideals using ", TT "gfan", ", ", TT "M2", ", ", "and ", TT "polymake", ".  Specifically, it computes ", ITALIC "State", SUB "m", "(I)", " or ", ITALIC "State", "(I)", " as defined in Sturmfels's book ", ITALIC "Grobner bases and convex polytopes", ", page 14.  We assume that the user has recent versions of ", TT "gfan", " and ", TT "polymake", " on his or her system in addition to ", TT "M2", "."
	}



document {
	Key => {isGfanCompatible, (isGfanCompatible,Ideal)},
	Headline => "checks whether an ideal is supported by gfan 0.3",
	Usage => "isGfanCompatible I",
	Inputs => {
		"I" => Ideal => {}
		},
	Outputs => {
		Boolean => {}
		},
	TT "gfan", " version 0.3 supports ideals in polynomial rings over ", TT "QQ", " or ", TT "ZZ/p", " with p < 32749.  Also, variables should not be named in a way such that one variable is a substring of another (e.g. x1 and x11).  Given an ideal I, this function tests ring(I) for these properties.",  
	   
	" This function is provided by the package ", TO StatePolytope, ".",
	EXAMPLE lines ///
		R = QQ[x,y];
		I = ideal(x^2,y);
		isGfanCompatible I
		S = ZZ[a,b];
		J = ideal(a^2,b);
		isGfanCompatible J
		///,
		}
	
   

   
   
   document {
	Key => {initialIdeals, (initialIdeals,Ideal,String,RR)},
	Headline => "calls gfan and returns the list of initial ideals",
	Usage =>  "see example",
	Inputs => {
		Ideal => "the ideal",
		String => "the command for running gfan on your system, preceded by an exclamation point and enclosed in quotation marks",
		RR => "the version of gfan you are running"
		},
	Outputs => {
		{"the list of initial ideals of I.  Note that each initial ideal is given as a list of monomial generators in ring(I), not as a monomial ideal."
		}},
	"Given an ideal ", TT "I", " in a ring ", TT "R", " in ", TT "M2", ", this function calls ", TT "gfan", " and returns the list of initial ideals of ", TT "I", ". ",
	   
	"This function is provided by the package ", TO StatePolytope, ".",
	EXAMPLE lines ///
	        R = QQ[a,b];
		I = ideal(a^2+b^2,a*b);
		initialIdeals(I,"!gfan",0.3)  
		///,
		}
	
   

   
   document {
	Key => {statePolytopeInDegreeM, (statePolytopeInDegreeM,Ideal,ZZ,String,RR)},
	Headline => ("computes the mth state polytope of an ideal"),
	Usage => "see example",
	Inputs => {
		Ideal => "the ideal",
		ZZ => "enter m to get the mth state polytope",
		String => "the command for running gfan on your system, preceded by an exclamation point and enclosed in quotation marks",
		RR => "the version of gfan you are running"
		},
	Outputs => {
		{"the list of vertices of the mth state polytope,  in polymake format."}
		},
	
	"See Sturmfels's book ", ITALIC "Grobner bases and convex polytopes", ", page 14 for the definition of " , ITALIC "State", SUB "m", "(I). ",
	   
	 "Note: currently it is necessary that the program ", TT " polymake", " is actually run by the command ", TT "polymake", " on your system.  Also, a file called temporarypolymakefile.txt is created in the process.",
	  
	   
	" This function is provided by the package ", TO StatePolytope, ".",
	EXAMPLE lines ///
	        R = QQ[a..d];
		I = ideal(a*c-b^2,a*d-b*c,b*d-c^2);
		statePolytopeInDegreeM(I,3,"!gfan",0.3)  
		///,
		}
	
   
   
   document {
	Key => {fullStatePolytope, (fullStatePolytope,Ideal,String,RR)},
	Headline => "computes State(I)",
	Usage =>  "see example",
	Inputs => {
		Ideal => "the ideal",
		String => "the command for running gfan on your system, preceded by an exclamation point and enclosed in quotation marks",
		RR => "the version of gfan you are running"
		},
	Outputs => {
		{"the list of vertices of the state polytope, in polymake format."}
		},
	"Given an ideal I, this function calls ", TT "gfan", " to get the list of initial ideals of I.  It then determines the largest degree d of any generator in any initial ideal, computes the state polytopes ", ITALIC "State", SUB "m", "(I)", " for m from 1 to d, and sums these, yielding ", ITALIC "State", "(I). ",

" Note: currently it is necessary that the program ", TT" polymake", " is actually run by the command ", TT "polymake", " on your system.  Also, a file called temporarypolymakefile.txt is created in the process.",

	   
	" This function is provided by the package ", TO StatePolytope, ".",
	EXAMPLE lines ///
	        R = QQ[a..d];
		I = ideal(a*c-b^2,a*d-b*c,b*d-c^2);
		fullStatePolytope(I,"!gfan",0.3)  
		///,
		}
	
	
	 document {
	Key => {testMHilbertStability, (testMHilbertStability,Ideal,ZZ,String,RR)},
	Headline => "determines whether the mth Hilbert point of I is GIT stable with respect to the action of the maximal torus",
	Usage =>  "see example",
	Inputs => {
		Ideal => "the ideal",
		ZZ => "specifies which Hilbert point to test",
		String => "the command for running gfan on your system, preceded by an exclamation point and enclosed in quotation marks",
		RR => "the version of gfan you are running"
		},
	Outputs => {
		{"a message saying (semi)stable or unstable"}
		},
	"Bayer and Morrison showed that GIT stability of the mth Hilbert point of I with respect to the maximal torus acting on a polynomial ring by scaling the variables can be tested by whether ", ITALIC "State", SUB "m", "(I) contains a certain point.",
	
	" Note: currently it is necessary that the program ", TT" polymake", " is actually run by the command ", TT "polymake", " on your system.",
	   
	" This function is provided by the package ", TO StatePolytope, ".",
	EXAMPLE lines ///
	        R = QQ[a..d];
		I = ideal(a*c-b^2,a*d-b*c,b*d-c^2);
		testMHilbertStability(I,3,"!gfan",0.3)
		I = ideal(a^2,b^2,b*c);
		testMHilbertStability(I,3,"!gfan",0.3) 
		///,
		}
   
   
   end
   
   
installPackage "StatePolytope"
installPackage("StatePolytope", RemakeAllDocumentation=>true)
check StatePolytope