(* ::Package:: *)

NewtonPolytope[poly_,vars_]:=Exponent[#,vars]&/@MonomialList[poly]
MinkowskiSum[x__,y_]:=DeleteDuplicates[Which[Length[{x}]>1,MinkowskiSum[MinkowskiSum[x],y],Length[{x}]==1,Flatten[Outer[Plus,x,y,1],1]]]
(*
prereducepoly[x_]:=If[Length[x]<=1,x,With[{min=First[MinimalBy[x,Identity]]},With[{hh=DeleteCases[x,min]},Append[prereducepoly[Select[hh,AnyTrue[#-min,#<0&]&]],min]]]]
reducepoly[poly_,vars_]:=With[{hh=prereducepoly[NewtonPolytope[poly,vars]]},Total[(Times@@(vars^#))&/@hh]]
explog[exp_]:=exp/.log[x_]:>Total[log[#[[1]]]#[[2]]&/@FactorList[x]]/. log[_?NumberQ]->0
*)
isfeasible[mat_,vec_]:=Quiet[Check[LinearProgramming[ConstantArray[0,Length[mat]],Transpose[mat],{#,0}&/@vec],False,LinearProgramming::lpsnf],{LinearProgramming::lpsnf}]

tranmat[exp_,vars_]:=With[{hh=Intersection[Variables[exp],vars]},Thread[hh->Transpose@NewtonPolytope[exp,hh]]]
newchoice[exp_,vars_]:=With[{hh=tranmat[exp,vars]},#->(Count[#,Min[#]]&[Total[#/.hh]])&/@MinimalBy[MaximalBy[Select[Subsets[First/@hh,{2,Max[2,Length[hh]]}],(exp/.Thread[#->0])===0&],(Count[#,Min[#]]&[Total[#/.hh]])&],Length]]
Options[l2choice]={"L2"->True};
l2choice[exp_,vars_,OptionsPattern[]]:=If[OptionValue["L2"],With[{hh=Select[Complement[Expand/@Cases[Collect[1+exp,_log],log[xx_]:>xx,Infinity],vars],Head[#]===Plus&&Length[#]===2&]},If[hh==={},{},First[First@Commonest[MinimalBy[Flatten[newchoice[#,vars]&/@hh,1],Last]]]]],{}]

(*prereducepoly[x_]:=If[Length[x]<=1,x,With[{min=First[MinimalBy[x,Identity]]},With[{hh=DeleteCases[x,min]},Append[prereducepoly[Select[hh,AnyTrue[#-min,#<0&]&]],min]]]]*)
prereducepolymat[x_]:=If[Length[x]<=1,x,Prepend[prereducepoly[Select[Rest[x],AnyTrue[#-First[x],#<0&]&]],First[x]]]
prereducepoly[x_]:=prereducepolymat[SortBy[x,Identity]]
reducepoly[poly_,vars_]:=Total[(Times@@(vars^#))&/@prereducepoly[NewtonPolytope[poly,vars]]]
explog[exp_]:=exp/.log[x_]:>Total[log[#[[1]]]#[[2]]&/@FactorList[x]]/. log[_?NumberQ]->0
droplog[exp_,vars_]:=Collect[Collect[explog[exp],_log]/.log[poly_]:>log[reducepoly[poly,vars]]/.log[1]->0,_log]//.log[x_^b_]:>b log[x]

findhp[mat_]:=Module[{$var$},With[{vars=Array[$var$,Length[mat[[1]]]]},vars/.First[FindInstance[Join[#>0&/@(mat . vars),{Total[vars]==0,Total[#^2&/@vars]!=0}],vars,Integers]]]]
unitvec[i_,len_]:=unitvec[i,len]=Normal@SparseArray[{{i}->1},{len}]
isout[nvec_,vec_]:=If[DuplicateFreeQ[vec],(isfeasible[Table[unitvec[#[[i+1]],Length@vec]-unitvec[#[[i]],Length@vec],{i,Length@Ordering[vec]-1}]&[Ordering[vec]],nvec])===False,False]
outvecs[mat_]:=SortBy[GatherBy[With[{hh=findhp[mat]},Select[mat,isout[hh,#]&]],Ordering],Length]
findlc[vec1_,vec2_,nvec_]:=Module[{n,k},n=0;k=0;While[k==0&&n<30,n++;Do[If[!isout[nvec,i vec1+(n+1-i) vec2],k=i;Break[];],{i,n}];];{k,n+1-k}];
getblowuplc[mat_]:=With[{pos=Position[mat,#][[1,1]]&/@(First/@outvecs[mat])},With[{hh=findlc[mat[[pos[[1]]]],mat[[pos[[2]]]],findhp[mat]]},{pos[[1;;2]],hh,Total[hh mat[[pos[[1;;2]]]]]}]];
varstep1[exp_done]:=exp
varstep2[exp_done]:=exp
varstep3[exp_done]:=exp
varstep1[vmat_]:=With[{hh=Transpose@DeleteCases[Transpose[{vmat[[1]],(#-Min[#]&/@(Transpose@prereducepoly@Transpose@vmat[[2]]))}],{_,{0..}}]},If[hh==={},done[vmat[[3]]],Append[hh,vmat[[3]]]]]
varstep2[vmat_]:=With[{mat=vmat[[2]]},With[{hh=isfeasible[mat,ConstantArray[1,Length[First@mat]]]},If[hh===False,{vmat},With[{jj=Complement[Range[Length[mat]],Flatten@Position[hh,0]],kk=hh/GCD@@hh},{Delete[vmat[[1]],#],Delete[mat,#],Append[vmat[[3]],blowup[{vmat[[1,jj]],vmat[[1,#]]},kk[[jj]]]]}&/@jj]]]]
varstep3[vmat_]:=With[{mat=vmat[[2]]},{{vmat[[1]],ReplacePart[mat,#[[1,1]]->#[[3]]],Append[vmat[[3]],blowup[{vmat[[1]][[#[[1]]]],vmat[[1,#[[1,1]]]]},#[[2]]]]},{vmat[[1]],ReplacePart[mat,#[[1,2]]->#[[3]]],Append[vmat[[3]],blowup[{vmat[[1]][[#[[1]]]],vmat[[1,#[[1,2]]]]},#[[2]]]]}}&@getblowuplc[mat]]
prevargo[vmats_]:=FixedPoint[varstep1/@Flatten[varstep2/@#,1]&,varstep1/@vmats]
vargo[vmat_]:=FixedPoint[prevargo[Flatten[varstep3/@#,1]]&,prevargo[{vmat}]]
Blowup[{vars_,var_},degs_]:={DeleteCases[#[[1]]->pow[#[[1]],#[[2]]]&/@Thread[{vars,degs}],_->pow[_,1]],#-># var&/@Complement[vars,{var}],DeleteCases[#[[1]]->pow[#[[1]],#[[2]]]&/@Thread[{vars,1/degs}],(_->pow[_,1])|(var->_)]}
transeq[exp_done]:=Flatten[List@@exp/.blowup->Blowup]

(*applyrep[logkn_,list_]:=applyrep[logkn,list]=If[Length[list]<1,logkn,applyrep[logkn,Most@list]/.(Last[list]/.pow[a_,b_]:>a^b//.log[x_^b_]:>b log[x])]*)
applyrep[logkn_,list_]:=If[Length[list]<1,logkn,applyrep[logkn,Most[list]]/. (Last[list]/. pow[a_,b_]:>a^b)//.log[x_^b_]:>b log[x]]
isdone[logkn_,vars_]:=(logkn/.Dispatch[log[#]->0&/@vars])===0;
(*pregenblowups[{logkn_,seq_,done_},vars_]:=With[{jj=Union[Cases[1+logkn,log[x_]:>x,Infinity]]},If[!done,With[{nn=droplog[applyrep[logkn,#],vars]},{nn,Join[seq,#],isdone[nn,vars]}]&/@(transeq/@vargo[{vars,Transpose@NewtonPolytope[First@MinimalBy[Complement[jj,vars],Times@@(Length/@Variables/@List@@#)&],vars],{}}]),{{logkn,seq,done}}]]*)
rescalelogkn[{logkn_,seq_,done_},vars_]:=With[{hh={#,GCD@@Exponent[Complement[Cases[1+logkn,log[xx_]:>xx,Infinity],vars],#]}&/@vars},With[{jj=DeleteCases[#[[1]]->pow[#[[1]],If[#[[2]]===0,1,1/#[[2]]]]&/@hh,(_->pow[_,1])]},With[{kk=applyrep[logkn,jj]},{kk,Join[seq,jj],isdone[kk,vars]}]]]
rescalelogkn[{logkn_, seq_, done_}, vars_] := {logkn, seq, done}
nl2genblowups[{logkn_,seq_,done_},vars_]:=If[!done,With[{jj=Union[Cases[1+logkn,log[x_]:>x,Infinity]]},With[{nn=droplog[applyrep[logkn,#],vars]},{nn,Join[seq,#],isdone[nn,vars]}]&/@(transeq/@vargo[{vars,Transpose@NewtonPolytope[First@MinimalBy[Complement[jj,vars],Times@@(Length/@Variables/@List@@#)&],vars],{}}])],{{logkn,seq,done}}]
pregenblowups[{ologkn_,oseq_,odone_},vars_]:=With[{pp=rescalelogkn[{ologkn//.log[Power[x_,b_]]:>b log[x],oseq,odone},vars]},With[{logkn=pp[[1]],seq=pp[[2]],done=pp[[3]]},If[!done,With[{hh=l2choice[logkn,vars]},If[hh==={},nl2genblowups[{logkn,seq,done},vars],With[{nn=droplog[applyrep[logkn,#],vars]},{nn,Join[seq,#],isdone[nn,vars]}]&/@{{First@hh->Times@@hh},{Last@hh->Times@@hh}}]],{{logkn,seq,done}}]]]
ppgenblowups[dones_,undones_,vars_]:=If[Length[undones]===0,dones,With[{hh=Flatten[pregenblowups[#,vars]&/@undones,1]},ppgenblowups[Join[dones,Select[hh,Last]],Select[hh,!Last[#]&],vars]]]
genblowups[logkn_,vars_]:=Most/@With[{hh=droplog[logkn,vars]},ppgenblowups[{},{{hh,{},isdone[hh,vars]}},vars]]

getleading[logkn_, 
  vars_] := {Total[First /@ #], 
    Flatten[Last /@ #] // 
     Union} &@({#[[2]]/Times @@ #[[1]], #[[1]]} & /@ 
    With[{hh = genblowups[droplog[logkn, vars], vars]},
     {Normal@Last@CoefficientArrays[#[[1]], log /@ vars], 
        Times @@ Cases[#, pow[x_, y_] :> y, Infinity]} & /@ hh])
allleading[logkn_, 
  vars_] := {Total[First /@ #], 
    Flatten[Last /@ #] // 
     Union} &@({#[[2]]/Times @@ #[[1]], #[[1]]} & /@ 
    With[{hh = 
       Flatten[genblowups[droplog[logkn /. #, vars], vars] & /@ 
         Subsets[# -> 1/# & /@ vars], 1]},
     {Normal@Last@CoefficientArrays[#[[1]], log /@ vars], 
        Times @@ Cases[#, pow[x_, y_] :> y, Infinity]} & /@ hh])
