--Package to compute some objects from the project of Milena Hering
-- and Diane Maclagan

--Code to compute whether overline{C(lambda)} cap C(mu) neq emptyset
--Question from Evain's paper 
--"Incidence relations among the Schubert cells
--of equivariant punctual Hilbert schemes "

--Todo:
--play with try and alarm to interrupt computation

newPackage(
     "TEdges",
     Version => "0.5",
     Date => "21 November, 2009",
     Authors => {{
	       Name => "Diane Maclagan",
	       Email => "D.Maclagan@warwick.ac.uk",
	       HomePage => "http://www.warwick.ac.uk/staff/D.Maclagan"}},
     Headline => "T-Edges in the Hilbert scheme of points",
     DebuggingMode => true
     )

export { 
     sigArrows,
     Clambda,
--     edgeIdeal,  --Deleted for now, as it has a bug.
     listOfMonomialIdeals,
     ComputeGraph,
     writeGraph,
     arrowMaps,
     findEdgeDirection,
--     findGradedEdgeDirection,
     posetLess,
     findTedge
    }

---------------------------------------------------------------------------
-- CODE
---------------------------------------------------------------------------



--Take monomial ideal in two generators (assumed x,y and return a list of 
-- generators in order x^a, x^by^c,....,y^d

sortLex=M->(
	R:=ring M;
	if numgens R >2 then error "this is just for Q[x,y]";
	exps:=reverse sort apply(flatten entries gens M,m->flatten exponents m);
	return apply(exps,s->(R_0^(s_0)*R_1^(s_1)));
);
		
--Take a pair of monomial ideals I, J, and return true if I < J in the poset
-- order.  The number bound is the bound on how high to check the poset order
--We assume that the ring of I and J is positively graded.

posetLess=method();

posetLess(Ideal, Ideal, ZZ) :=  (I,J,bound) ->(
     local basI; local basJ;
     R:=ring I;
     less:=true;
     i:=1;
     while less and i<bound+1 do (
	  basI=flatten entries gens image basis(i,I);
	  basJ=flatten entries gens image basis(i,J);
          if (#basI)!=(#basJ) then less=false;
     	  if #basI>0 and less then (
	       scan(#basI,j->(if basI#j>basJ#j then less=false;));
	  );
	  i=i+1;
     );
     return(less);   
);     
     
--This one assumes that the ideals are artinian     
-- actually assuming two variables right now.
posetLess(Ideal, Ideal, List, MutableHashTable) := (I,J,v,Ibasis) ->(
     local d; local Id; local Jd; local yy;
     R:=ring I;
     Imonos:=new MutableHashTable;
     Jmonos:=new MutableHashTable;
     a:=-v#1; b:=v#0;
     if v#0<0 then (a=-a; b=-b;);
     yy:=symbol yy;
     S:=QQ[ yy_1, yy_2,Degrees=>{a,b}];
     f:=map(S,R,{yy_1,yy_2});
     Is:=f I;
     Js:=f J;
     Ibas:=apply(Ibasis#I,m->(f(m)));
     Jbas:=apply(Ibasis#J,m->(f(m)));
     scan(Ibas,m->(
	       d=(degree(m))_0;
	       if Imonos#?d then Imonos#d=append(Imonos#d,m)
	       	    else Imonos#d = {m};
     ));      
     scan(Jbas,m->(
	       d=(degree(m))_0;
	       if Jmonos#?d then Jmonos#d=append(Jmonos#d,m)
	       	    else Jmonos#d = {m};
     ));      
     less := true;
     keysI:=sort keys Imonos;
     keysJ:=sort keys Jmonos;
     if not(keysI==keysJ) then less = false;
     i:=0;
     while less and i<#keysI do (
--     	  if keysI#i!=keysJ#i then less = false
--	  else (
	       Id=sort(Imonos#(keysI#i));
       	       Jd=sort(Jmonos#(keysI#i));
--<<Id<<endl<<Jd<<endl<<endl;	       
     	       if #Id!=#Jd then less=false else 
       	           scan(#Id,j->(if Id#j<Jd#j then less=false;));
	       i=i+1;
--	  );
     );
     return less;	                           
);

--Code to work out the significant arrows of a monomial ideal 
--Input: A monomial ideal (in type Ideal), string "+", "-", "all"
--     	   Assumption:  the ideal lives in a polynomial
--     	    	 ring in two variables [x,y] with degrees of 
--     	    	 x and y given (which implies
--     	    	 that the term order is lex x>y)
--Output: A MutableHashTable with keys the generators of the ideal and
--     entries the list of monomials corresponding to heads of the
--     (positive, negative, all) significant arrows starting at that
--     generator.

--"Feature"  Can get nonemptylength empty lists { , , }

sigArrows=(I,dir)->(
  local Spair;
  R:=ring I;    
  SigArrows:=new MutableHashTable;
  --This next command should have the generators in decreasing
  --order (so x^* comes first...)
  gensI:=sortLex I;
  for i from 0 to #gensI-1 do (
     m:= gensI#i;       	    
     Basm:=flatten entries basis(degree(m),R/I);
     SigArrows#m={};
     SigArrows#m = for std in Basm list (
	  std=lift(std,R);
	  if m<std then (
	     if (dir=="+" or dir=="all") then (
	     	  SPair:=lift(std*lcm(gensI#(i-1),m)/m,R);
		  if (SPair % I)==0 then std
	     )
	  ) 
     	  else (
	       if (dir=="-" or dir=="all") then (
     	    	  SPair=lift(std*lcm(gensI#(i+1),m)/m,R);
 		  if (SPair % I)==0 then std
	       )	  
     	  )
     );	       
     SigArrows#m=delete(null,SigArrows#m);
  );
  return(SigArrows);
);     

--Then build the Evain equations for C(lambda)
--Input: Monomial ideal in k[x,y,Degrees=>{{a},{b}}], string "+", or  "-".
--Output: Ideal in variables x,y, z_1...z_s being the
-- 	     universal family over C(M) = { I in H_ab^H : in_{y>x}(I) = M
--     	    	      	   
Clambda=(M,sign)->(
     local varsS;
     R:=ring M;
     a:=degree(R_0);
     b:=degree(R_1); 
     --Need to check that this will always be in the right order
     GensM:=sortLex M;
     if sign=="-" then GensM=reverse GensM;
     Sig:=sigArrows(M,sign);
     Igens:=new MutableHashTable;
     m:=GensM#0;
     N:= #(flatten for m in keys Sig list Sig#m);
--Warning!  Seems to be assuming variables are x and y - document!! 
    if sign=="+" then (
    	 varsS=gens R;
    )
    else varsS=reverse gens R;
    z:=symbol z;
    S:=(coefficientRing R)[varsS, z_1..z_N,MonomialOrder=>Eliminate 2];
     Igens#m=substitute(m,S);
     alphacount:=0;
     for i from 1 to #GensM-1 do (
	  m=GensM#i;
	  mprev:=GensM#(i-1);
	  f:=lift ((substitute(m,S)*Igens#mprev)/substitute(mprev,S) ,S);
	  xdiff:=numerator(mprev/m);
	  if #Sig#m>0 then 
     	  for alpha in Sig#m do (
	       	  j:=position(GensM, m->((alpha*xdiff % m) ==0));
		  msyz:=GensM#j;
     	       	  newterm:=S_(alphacount+2)*substitute(alpha,S)*(Igens#(GensM#j))/substitute(GensM#j,S);
		  newterm=lift(newterm,S);
		  f=f+newterm;
		  alphacount=alphacount+1;		  
     	  );
     	  Igens#m=f;
      );
     I:=ideal(apply(keys Igens,m->Igens#m));
     return(I);  	  
);   	     


--Compute equations for the edge scheme
--Input: Monomial ideals I,J in k[x,y,Degrees=>{a,b}], with I<J in the poset order
--???get poset order correct?
edgeIdeal=(I,J)->(
  R:=ring I;
--  if not posetLess(I,J) then error "I should be less than J";
  Iideal:=Clambda(I,"-");
  Jideal:=Clambda(J,"+");
  NI:=numgens ring Iideal;
  NJ:=numgens ring Jideal;
  N:=NI+NJ-4;
  z:=symbol z;
  S:=coefficientRing R[gens R,symbol z_1..symbol z_N,MonomialOrder=>Lex];
  fI:=map(S,ring Iideal,apply(NI,i->S_i));
  IS:=fI Iideal;
  --Might be bug in the next line
  Jvars:=flatten append({S_1,S_0},apply(NJ-2,i->(S_(NI+i))));
  fJ:=map(S,ring Jideal,Jvars);
  JS:=fJ Jideal;
  GI:=forceGB gens IS;
  reducedJ:=flatten entries ((gens JS) % GI);
--<<reducedJ<<endl;  
  coeffs:=flatten apply(reducedJ, f->(
	    flatten entries((coefficients(f,Variables=>{S_0,S_1}))_1)
  ));	    
--  coeffs = ideal mingens ideal coeffs;
  return {IS, JS, coeffs};
);     


---------------------------------------------------------------
-- Produces a list of all monomial ideals with d standard 
-- monomials in n variables.
---------------------------------------------------------------
listOfMonomialIdeals = method();

listOfMonomialIdeals (ZZ,ZZ,Ring) := (d,n,R) -> (
     return(listOfMonomialIdealsIF(ideal(gens R),d,{0_R}));
);     

listOfMonomialIdeals (ZZ,ZZ) := (d,n) -> (
     x:=symbol x;
     R:=QQ[x_1..x_n];
     return(listOfMonomialIdeals(d,n,R))
);

--Produces all monomial ideals with the same Hilbert function as I
listOfMonomialIdeals (Ideal) := (I) -> (
     R:= ring I;
     --Assumes Artinian
     d:=rank source basis(R/I);
--     J:= ideal leadTerm I;
--     Jgens := sortLex J;
     Bound := 1;
     for i from 0 to numgens R -1 do 
     	  Bound = Bound*d*(degree(R_i))_0;
--     Bound := (degree(Jgens_0))_0 * (degree(Jgens_(-1)))_0;
     return(listOfMonomialIdeals(I,Bound));
--    h:=new MutableHashTable;
--     for i from 0 to Bound do
--     	 h#i=hilbertFunction(i,R/J);
--     return(listOfMonomialIdealHilbFcnH(ideal gens R, h, {0_R}, Bound));
);

listOfMonomialIdeals (Ideal, ZZ) := (I,Bound) -> (
     h:=new MutableHashTable;
     R:=ring I;
     for i from 0 to Bound do
     	 h#i=hilbertFunction(i,R/I);
     return(listOfMonomialIdealHilbFcnH(ideal gens R, h, {0_R}, Bound));
);
     
---------------------------------------------------------------
-- Produces a list of monomial ideals with a given number of 
-- standard monomials contained in  a given ideal
-- Input: ideal I that the monomial ideals should be contained in, d=#standard monomials
--     	    F=list of monomials representing boxes not to remove
-- Output:  List L of monomial ideals
---------------------------------------------------------------
listOfMonomialIdealsIF = (I,d,F) -> (
     R := ring I;
     r := #(flatten entries basis(R/I));
     if r > d then error "input ideal is too small"
     else if r == d then {I}
     else (
	  L := {};
	  G := flatten entries gens I;
	  G = select(G, g -> not member(g,F));
	  local J;
	  if #G == 1 then (
	       J = monomialIdeal mingens(ideal(F) + sum(gens R, r -> ideal(r*G#0)));
	       L = L | listOfMonomialIdealsIF(J,d,F);
	       F = F | {G#0})
	  else scan(G, g -> (
		    J := monomialIdeal mingens(ideal delete(g,G) + ideal(F) + sum(gens R, r -> ideal(r*g)));
	       	    L = L | listOfMonomialIdealsIF(J,d,F);
	       	    F = F | {g}));
	  L));


---------------------------------------------------------------
--Produce a list of monomial ideals containing a given ideal with a fixed Hilbert function h
--Input:  Ideal I that the monomial ideals should be contained in,
--     	    mutable Hash table h consisting of the Hilbert function
--     	    F=list of monomials representing boxes not to remove
--     	    Bound = upper bound on how far to check Hilbert function
--     	    Important Note:  
--Assumptions:  The ring R of I has two variables, and is graded 
--Output:
---------------------------------------------------------------
listOfMonomialIdealHilbFcnH = (I, h, F,Bound) -> (
--<<gens I<<endl;     
     local diffH;
     R := ring I;
--<<R<<endl;
     --First check that Hilbert function is ok
     hI:=new MutableHashTable;
     scan(Bound+1,i->(hI#i=hilbertFunction(i,R/I)));
--scan(Bound+1,i->(hI#i=hilbertFunction(i,R/I); <<hI#i<<"  ")); <<endl;
     ItooSmall:=false;
     HtheSame:=true;
     scan(Bound+1, i->(diffH=hI#i-h#i; if diffH>0 then ItooSmall=true; if diffH<0 then HtheSame=false));
     if ItooSmall then error "Input ideal is too small"
   --If Hilbert function same return I
     else if HtheSame then (
	  gensI:=monomialIdeal select(flatten entries gens I, m->(
		    	 (degree(m))_0<Bound+1));
--	  <<gens gensI<<endl;
--<<endl;scan(Bound+1,i-><<hI#i<<"   "); <<endl;	  <<endl;
	  return({gensI});
	  )
     else (
	--if Hilbert function of R/I smaller create smaller ideal.
 	L:={};
     	G := flatten entries gens I;
     	G = select(G, g -> not member(g,F));
	  if (#G == 1 and (degree(G#0))_0<Bound+1) then (
	       if h#((degree(G#0))_0) > hI#((degree(G#0))_0) then (
	       	    J:= monomialIdeal mingens(ideal(F) + sum(gens R, r -> ideal(r*G#0)));
	       	    newL:= listOfMonomialIdealHilbFcnH(J,h,F,Bound);
		    if #newL>0 then L= L | newL;
	       	    F = F | {G#0}
		    );
	       )
	  else scan(G, g -> (
		if (degree(g))_0<Bound+1 then (    
    	       if h#((degree(g))_0) > hI#((degree(g))_0)then (
		    J := monomialIdeal mingens(ideal(delete(g,G)) + ideal(F) + sum(gens R, r -> ideal(r*g)));
	       	    newL:= listOfMonomialIdealHilbFcnH(J,h,F,Bound);
		    if #newL>0 then L=L | newL;
	       	    F = F | {g}
		    );
	       );
	       ));
	  );        
     	  return(L);     
);


--takes monomial ideals I and J and returns the the unique v
--for which both I and J are Z^n/v graded.
--Assumes that I and J are both artinian

findEdgeDirection=(II,JJ,Ibasis,Graded)->(
       	local Ibas; local Jbas; local Ivect; local Jvect; local vectIJ, local vIJ;
	R:=ring II;
	if not(Graded) then (
	  Ibas=Ibasis#II;
     	  Jbas=Ibasis#JJ;
	  Ivect=apply(gens R,i->0);
	  Jvect=apply(gens R,i->0);
	  scan(Ibas,m->(m=lift(m,R); Ivect=Ivect+((exponents(m))_0)));
	  scan(Jbas,m->(m=lift(m,R); Jvect=Jvect+((exponents(m))_0)));
	  vectIJ=Ivect-Jvect;
	  if gcd(vectIJ)>0 then vectIJ=vectIJ/gcd(vectIJ);
	  vIJ=apply(vectIJ,i->lift(i,ZZ));
	)
      	else (
	n:=numgens R;
	v0:=(exponents(1_R))_0;
	vv:=v0;
	NI:=max(flatten(apply(flatten entries gens II,m->degree(m))));
	NJ:=max(flatten(apply(flatten entries gens JJ,m->degree(m))));
	N:=max(NI,NJ);
	Direction:=true;
	scan(N+n,d->(
		Ibas=flatten entries basis(d,R/II);
		Jbas=flatten entries basis(d,R/JJ);
		Ivect=apply(gens R,i->0);
		Jvect=apply(gens R,i->0);
		scan(Ibas,m->(m=lift(m,R); Ivect=Ivect+((exponents(m))_0)));
		scan(Jbas,m->(m=lift(m,R); Jvect=Jvect+((exponents(m))_0)));
		vIJ=Ivect-Jvect;
		if gcd(vIJ)>0 then vIJ=vIJ/gcd(vIJ);
		if vv==v0 then vv=vIJ;
		if vv!=v0 and vv!=vIJ then Direction=false;
	));
      	if not(Direction) then vIJ={};
      );
    return(vIJ);
);



-- --takes monomial ideals I and J and returns the the unique v
-- --for which both I and J are Z^n/v graded.
-- --This is the non-artinian, but graded, case
-- findGradedEdgeDirection=(I,J)->(
--      	local Ibas; local Jbas; local Ivect; local Jvect; local vectIJ;
-- 	R:=ring I;
-- 	n:=numgens R;
-- 	v0:=(exponents(1_R))_0;
-- 	vv:=v0;
-- 	NI:=max(flatten(apply(flatten entries gens I,m->degree(m))));
-- 	NJ:=max(flatten(apply(flatten entries gens J,m->degree(m))));
-- 	N:=max(NI,NJ);
-- 	Direction:=true;
-- --???WORK OUT HOW HIGH TO GO???
-- 	scan(N+n,d->(
-- 		Ibas=flatten entries basis(d,R/I);
-- 		Jbas=flatten entries basis(d,R/J);
-- 		Ivect=apply(gens R,i->0);
-- 		Jvect=apply(gens R,i->0);
-- 		scan(Ibas,m->(m=lift(m,R); Ivect=Ivect+((exponents(m))_0)));
-- 		scan(Jbas,m->(m=lift(m,R); Jvect=Jvect+((exponents(m))_0)));
-- 		vectIJ=Ivect-Jvect;
-- 		if gcd(vectIJ)>0 then vectIJ=vectIJ/gcd(vectIJ);
-- 		if vv==v0 then vv=vectIJ;
-- 		if vv!=v0 and vv!=vectIJ then Direction=false;
-- 	));
-- 	if Direction then return(vectIJ)
-- 	else return({});
-- 	--changed from else return
-- );



--Find the equations for ideals on the T-edge joining I and J
--We'll assume I and J are ideals (not monomial ideals)
--Input: Ideals I, J, true/false (true if I, J are homogenous, false if they
-- 		are artinian)
--Output: ???
findTedge=(I,J,Ibasis,Graded,vect)->(
     	local f; local monom; local coeffs;
	R:=ring I;
	K:=coefficientRing R;
	n:=numgens R;
	--Get edge direction
	if vect=={} then (
	   if Graded then vect=findEdgeDirection(I,J,{},true)
		else vect=findEdgeDirection(I,J,Ibasis,false);
	);
	uvect:=unique vect;
	numvars:=0;
	if #uvect>1 or uvect_0!=0 then (
	  posm:=1_R;   
 	  negm:=1_R;
 	  scan(numgens R, i->(if vect_i>0 then posm=posm*(R_i)^(vect_i); 
		if vect_i<0 then negm=negm*(R_i)^(-vect_i);
		));
	  if posm!=1_R and negm!=1_R then (
	    --Work out how many variables we need
	    scan(flatten entries mingens I, mono->(
		while not(mono//negm == 0_R) do (
			mono=(mono //negm)*posm;
			if not(mono % I == 0_R) then 
				numvars=numvars+1;
		);
	    ));
	    scan(flatten entries mingens J, mono->(
		while not(mono//posm == 0_R) do (
			mono=(mono //posm)*negm;
			if not(mono % J == 0_R) then 
				numvars=numvars+1;
		);
	    ));
	  );
	);
	if numvars==0 then (
		return(ideal(1_R));
	);
	--Create ring and ideals.
	if numvars>0 then (
	yy:=symbol yy;
	aa:=symbol aa;     
	S:=K[yy_1..yy_n,aa_1..aa_(numvars),MonomialOrder=>Lex];
	sofar:=0;
	RinS:=apply(numgens R,i->S_i);
	phi:=map(S,R,RinS);
	posm=phi posm;
	negm=phi negm;
	I2:=ideal(0_S);
	scan(flatten entries mingens I, mono->(
		mono=phi mono;
		f=mono;
		while not(mono//negm==0_S) do (
			mono=(mono//negm)*posm;
			if not(mono % (phi I) == 0_S) then (
				f=f+S_(sofar+n)*mono;
				sofar=sofar+1;
			);
		);
		I2=I2+ideal(f);
	));
	J2:=ideal(0_S);
	scan(flatten entries mingens J, mono->(
		mono=phi mono;
		f=mono;
		while not(mono//posm==0_S) do (
			mono=(mono//posm)*negm;
			if not(mono % (phi J) == 0_S) then (
				f=f+S_(sofar+n)*mono;
				sofar=sofar+1;
			);		
		);
		J2=J2+ideal(f);
	));
      	I2=ideal mingens I2;
	J2=ideal mingens J2;
	--Switch them around if necessary so that the initial ideal
	--of the edge ideal in the standard lex order is I
	pos1:=min(positions((exponents(posm))_0, i-> i>0));
	neg1:=min(positions((exponents(negm))_0, i-> i>0));
	if pos1<neg1 then (
		Ktemp:=I2;
		I2=J2;
		J2=Ktemp;
		Ktemp=I;
		I=J;
		J=Ktemp;
	);
	--Now work out the equations on the edge ideals.
	T:=K[aa_1..aa_(numvars)];
	Mideal:=ideal(0_T);
	--First reduce I2 by itself
	M:=flatten entries gens gb I2;
	N:=select(M,m->((m % phi(I))==m));
	if N!={} then (
		scan(N,m->(
			{monom,coeffs}:=coefficients(m,Variables=>RinS);
			Mideal=Mideal+ideal substitute(coeffs,T);
		));
	);
	--Then reduce J2 by I2
	gb I2;
	M=flatten entries ((mingens J2) % I2);
	scan(M,m->(
		{monom,coeffs}=coefficients(m,Variables=>RinS);
		Mideal=Mideal+ideal substitute(coeffs,T);
	));
	Mideal=ideal mingens Mideal;
	return({vect,I2,J2,Mideal});
	);
);


---------------------------------------------------------------
--Main procedure - computes the graph of T-edges in the Hilbert
--scheme of d points in A^n.
--Warning - probably n=2 here again!
---------------------------------------------------------------
ComputeGraph=method();

--Input: d=number of points, n=dimension, K=field
--Output: graph on all monomial ideals with d standard monomials
-- 	     in n variables.

ComputeGraph (ZZ,ZZ) := (d,n)->(
	return(ComputeGraph(d,n,QQ));
);	

ComputeGraph (ZZ,ZZ,Ring) := (d,n,K)->(
        x:=symbol x;
	R:=K[x_1..x_n];
	L:=listOfMonomialIdeals(d,n);
	return(ComputeGraph(L,false,false));
);	

--Input:  Ideal I.     
--Output:  graph on all monomial ideals with the same Hilbert
--     	    function as I.


ComputeGraph (Ideal, Boolean) := (I,graded) ->(
     L:=listOfMonomialIdeals(I);
     return(ComputeGraph(L,false,graded));
);  
     
ComputeGraph (List, Boolean) := (L,graded)->(
     return(ComputeGraph(L,false,graded));
);     
     
--Input the list of monomial ideals we want the graph on	
--Slow is a boolean that says the code should write the output as it goes
-- which is useful for large examples
ComputeGraph (List, Boolean, Boolean) := (L,Slow, graded )->(
     	local I; local J; local K; local II; local JJ; local R;
	local Edge; local vectpos;
	R=ring L#0;
     	Ibasis:=new MutableHashTable;
	if not(graded) then (
	     scan(L, I->(Ibasis#I=apply(flatten entries basis(R/I), m->lift(m,R))));
	)
   	else (
	     degreesGens:=sort unique flatten apply(L, I->(apply(flatten entries gens I,m->((degree(m))_0))));
	     scan(L, I->(Ibasis#I=flatten(apply(degreesGens,d->(apply(flatten entries basis(d,R/I),m->lift(m,R)))))));
	);
	if Slow then (
	      fp:=openOut("EdgeList.out");
	      fp<<"Edges = { ";
	)
   	else EdgeLists:=new MutableHashTable;
	existEdge:={};
	scan(#L,i->(
	   I= L#i;
	   scan(#L-i-1,j->(
     	       if Slow then <<i<<"   "<<i+j+1<<endl;		     
		J=L#(i+j+1);
		vect:=findEdgeDirection(I,J,Ibasis,false);
--		vect:=findGradedEdgeDirection(I,J);		
--This is meant to check for poset-less - but hasn't been double-checked!!
--????FIX HERE!!!!
--     	        if vect#0>0 then (
--     	       	    II=J;
--		    JJ=I;
--		   )
     	        if not(unique(vect) == {0}) then (
       	        vectpos=0;
		while (vect#vectpos)==0 do vectpos=vectpos+1;
		if vect#vectpos<0 then (
    	       	    II=J;
		    JJ=I;
		)
	   	else (II=I; JJ=J);
     	        if not(graded) then Edge=findTedge(II,JJ,Ibasis,false,vect)
     	       	else Edge=findTedge(II,JJ,Ibasis,true,vect);		
--<<Edge<<endl;		
		if Slow then (
		     fp<<" { {" <<i<<", "<<i+j+1<<"}, "<<Edge << " }, "<<endl;
     	       	)		     
		else EdgeLists#{i,i+j+1}=Edge;
		if #Edge>3 and dim(Edge_3)>-1 then existEdge = append(existEdge,{i,i+j+1});
--		EdgeLists#{i,i+j+1}={vect} |Edge;
     	        )
	   	else EdgeLists#{i,i+j+1}=ideal(1_R);
	   ));
	));
    	if Slow then (
	     <<"There are "<<#existEdge<<" edges."<<endl;
	     close(fp);
	     return(existEdge);
        )	    
	else return({L,EdgeLists});
);



---------------------------------------------------------------
--Write GraphViz output for the graph of all monomial ideals
---------------------------------------------------------------

writeGraph = method();
writeGraph (ZZ, ZZ) := (d,n) ->(
     return(writeGraph(d,n, QQ, false));
);

writeGraph( ZZ, ZZ, Ring, Boolean) := (d,n,K,verbose) ->(
     L:=ComputeGraph(d,n,K);
     return(writeGraph(L,verbose,false));
);     
    

--The list has the form {list of ideals, Edges},
--where Edges is a Hash table.
--This is the format output by ComputeGraph
--verbose is a boolean indicating how much output we want on the screen.
writeGraph (List, Boolean,Boolean) := (LL, verbose,graded) -> (
     	L:=LL_0;
	Edges:=LL_1;
	R:=ring L#0;
	--In the newer edgeIdeal the ideal lives in a different polynomial ring than the old (differing
	--by the generators of R) which changes dimension computations
--	n:=numgens R;
     	n:=0;
	<<"The vertices:  "<<endl;
	scan(#L,i->(
	 	<<i<<"     "<<gens ideal L#i<<endl;
	));
	--Open file for output to Graphviz
	fp:=openOut "Tgraph.dot";
	fp<<"digraph G"<<endl<<"{"<<endl;
	<<endl<<endl;
        ValidEdges:={};
	FatEdges:={};
	ColorList:=new MutableHashTable;
	ColorList#2="red";
	ColorList#3="green";
	ColorList#4="blue";
	ColorList#5="purple";
	ColorList#6="orange";
	scan(#L,i->(
	   scan(#L-i-1,j->(
		if not(isIdeal Edges#{i,i+j+1}) and #Edges#{i,i+j+1}>0 then (
 		   Eset:=Edges#{i,i+j+1};
--<<Eset<<endl;		   
   		   dimI:=dim(Eset#3);
		   --Note that Eset#3 lives in a ring with x and z variables, so the dimension are off.
		   if dimI>n then (
			if verbose then ( <<"("<<i<<", "<<i+j+1<<"): Direction "<<toString Eset_0<<endl<<toString Eset_1<<endl<<toString Eset_2<<endl<<toString Eset_3<<endl;
				<<"Dimension: "<<dimI-n<<endl<<endl;
if dimI<n then 
<<"("<<i<<", "<<i+j+1<<") : "<<toString L#i<<"    "<<toString L#(i+j+1)<<endl<<endl<<endl;
			);
			ValidEdges=append(ValidEdges,{i,i+j+1});
			fp<<i<<" -> "<<i+j+1<<"    [label = \"";
			if dimI>n+1 then fp<<dimI;
			fp<<"("<<(Eset_0)_0<<","<<(Eset_0)_1<<")\"";
			if dimI>n+1 and dimI<n+7 then fp<<" color = "<<ColorList#(dimI-n);
			if dimI>n+1 and dimI>n+6 then fp<<" color = "<<ColorList#6;
			fp<<"]; "<<endl;
			if dimI>n+1 then FatEdges=append(FatEdges,{i,i+j+1});
		   );
		);
	   ));
	));
	fp<<"}"<<endl;
	close fp;
	<<endl<<"Edge summary: "<<toString ValidEdges<<endl;
	if #FatEdges>0 then (<<"Fat Edges: "<<toString FatEdges<<endl )
	else <<"No fat edges."<<endl;
	return(ValidEdges);
);



--Work them together
-- (this should work as for the edges - we don't need to compute
-- anything, just write them down... )(plus the determinantal
-- border basis ones)


-- Then intersect the the U_mu variables.



--Check if 0 lies in the variety.

---------------------------------------------------------------
--Procedure used in arrowMaps that works out the conditions in a given degree  Input: HashTables Mlength,
--Nlength defined in arrowMaps, and Mbasi, Nbasi, which are subsets of
--their keys, a list mapList that is a partial arrowMap, and m in M of
--degree i to start with 
--Output: A list, each element of which is a list of possible degree i parts of
--the arrow-maps.
---------------------------------------------------------------

allDegiArrows=(Mlength,Nlength,Mbasi,Nbasi,mapList,m)->(
--Need to fix here to demand that the arrow map is a bijection!!!     
     alreadyChosen:=apply(mapList,p->p_1);
--<<"Showing mapList   "<<mapList<<endl;     
     possibleNs:=select(Nbasi,n->(not(n>Mlength#m) and
     	       not(Nlength#n>m) and not(member(n,alreadyChosen)) and not(n<m)));
--    <<m<<"   "<<possibleNs<<endl;
     newmapList:=apply(possibleNs,  nn->(
     	  	append(mapList,{m,nn})
      ));
      if not(m==last(Mbasi)) then (
      	   mnew:=select(Mbasi, mm->(mm<m));
     	   mnew=mnew_0;
	   newermapList:={};
      	   scan(newmapList,mList->(
	     	     newermapList=newermapList | allDegiArrows(Mlength,Nlength,Mbasi,Nbasi,mList,mnew)));
	   newmapList=newermapList;
      );
--      else newmapList ={newmapList}
--      newmapList=flatten newmapList;
-- <<"newmapList   "<<newmapList<<endl;
      return(newmapList);
);      


---------------------------------------------------------------
--List all arrow maps between a pair of monomial ideals 
--  in the sense of Hering-Maclagan
---------------------------------------------------------------

--Possible bugs:
--     Check how it handles nonstandard gradings
--   ???Would it be better to rewrite this to work just on the std monomials?

arrowMaps=method();

arrowMaps(Ideal, Ideal) := (M,N) ->(
     R:=ring M;
     if numgens R != 2 then error("This command only works for ideals in two variables");
     Ibasis := new MutableHashTable;
     Ibasis#M=apply(flatten entries basis(R/M), m->lift(m,R));
     Ibasis#N=apply(flatten entries basis(R/N), m->lift(m,R));
     v:=findEdgeDirection(M,N,Ibasis,false);
     return(arrowMaps(M,N,v,Ibasis));
);

  
arrowMaps (Ideal, Ideal, List, MutableHashTable) := (M,N,v,Ibasis) -> (
     local pset; local j; local d; 
     R:=ring M;
     --This assumes two variables.  Make assume artinian a well.
     Mbasis:=Ibasis#M;
     Nbasis:=Ibasis#N;
     a:=-v#1; b:=v#0;
     if v#0<0 then (a=-a; b=-b;);
     if (a==0 or b==0 or a*b<0) then return({});
     yy:=symbol yy;
     S:=QQ[yy_1,yy_2,Degrees=>{a,b}];
     f:=map(S,R,{yy_1,yy_2});
     M=f M;
     N=f N;
     Mbas:=new MutableHashTable;
     Nbas:=new MutableHashTable;
     Mbasis=apply(Mbasis,m->(f(m)));
     Nbasis=apply(Nbasis,m->(f(m)));
     scan(Mbasis,m->(
	       d=(degree(m))_0;
	       if Mbas#?d then Mbas#d=append(Mbas#d,m)
	       	    else Mbas#d = {m};
     ));      
     scan(Nbasis,m->(
	       d=(degree(m))_0;
	       if Nbas#?d then Nbas#d=append(Nbas#d,m)
	       	    else Nbas#d = {m};
     ));      
    --Check for same Hilbert function and replace standard monomials
    -- by elements in the ideal.
    sameHilbertFunction:=true;
    Sbas:=new MutableHashTable;
    scan(keys Mbas, d->(if not(Nbas#?d) then sameHilbertFunction=false;
	                if sameHilbertFunction then 
			     if #(Mbas#d)!=#(Nbas#d) then sameHilbertFunction=false;
			if sameHilbertFunction then (
			     Sbas#d=flatten entries basis(d,S);
			     Mbas#d = select(Sbas#d,m->(not(member(m,Mbas#d))));
     			     Nbas#d = select(Sbas#d,m->(not(member(m,Nbas#d))));
			);
   ));			
   if not(sameHilbertFunction) then return({});
   minDegree:=min(apply(flatten entries gens M, m->((degree(m))_0)));
   Akeys:=sort apply(max keys Mbas - minDegree,i->(i+minDegree));
   --Any extra degree we've added has no standard monomial in that degree
   scan(Akeys,d->(
	     if not(Mbas#?d) then (
		  Sbas#d=flatten entries basis(d,S);
		  Mbas#d=Sbas#d;
		  Nbas#d=Sbas#d;
	     );
   ));
   scan(keys Mbas, d->(
	     if #(Mbas#d)==0 then (
		  remove(Mbas,d);
		  remove(Nbas,d);
	     );
   ));	
   Akeys=sort keys Mbas;
   Alist:={new MutableHashTable};
    --Next bit may change - would we rather have an error here?
    Mlength:= new MutableHashTable;
    Nlength:= new MutableHashTable;
--Don't forget to check if keys exist before looking them up!!!???    
    --Now for each degree work out the possible pairs
    --We start at the lowest degree and work up.  Note that Akeys is sorted from smallest to largest.
    --This is a breadth first search through all arrowmaps
--<<Akeys<<endl;
    for i in Akeys do (
 	 newAlist:={};
      	 Mbasi:=reverse sort Mbas#i;
--<<"Checking the order is decreasing wrt x>y"<<endl<<Mbasi<<endl;	 
     	 Nbasi:=sort Nbas#i;
	 for A in Alist do (
	      --First work out the lower and upper bound
 --Need to add in the order if there isn't a lower arrow to do it for us!!!???	      
              for m in Mbasi do (
     	      	   mexps:=flatten exponents m;
		   mmin:=last(Nbasi);
	      	   if mexps_0>0 then (
		   	moverxdeg:=(degree(lift(m/S_0,S)))_0;
   		   	if member(moverxdeg,Akeys) then (
			     pset=select(A#moverxdeg,p->(p_0==lift(m/S_0,S)));
			     if #pset>0 then mmin=min(mmin,(pset_0)_1*S_0);
       	      	   	);
        	   );
	      	   if mexps_1>0 then (
			moverydeg:=(degree(lift(m/S_1,S)))_0;
   		   	if member(moverydeg,Akeys) then (
			     pset=select(A#moverydeg,p->(p_0==lift(m/S_1,S)));
    			     if #pset>0 then mmin=min(mmin,(pset_0)_1*S_1);
       	      	        );
       	      	   );
		   Mlength#m=mmin;
	      );
	      for n in Nbasi do (
     	      	   nexps:=flatten exponents n;
		   nmin:=last(Mbasi);
	      	   if nexps_0>0 then (
		   	noverxdeg:=(degree(lift(n/S_0,S)))_0;
   		   	if member(noverxdeg,Akeys) then (
                             pset=select(A#noverxdeg,p->(p_1==lift(n/S_0,S)));
			     if #pset>0 then nmin=max(nmin,((pset_0)_0)*S_0);
       	      	   	     );
	      		);
	      	   if nexps_1>0 then (
			noverydeg:=(degree(lift(n/S_1,S)))_0;
			if member(noverydeg,Akeys) then (
			     pset=select(A#noverydeg,p->(p_1==lift(n/S_1,S)));
			     if #pset>0 then nmin=max(nmin,(pset_0)_0*S_1);
       	      	   	     );
       	      		);
     	           Nlength#n=nmin;
	      );
     	   --Now add the arrows.
              arrowList := allDegiArrows(Mlength,Nlength,Mbasi,Nbasi,{},first(Mbasi));
              for arrows in arrowList do (
		if #arrows>0 then (
		     Anew:=new MutableHashTable;
		     j=0;
		     while Akeys#j<i do (
		     	  Anew#(Akeys#j)=A#(Akeys#j);
		     	  j=j+1;
			  );     
		     Anew#i=arrows;
		     newAlist=append(newAlist,Anew);
	       );
    	      );
     	 );     
         Alist=newAlist;
     );	  
     if #Alist==1 and #(keys (Alist#0))==0 then Alist={};
     return(Alist);
     
);     


---------------------------------------------------------------------------
-- DOCUMENTATION
---------------------------------------------------------------------------

beginDocumentation()

document { 
     Key => TEdges,
     Headline => "T-Edges in multigraded Hilbert schemes",
     "This package contains routines for the work of Milena Hering and
Diane Maclagan on T-Edges in multigraded Hilbert schemes.  The
n-dimensional torus action T on A^n extends to an action on the
multigraded Hilbert scheme.  The fixed points of this T-action are the
monomial ideals M.  The closures of the one-dimensional T-orbits are
multigraded Hilbert schemes parameterizing ideals with fixed Hilbert
function with respect to some grading Z^n/Zc.  Some commands are only
for the case n=2, where the multigraded Hilbert scheme is the Hilbert
scheme of points in the plane.",
Subnodes => {
            TO "How to use this package",
	              TO "Examples"
}
}
			      

			      
-- PARA{},
-- " In ", EM "Incidence relations among the Schubert cells of
-- equivariant punctual Hilbert schemes"," Math. Z. 242, 743-759 (2002),
-- Evain defines C(lambda) to be the set {I in H : in_{y>x}(I)=M_lambda},
-- and gives necessary conditions for overline{C(lambda)} cap C(mu) to be
-- nonempty.  In ongoing work we extend this to an analogous condition
-- for M_lambda, M_mu to live in the closure of a T-edge.",
-- PARA{},
-- "This code computes whether there is an edge joining two partitions,
-- and whether M_{mu} lives in the closure of C(lambda)."
--           }  


document {
  Key => "How to use this package",
  "The main commands are ", TO ComputeGraph, " and ", TO writeGraph,
  ".", PARA{}, "ComputeGraph will compute the T-graph of a multigraded
  Hilbert scheme.  The command writegraph calls ComputeGraph, and then
  turns the output into a GraphViz input file for easier
  visualization.  Thus if just the data is required, use ComputeGraph,
  while if you want the T-graph and the GraphViz file use
  writeGraph.", "The package also has a command, ", TO arrowMaps, " to
  compute all arrow maps between a pair of monomial ideals, and a
  command ,", TO posetLess, " to determine the relation between two
  monomial ideals in the poset order.", PARA{},

"The last part of the package consists of some data files for examples
  in the Hering-Maclagan paper.  See ", TO "Examples", " for more
  information about this."
  }


document { 
     Key => sigArrows,
     Headline => "compute the significant arrows",
     Usage => "sigArrows(I,dir)",
     Inputs => { "I"=>Ideal,
	  String=>"dir, which is  '+', '-', or 'all' giving the direction of the
significant arrows"},
     Outputs=> { MutableHashTable => " with keys
the generators of I, and entries a list of the possible heads of
significant arrows with that tail."},
     "For a monomial ideal M in a multigraded Hilbert scheme H of
points in the plane, the set of significant arrows index a basis for
the tangent space to H at M.  They were introduced by Evain in ", 
EM "Irreducible components of the equivariant punctual Hilbert schemes ",
" but we follow the naming convention and notation of Maclagan and
Smith,",  EM "Smooth and irreducible multigraded Hilbert schemes.",
PARA{},
"See Definition 4.1 of the Hering/Maclagan paper for the precise definition.",
PARA{},

"The direction + means that the f(m)<m in the order on the polynomial
ring, where the arrow points from m to f(m).  The direction - means
that f(m)>m in the order on the polynomial ring.  Note that + or - is
entered as a string.", PARA{},
"The output of an arrow map is a hash table with keys the minimal
generators of the monomial ideal I.  The entry corresponding to a
monomial generator m is a list of all monomials n for which f(m)=n is
a significant arrow.",
     EXAMPLE lines ///
     R=QQ[x,y];
     M=ideal(x^4,x^2*y,y^2);
     Sig=sigArrows(M,"+")
     peek Sig
     ///,
}       

document { 
     Key => Clambda,
     Headline => "compute the Groebner basin C(lambda)",
     Usage => "Clambda(M,+)",
     Inputs => { Ideal => "M, in a polynomial ring R in two variables
     with finitely many standard monomials",
     	  String=>"dir, which is + or -.  This determines the term
     	  order on R.  If the variables in R are R_0 and R_1, and
     	  deg(R_0)=a, deg(R_1)=b, then sign + means R_0^b<R_1^a, while
     	  - means R_0^b>R_1^b.  If R=QQ[x,y] with the standard
     	  grading, then + means y>x, while - means y<x.",},
     Outputs => { Ideal=> "J, which is the ideal of the universal
     family over the Grobner basin (all ideals with initial ideal M)"},
       "This procedure is only for two variables.", PARA{},
          "The ideal J is the ideal I_{prec}(M) of Definition 4.2 of
     the Hering-Maclagan paper, where prec is determined by + or - as
     above.  The variables z_i correspond to the positive significant
     arrows.", PARA{},
          "Warning: This does not naturally construct the reduced
          Groebner basis, but instead the generators f_i of Definition
          4.2.",
     EXAMPLE lines ///
     R=QQ[x,y]
     M=ideal(x^4,x*y,y^3)
     K=gens Clambda(M,"+")
     ///,
     "By contrast, the reduced Groebner basis is the following:",
     EXAMPLE lines ///
     S=QQ[y,x,z_1,z_2,MonomialOrder=>Lex]
     K2=substitute(K,S)
     gens gb K2
     ///,
     "Note the different polynomial with leading term y^3"
}

-- document { 
--      Key => edgeIdeal,
--      Headline => "compute the equations for the edge joining two ideals",
--      Usage => "edgeIdeal(I,J)",
--      Inputs => { "I" => Ideal, "J" => Ideal},
--      Outputs => { Set=>"of polynomials in auxiliary variables z"},
--      "The ideal I should be smaller in the poset order than J",
--      "CAVEAT!!!.  For now I'm assuming that the ideal being 'more
--      lexward' makes it smaller in the poset order", " so (x,y^3) <
--      (x^3,y)",
--      EXAMPLE lines ///
--      R=QQ[x,y]
--      I=ideal(x^2,y^3)
--      J=ideal(x^3,y^2)
--      edgeIdeal(I,J)
--      ///
-- }

document { 
          Key => {listOfMonomialIdeals, (listOfMonomialIdeals,ZZ,ZZ),
     (listOfMonomialIdeals,ZZ,ZZ,Ring), (listOfMonomialIdeals,
     Ideal), (listOfMonomialIdeals, Ideal, ZZ)},
          Headline => "compute all monomial ideals with d standard monomials",
     Usage => "listOfMonomialIdeals(d,n)",
     Inputs => {ZZ=>"d, being the number of standard monomials",
	  ZZ=>"n, being the number of variables"},
     Outputs=> { List => "of all ideals with d standard monomials in n
generators"},
     EXAMPLE lines ///
     listOfMonomialIdeals(4,2)
     ///,
     "An optional argument lets the user determine the ring in which
     the ideals live.",
     EXAMPLE lines ///
     R=QQ[x,y]
     listOfMonomialIdeals(4,2,R)
     ///,
     "An alternative usage is listOfMonomialIdeals(I), where I is an
ideal.  This will compute all monomial ideals with the same Hilbert function as 
I.  This assumes that the ideal I is Artinian.",
    EXAMPLE lines ///
    use R;
    I=ideal(x^3,y^2)
   listOfMonomialIdeals(I)
    ///,
        PARA{},"If the ideal I is not Artinian, but is homogeneous with
    respect to a positive grading, then listOfMonomialIdeals(I,d) will
    compute all monomial ideals with the same Hilbert function as I in
    degrees up to d.",
    EXAMPLE lines ///
    use R;
    I=ideal(x^2*y);
    listOfMonomialIdeals(I,4)    
    ///
  }

document {
         Key => {ComputeGraph, (ComputeGraph, ZZ, ZZ), (ComputeGraph, ZZ,
     ZZ, Ring), (ComputeGraph, Ideal, Boolean), (ComputeGraph, List,
     Boolean), (ComputeGraph, List, Boolean, Boolean)},
     Headline => "Compute the T-graph on monomial ideals",
     Usage => "ComputeGraph(d,n)",
     Inputs =>{ZZ=>"d, being the number of points",
	  ZZ=>"n, being the number of variables",
	  Ring=>"R, being the desired coefficient ring of the ideals"},
     Outputs => {List=>"of the form {List of monomial ideals, MutableHashTable}"},
     "This is the main procedure of the package.  The standard version
     takes integers d and n, and returns a list {L, EdgeLists}, where
     L is a list of all monomial ideals in n variables with d standard
     monomials, and EdgeLists is a MutableHashTable with keys pairs
     {i,j} indexing the pairs of monomial ideals in L.  The entry
     EdgeLists#{i,j} is either ideal(1), or a list with four
     entries.  The first of these is a list c giving a vector for
     which L#i and L#j share the same Hilbert function with respect to
     the Z^n/c grading.  The second and third are the ideals in
     R[yy,aa] that have initial ideals L#i and L#j, and the fourth
     ideal is the ideal in the aa variables recording the relations
     that must be satisfied for these ideals to have the correct
     Hilbert function.  Thus the fourth ideal defines the edge-scheme
     E(M,N).",
     EXAMPLE lines ///
      A=ComputeGraph(4,2) 
      A#0 
      peek A#1 
      ///, 
      PARA{},
      "The output above shows that in this graph there is an edge
      joining any pair of ideals except (0,3), and (1,4).  Compare
      Figure 5 of the Hering-Maclagan paper. ", PARA{},
      "Other formats are: ",
     UL { 
              "ComputeGraph(d,n,R), where d and n are integers, and R
      is the coefficient of the polynomial ring, so the ideals live in
      R[x_1,...,x_n].  ComputeGraph(d,n,QQ) is the same as
      ComputeGraph(d,n).",
       "ComputeGraph(I,graded), where I is an ideal in a polynomial
       ring, and graded is a boolean indicating whether the ideal is
       homogeneous with respect to the grading on the polynomial ring.
       This returns the graph on the set of ideals with either the
       same Hilbert function as the ideal I (if graded is true) or the
       same number of standard monomials (if graded is false).  Note
       that if graded is false, I must be an artinian ideal.",
              "ComputeGraph(L,graded), where L is is a list of monomial
        ideals, and graded is a boolean.  This computes the graph on
        this list of ideals, subject to the same constraint on Hilbert
        function or number of standard monomials as above.",

	"ComputeGraph(L,Slow,graded), where L is a list, and Slow and
	graded are booleans.  ComputeGraph(L,false, graded) is the
	same as ComputeGraph(L,graded).  If the boolean Slow is true,
	then the output is printed to a file EdgeList.out.  This
	allows restarting large computations that crash for lack of
	memory.",
	  },
       
           }

document {Key => findTedge,
  Headline => "find the edge joining a pair of monomial ideals",
  Usage => " findTedge(I,J,Ibasis,Graded, vect)",
  Inputs => {Ideal => "I", Ideal=>"J",
  MutableHashTable=>"Ibasis",Boolean=>"Graded", List=>"vect"},
  Outputs => {"Either ideal(1), or a list {vect, I2, J2, Mideal}" },
"The ideals I and J are monomial ideals.  The MutableHashTable Ibasis
has keys including I and J.  If I is artinian, then Ibasis#I is the
list of standard monomials of I.  If not, then Ibasis#I can be taken
to be {} (when called in the course of other code Ibasis#I will be the
list of the standard monomials of I in the same degrees as the
generators of I). The input vect is a list representing a vector c for
which I and J are both homogeneous with respect to the Z^n/c grading
with the same Hilbert function.  This can also be left as an empty
list (see example below), in which case the code will work out such a
vector.", PARA{},
"If the output has size four, the first entry, vect, is a list
     representing a vector c for which I and J  share
     the same Hilbert function with respect to the Z^n/c grading.  The
     second and third, I2 and J2,  are the ideals in R[yy,aa] that have initial
     ideals I and J, and the fourth ideal is the ideal in the aa
     variables recording the relations that must be satisfied for
     these ideals to have the correct Hilbert function.  Thus the
     fourth ideal defines the edge-scheme E(M,N).", 
     EXAMPLE lines ///
     R=QQ[x,y]
     I=ideal(x^3,y^2);
     J=ideal(x^2,y^3);
     Ibasis = new MutableHashTable;
     Ibasis#I=apply(flatten entries basis(R/I), m->lift(m,R));
     Ibasis#J=apply(flatten entries basis(R/J), m->lift(m,R));
     findTedge(I,J,Ibasis,false,{})
     /// 
     }
   
  
document { Key => findEdgeDirection,
     Headline => "find the possible edge direction between two monomial ideals",
     Usage => "findEdgeDirection(I,J)",
     Inputs =>{Ideal =>"I", Ideal=>"J",MutableHashTable=>"Ibasis",Boolean=>"Graded"},

     Outputs =>{List=>"c that is the vector for which I and J might have the
     same Hilbert function with respect to the Z^n/c grading, if that exists"},

     "The ideals I and J should both be monomial ideals.  The
     MutableHashTable Ibasis has among its keys I and J.  If
     Graded=false, then Ibasis#I is the list of standard monomials of
     I.  If Graded=true, then Ibasis#I can be taken to be {} (when called in the course of other code Ibasis#I will be  the list of the standard
     monomials of I in the same degrees as the generators of I. ", PARA{},
     "The first example is ungraded:",
     EXAMPLE lines ///
     R = QQ[x,y]
     I = ideal(x^3,y)
     J = ideal(x^2,x*y,y^2)
     Ibasis = new MutableHashTable;
     Ibasis#I = apply(flatten entries basis(R/I), m->lift(m,R));
     Ibasis#J = apply(flatten entries basis(R/J), m->lift(m,R));
     findEdgeDirection(I,J,Ibasis,false)
     ///, PARA{},
     "We now consider a graded example:",
     EXAMPLE lines ///
     R = QQ[x,y,z]
     I = ideal(x*y)
     J = ideal(x*z)
     findEdgeDirection(I,J,{},true)
     ///
     }
     

document {
     Key => {writeGraph, (writeGraph, ZZ, ZZ), (writeGraph, ZZ,ZZ, Ring, Boolean), (writeGraph, List, Boolean,Boolean) },
     Headline => "write GraphViz output for the T-graph on monomial ideals",
     Usage => "writeGraph(d,n)",
     Inputs =>{ZZ=>"d, being the number of points",
	  ZZ=>"n, being the number of variables"},
     Outputs => {"This creates a file Tgraph.dot containing GraphViz
     input illustrating the graph, and writes the edges to the screen
     if the verbose option (discussed below) is chosen"},
     EXAMPLE lines ///
     writeGraph(4,2)
     ///,
     PARA{},
     "Other formats are: ",
     UL { 
       "writeGraph(d,n,K,verbose), where d and n are integers as
     above, K is the desired coefficient ring, so the polynomial ring
     is K[x_1,...,x_n], and verbose is a boolean indicating how much
     output should be printed to the screen (true means more).  The
     case writeGraph(d,n,QQ,false) is the same as writeGraph(d,n).",
    "writeGraph(L,verbose,graded), where L is a list of monomial
    ideals, verbose is a boolean indicating whether more output should
    be printed, and graded is a boolean indicating whether the ideals
    in L are homogeneous with respect to the grading of the polynomial
    ring.",
     },
   "For comparison, here is the output of the more verbose choice:",
   EXAMPLE lines ///
   writeGraph(4,2,QQ,true)
   ///
   }


document { 
    Key => {posetLess, (posetLess, Ideal, Ideal, ZZ), (posetLess, Ideal,
     Ideal, List, MutableHashTable)},
     Headline => "determines the poset order on monomial ideals",
     Usage => "posetLess(I,J,bound)",
     Inputs=> {"I" =>Ideal, "J"=>Ideal, "bound" => ZZ, "or",
       "I" =>Ideal, "J"=>Ideal, "v"=>List, "Ibasis"=>MutableHashTable }, 
     Outputs => {Boolean},
     "This function takes two monomial ideals I and J in the same
     (positively) graded polynomial ring, and returns true if I is
     less than J in the poset order.  This is the order where I is
     less than J if there is a degree preserving bijection f from the
     monomials of J to the monomials of I with the property that f(m)
     leq m in the term order on the ring.", PARA{},
     "While the input ideals I
     and J are monomial, we assume they are given as ideals, not as
     monomialIdeals.",PARA{},
     "The first version of this function, which has syntax
     posetLess(I,J,bound), takes an additional integer input, and only
     checks up to the degree given by this bound.", PARA{},
     "The second version of this function, which has syntax
     posetLess(I,J,v,Ibasis), where v is a list representing a vector
     in Z^2, assumes that the polynomial ring has only two variables,
     that I and J are artinian, and that I and J have the same Hilbert
     function with respect to the Z^2/v grading.  Ibasis is a
     MutableHashTable whose keys include I and J.  The entries of
     Ibasis#I are the standard monomials of I, and similarly for
     Ibasis#J.",
     EXAMPLE lines ///
     R=QQ[x,y]
     I=ideal(x^4,y^2)
     J=ideal(x^2,y^4)
     posetLess(I,J,8)     
     Ibasis = new MutableHashTable;
     Ibasis#I=apply(flatten entries basis(R/I), m->lift(m,R));
     Ibasis#J=apply(flatten entries basis(R/J), m->lift(m,R));
     posetLess(I,J,{1,-1},Ibasis)
     ///
}     
     
document {
     Key => {arrowMaps, (arrowMaps, Ideal, Ideal), (arrowMaps, Ideal, Ideal, List, MutableHashTable)},
     Headline => "computes the list of arrow maps joining two monomial ideals",
     Usage => "arrowMaps(M,N)",
     Inputs =>{Ideal=>"M and ", Ideal=>"N, both of which are monomial
     ideals with the same (multigraded) Hilbert function in a
     polynomial ring in two variables."},
     Outputs => {List=>"of MutableHashTables"}, 
     "Each hash table A represents an arrow map from M to N.  It has
     as keys the degrees in which M and N differ.  The entry of A in
     degree k is a list of pairs {m,n} of monomials, indicating that the monomial
     m of degree k in M is taken to the monomial n of
     degree k in N.",
     EXAMPLE lines /// 
     R=QQ[x,y]
     M=ideal(x^3,x*y,y^2)
     N=ideal(x^2,x*y,y^3)
     A=arrowMaps(M,N);
     peek A_0
     peek A_1
     ///
}     


document { 
     Key => "Examples",
     "We illustrate here some of the computations that are discussed
     in the companion Hering/Maclagan paper.",
     Subnodes => {
       	  TO "4 points in A^2",
	  TO "2 points in P^2",
	  TO "8 points in A^4",
	  TO "Points in the plane"
     },
}

document {
     Key =>"4 points in A^2",
     "This illustrates the example of Section 5.1 of the Hering-Maclagan paper.",
     EXAMPLE lines ///
     writeGraph(4,2)
     ///,
     PARA{},
     "The output 'fat-edge' indicates that the dimension of the edge
     is greater than one",
}

document {
     Key =>"2 points in P^2",
     "This illustrates the example of Section 5.2 of the
     Hering-Maclagan paper.  Note that the ideals are not given in
     their saturated form here.",
     EXAMPLE lines ///
     R=QQ[x_0,x_1,x_2];
     I=ideal(gens image basis(2, ideal(x_0,x_1*x_2)));
     L=listOfMonomialIdeals(I,2);
     G=ComputeGraph(L,true)
     keys G_1
     ///  
}

document {
     Key =>"8 points in A^4",
     "There are 684 monomial ideals in 4 variables of colength 8.
     These are found in the file ", TT "ideals8pointsin4vars.dat", "
     available from ",PARA{},
     "www.warwick.ac.uk/staff/D.Maclagan/papers/TEdges.html", PARA{},
     "The list of 9278 edges are found in the file ", TT  "edges8pointsin4vars.dat", ".", PARA{}, "The first of these can
     be recreated by the command ", TT "listOfMonomialIdeals(8,4)", ".
     With limitless memory the second could be computed with ", TT  "ComputeGraph(8,4)", "."
}

document { Key =>"Points in the plane", 
    "The file ", TT "PointsinPlaneData.tar.gz", " available from ",
     PARA{}, "www.warwick.ac.uk/staff/D.Maclagan/papers/TEdges.html",
     PARA{}, " contains the data for the table in section 5.4 of the
     Hering/Maclagan paper.", PARA{}, "This is a set of files with
     names ", UL { "ideals10pointsIn2variables.dat",
     "arrowmap10pointsin2vars.dat",
     "arrowmapduals10pointsin2vars.dat",
     "posetorder10pointsIn2vars.dat", "edges10pointsin2vars.dat" },
     "where 10 can be replaced by any number up to 16.", PARA{}, "To
     unpack at the command-line (on unix/macos) type ", TT "gunzip
     PointsinPlaneData.tar.gz", " then ", TT "tar xvf
     PointsinPlaneData.tar", ".  This will create a directory called
     ", TT "Data", " containing these files.  These can be loaded into
     M2 using the command ", TT "load deals10pointsIn2variables.dat",
     "."
         }
 


--     R=QQ[x_1..x_4]   
--     I=ideal(x_1^8,x_2,x_3,x_4)
--     listOfMonomialIdeals(I)
     
---------------------------------------------------------------------------
-- TEST
---------------------------------------------------------------------------

TEST ///
R=QQ[x,y]
M=ideal(x^5,x*y^2,y^4);
N=ideal(x^3,y^4);
assert(#arrowMaps(M,N)==2)
///

end