comphp = (m,J) -> (R:=ring(J);
     N:= numgens ring(J);    
     s:=apply(N,i->0);
     L:= flatten entries basis(m,R/J);
     for i from 0 to #L-1 do s = s + sum(exponents(L_i));
     return s
);
subhp = (m,I) -> (
     N:= numgens ring(I);
     totdeg:= ((m+N-1)!*m / (m!*N!));
     st:=comphp(m,I);
     apply(N, k -> totdeg - st_k));



--functions more or less copied from StatePolytope.m2
--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
--     )

printHilbPt = (L,I) -> ( 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,I) -> (
openOut "temporarypolymakefile.txt";
polymakeinput :=  concatenate("POINTS",newline);
for i from 0 to (#statePolytopePointsList - 1) do polymakeinput = concatenate(polymakeinput,concatenate(printHilbPt(statePolytopePointsList_i,I),newline));
"temporarypolymakefile.txt" << polymakeinput << closeOut
);

polymakeToM2 = (st) -> (
p := concatenate("VERTICES",concatenate(newline,"1 "));
st = replace(p,"{{",st);
p = concatenate(newline,"1 ");
st = replace(p,"},{",st);
st = replace(" ",", ",st);
st = concatenate(st,"}}");
return value st
);
----------------------------------------------------------------

--new functions

randInId = (I) -> (
R:= ring I;
K:=coefficientRing R;
n:=numgens R;
--calculate initial ideal wrt a random weight vector
w:=apply(n,i-> random(-50,50));
Sw := K[gens R, Weights => w, MonomialOrder => GLex, Global=>false];
W := map(Sw,R,vars Sw);
inwI:= ideal leadTerm W(I);
gensinwI := gens(inwI);
V:= map(R,Sw, vars R);
L:= apply(numColumns(gensinwI), i-> V(gensinwI_(0,i)) );
output := {{L,w}};
--also compute the initial ideal wrt the opposite weight vector, for antithetical variance reduction
w=-w;
Sw = K[gens R, Weights => w, MonomialOrder => GLex, Global=>false];
W = map(Sw,R,vars Sw);
inwI= ideal leadTerm W(I);
gensinwI = gens(inwI);
V= map(R,Sw, vars R);
L= apply(numColumns(gensinwI), i-> V(gensinwI_(0,i)));
return append(output,{L,w})
);


mcStable = (m,I, maxrounds) -> (
     n := numgens ring(I);
   --  monomialsList := flatten entries(basis(m,ring(I)));
     hps:={};
     wts:={};
     r:={};
     polymakesession:="";
     st:="";
     str:="";
     barycenter:={};
     barycenterstring:="";
     augmentedstatepolytope:="";
     rounds:=0;
     while rounds <= maxrounds do (
	  for t from 0 to 2*n-1 do (
     r= randInId(I);
   --  hps = append(hps,hilbertPt(n,m,(r_0)_0,monomialsList));
   --  hps = append(hps,hilbertPt(n,m,(r_1)_0,monomialsList));
     hps = append(hps, subhp(m,ideal((r_0)_0)));
     hps = append(hps, subhp(m,ideal((r_1)_0)));
     wts = append(append(wts,(r_0)_1),(r_1)_1) );
createPolymakeInputFile(hps,I);
polymakesession = "!polymake temporarypolymakefile.txt VERTICES";
polymakesession << closeOut;
st = get polymakesession;
    str= st;
    if hps == {} then error "hps is empty!";
     if rounds == 0 then (barycenter = toString (  (sum (hps_0)) / (n ));
     if rounds == 0 then barycenterstring= concatenate("POINTS",concatenate(newline,"1 "));
     for i from 0 to (numgens(ring(I))-1) do barycenterstring = concatenate(barycenterstring, concatenate(barycenter," "));
     print "barycenter is ";
     print barycenter );
     str = replace("VERTICES",barycenterstring,str);
     openOut "augmentedtemporarypolymakefile.txt";
     "augmentedtemporarypolymakefile.txt" << str << closeOut;
     polymakesession2 = "!polymake augmentedtemporarypolymakefile.txt VERTICES";
     polymakesession2 << closeOut;
     augmentedstatepolytope = get polymakesession2;
     if set polymakeToM2(st) === set polymakeToM2(augmentedstatepolytope)   then return (true, hps,wts) else rounds = rounds +1
     );
     if rounds > maxrounds then return (false, hps, wts)
     )