-- Implementation in Macaulay 2 of algorithm 4.4.2 
-- of SST "Groebner deformations..."
--
--
-- ROUTINE: biggestMonIdeal 
-- input: an ideal in polynomial ring with n variables
-- output: the largest monomial ideal in that ideal
-- 
-- (non-trivial fact: if I is a toric initial ideal, possibly 
-- non-monomial if c non-generic, then this largest monomial ideal describes 
-- the set of non-optimal solutions to IP_{A,c})
--
--
-- first: a short preliminary procedure. Given a list of polynomials, the 
-- keepMonsOnly procedure keeps the monomials from that list only.
--

keepMonsOnly = (listOfPolys) -> (
  monsOnly := {};
  m = #listOfPolys;
  scan (0..(m-1), j ->(
	     polynom = listOfPolys_j; 
	     if # exponents polynom == 1 then monsOnly = append (monsOnly, polynom);
	     ));	
     monsOnly
     );


biggestMonIdeal = (I)  -> (
	S = ring I; 
	n := # gens S; 
	multGensList := {};
	bigS = QQ[(vars(53..52+n), gens S), MonomialOrder => Lex];
	genbigS = gens bigS;
	irrM = toList(apply(n, q -> genbigS_{n-1+q}));
	preH = substitute(I, bigS);
	genList = flatten entries gens preH;
	t := # genList;	
	scan(t, i -> (
		newPol := genList_i;
		apply(n, j -> (
		half = toList (apply(n, k -> (if j==k then 1 else 0))); 
		full = join(half, half);
		newPol = homogenize(newPol, genbigS_j, full)
		));
	multGensList = append(multGensList, newPol)
	));
	newI = ideal multGensList; 
	scan(n, l -> (newI = saturate(newI, genbigS_l );));
	bigGBList = flatten entries gens gb newI;
	listOfMonomials = keepMonsOnly bigGBList;
	monomialIdeal substitute(monomialIdeal(listOfMonomials), S)
	);