<< ComputationalGeometry`
Clear[CAPlot,AmPlot, LopsidedCAPlot, GKZVol, ZonotopeVol, MellinBarnesQ, VerticesQ, PrimitiveEdgesQ, PlotShell];

CAPlot[func_]:= Module[{rmax,rstep,tstep,lattice,solutions,internSolver,r,t,z,w,ExponentList,infunc,minExp},

 ExponentList = Map[Function[{x},{Exponent[Expand[func[z,w]][[x]],z],Exponent[Expand[func[z,w]][[x]],w]}],Table[n,{n,1,Length[Expand[func[z,w]]]}]];
 minExp = Min/@ (ExponentList\[Transpose]);


 infunc[z_,w_]:= Expand[func[z,w]/(z^minExp[[1]] w^minExp[[2]])];

 rmax = 4;
 rstep = .2;
 tstep = Pi/60;
 
 lattice = Flatten[Table[Exp[r + I t],{r,-rmax,rmax,rstep},{t,0,2Pi,tstep}]];
 
 internSolver[pair_]:= Module[{sols},
   sols=z/.NSolve[pair[[1]]==0,z];
   Return[{sols,Table[pair[[2]],{n,1,Length[sols]}]}\[Transpose]];
 ];

 solutions = Select[Join[Flatten[internSolver/@({infunc[z,lattice],lattice}\[Transpose]),1],
 Reverse/@Flatten[internSolver/@({infunc[lattice,z],lattice}\[Transpose]),1]],#[[1]]≠0&&#[[2]]≠0&];

 Return[ListPlot[Function[{x},Arg/@ x]/@solutions,AspectRatio->1,Axes->False]];

];

AmPlot[func_]:= Module[{rmax,rstep,tstep,lattice,solutions,internSolver,r,t,z,w,ExponentList,infunc,minExp, LogAbs},

ExponentList=Map[Function[{x},{Exponent[Expand[func[z,w]][[x]],z],Exponent[Expand[func[z,w]][[x]],w]}],Table[n,{n,1,Length[Expand[func[z,w]]]}]];
minExp = Min/@ (ExponentList\[Transpose]);

infunc[z_,w_]:=Expand[func[z,w]/(z^minExp[[1]] w^minExp[[2]])];

rmax = 4;
rstep=.2;
tstep=Pi/60;

lattice = Flatten[Table[Exp[r + I t],{r,-rmax,rmax,rstep},{t,0,2Pi,tstep}]];

internSolver[pair_]:=Module[{sols},
sols=z/.NSolve[pair[[1]]==0,z];
Return[{sols,Table[pair[[2]],{n,1,Length[sols]}]}\[Transpose]];];

solutions = Select[Join[Flatten[internSolver/@({infunc[z,lattice],lattice}\[Transpose]),1],
Reverse/@Flatten[internSolver/@({infunc[lattice,z],lattice}\[Transpose]),1]],#[[1]]≠0&&#[[2]]≠0&];

LogAbs[x_]:=Log[Abs[x]];

Return[ListPlot[Function[{x},LogAbs/@ x]/@solutions,AspectRatio->1,Axes->False]];

];

LopsidedCAPlot[g_]:=Module[{GetPart,z,w,myRules,allPieces},

(* Input: Trinom as coeffRules *)
(* Returns Null if there is no volume by matrix *)
GetPart[crules_]:=Module[{gobj,az,aw,theta1,theta2,exps,geomTrans,cArgs,rep,matrix,constant,triangle1,triangle2,triangles},

exps = First/@crules;
cArgs = Arg[Last/@crules];

theta1=cArgs[[2]]-cArgs[[1]]+exps[[2]].{az,aw}-exps[[1]].{az,aw};
theta2=cArgs[[3]]-cArgs[[1]]+exps[[3]].{az,aw}-exps[[1]].{az,aw};

matrix={{Coefficient[theta1,az],Coefficient[theta1,aw]},
{Coefficient[theta2,az],Coefficient[theta2,aw]}};
constant =Mod[({theta1,theta2}/.{az->0,aw->0})+100Pi,2Pi];

If[Det@matrix==0,Return[Null]];

geomTrans[gobj_]:=GeometricTransformation[gobj,InverseFunction[AffineTransform[{matrix,constant}]]];

rep=Floor[2+Max[Abs[Eigenvalues[matrix]]]];

triangles=Polygon[Flatten[Table[{
{{x-Pi,y},{x,y+Pi},{x-Pi,y+Pi}},
{{x,y-Pi},{x+Pi,y},{x+Pi,y-Pi}}},{x,-2Pi rep,2Pi rep,2Pi},
{y,-2Pi rep,2Pi rep,2Pi}],2]];

Graphics[{Purple,Opacity[0.1],geomTrans@triangles},PlotRange->Pi]//Return;

];

(* Create trinoms *)
myRules=CoefficientRules[g[z,w],{z,w}];
Return[Show[DeleteCases[GetPart/@Subsets[myRules,{3}],Null]]];

];

(*
GKZVol[A_]:=Module[{B, DetList, VectArg, IList},

 VectArg[b_]:= Arg[b.{1,I}];
 IList = Subsets[Range[5],{2}];

 If[Length[A] ≠ 3 || Length[A\[Transpose]] ≠ 5, Print["GKZVol: Error, A of wrong dimensions."]; Return[{}];,

 B = Sort[NullSpace[A]\[Transpose], VectArg[#1] ≥ VectArg[#2]&];
 DetList = Det[B[[#]]]&/@IList;
 Return[Total[Select[DetList, # > 0&]]/(GCD@@DetList)]
]];
*)

GKZVol[A_]:=Module[{AA, NN},
If[Length[A] ≠ 2, Print["GKZVol: Error, A should be a matrix with two rows."]; Return[{}];,
 AA = A\[Transpose][[ConvexHull[A\[Transpose]]]];
 NN = Length[AA];
 Return[Total[Table[Det[(AA[[{Mod[i, NN]+1,Mod[i+1, NN]+1}]])], {i, 1, NN}]]];
]];



ZonotopeVol[A_]:= Module[{B, DetList, VectArg, IList},

 VectArg[b_]:= Arg[b.{1,I}];
 IList = Subsets[Range[5],{2}];

 If[Length[A] ≠ 3 || Length[A\[Transpose]] ≠ 5, Print["ZonotopeVol: Error, A of wrong dimensions."]; Return[{}];,

 B = Sort[NullSpace[A]\[Transpose], VectArg[#1] ≥ VectArg[#2]&];
 DetList = Det[B[[#]]]&/@IList;
 DetList = DetList/(GCD@@DetList);
 Total[Abs[DetList]]
]];



MellinBarnesQ[A_]:= Module[{B, DetList, VectArg, IList},

 VectArg[b_]:= Arg[b.{1,I}];
 IList = Subsets[Range[5],{2}];

 If[Length[A] ≠ 3 || Length[A\[Transpose]] ≠ 5, Print["MellinBarnesQ: Error, A of wrong dimensions."]; Return[{}];,

 B = Sort[NullSpace[A]\[Transpose], VectArg[#1] ≥ VectArg[#2]&];
 DetList = Det[B[[#]]]&/@IList;
 DetList = DetList/(GCD@@DetList);
 If[Ceiling[Total[Abs[DetList]]/4] ≥ Total[Select[DetList, # > 0&]], Return[True], Return[False]]
]];



VerticesQ[A_] := Module[{dList, AA},

 If[Length[A] ≠ 2, Print["VerticesQ: Error, Input should be a matrix with two rows."]; Return[{}];,
   AA = A\[Transpose][[ConvexHull[A\[Transpose]]]];
   If[Length[AA]≠ Length[A\[Transpose]], Return[False], 
    dList = Differences[Append[AA, First[AA]]];
    If[Min[Abs[Det/@({dList,Drop[Append[dList, First[dList]],1]}\[Transpose])]]==0,Return[False], Return[True]];

]]];



(* Assuming A sorted *)
PrimitiveEdgesQ[A_]:= Module[{AA},

 If[Length[A]≠ 2, Print["IrreducibleEdgesQ: Error, Input should be a matrix with two rows."]; Return[{}];,
 AA = A\[Transpose][[ConvexHull[A\[Transpose]]]];
 If[Not[VerticesQ[AA\[Transpose]]], Return[False], 
   If[Max[GCD@@#&/@ Differences[Append[AA, First[AA]]]]==1, Return[True], Return[False]]

]]];

PlotShell[func_]:= Module[{Ac,A, M,a,b, PreEQ, t1,t2,a1,a2, EQ,z,w, LFormTranslations, FixedEQ},

Ac = ({First[#], Arg[Last[#]]}& /@ CoefficientRules[func[z,w], {z,w}]);
A = First/@Ac;
Ac = Ac[[ConvexHull[A]]];
A = A[[ConvexHull[A]]];

If[Not[VerticesQ[(A\[Transpose])]], Print["PlotShell: Error, polynomial is not sparse."]; Return[{}];,

M = Length[Ac];
PreEQ=Evaluate[ (Subtract@@#&/@Table[Ac[[{Mod[i,M]+1, Mod[i+1,M]+1}]], {i,1,M}])];

EQ[t1_,t2_]:= Evaluate[First[#].{t1,t2} + Last[#] + Pi&/@ PreEQ];

LFormTranslations[expr_]:=Module[{mi,ma, NN, inexpr},
NN=LCM[Denominator[D[expr, t1]], Denominator[D[expr, t2]]];
inexpr =Evaluate[Expand[ NN expr]];
{mi,ma}={Floor[Min[#]], Ceiling[Max[#]]}&@@{{inexpr/.{t1->-Pi, t2->-Pi},inexpr/.{t1->-Pi, t2->Pi}, inexpr/.{t1->Pi, t2->-Pi}, inexpr/.{t1->Pi, t2->-Pi}}/(2Pi)};
Return[Table[inexpr- k 2Pi, {k,mi-1,ma+2}]];
];

FixedEQ = Flatten[LFormTranslations/@ EQ[t1,t2]];

Return[Show[ContourPlot[#==0, {t1,-3.15,3.15},{t2,-3.15,3.15}, ContourStyle->Purple]&/@FixedEQ]];

]];

JFCoAmoebaHelp := Print["JFCoAmoebas -- included functions:
  CAPlot : Polynomial -> Approximate coamoeba picture.
  AMPlot : Polynomial -> Approximate amoeba picture.
  LopsidedCAPlot : Polynomial -> Lopsided coamoeba plot.
  PlotShell : Polynomial -> Plot of the shell of the coamoeba.
  GKZVol : GKZ matrix in dimension and codimension two -> Normalized volume of polytope.
  ZonotopeVol : GKZ matrix in dimension and codimension two -> Volume of dual zonotope.
  MellinBarnesQ : GKZ matrix in dimension and codimension two -> True if Mellin Bernes basis known to exist.
  VerticesQ : Matrix of dimension two support -> True if all points are vertices of Newton polygon.
  PrimitiveEdgesQ : Matrix of dimension two support -> True if all edges of Newton polygon are primitive."]

DumpSave["JFCoAmoebas.mx",{CAPlot,AmPlot,LopsidedCAPlot, GKZVol, ZonotopeVol, MellinBarnesQ, VerticesQ, PrimitiveEdgesQ, PlotShell, JFCoAmoebaHelp}];

(* OLD FUNCTIONS *)

(*

CoAmoebaCriticalPointsPlot[func_]:=Module[{points,z,w,RemoveMonomialFactor,f1,f2},

RemoveMonomialFactor[ifunc_]:=Module[{exponentList,mz,mw,zz,ww,g},
exponentList=Map[Function[{x},{Exponent[Expand[ifunc[zz,ww]][[x]],zz],Exponent[Expand[ifunc[zz,ww]][[x]],ww]}],Table[n,{n,1,Length[Expand[ifunc[zz,ww]]]}]];
{mz,mw} = Min/@(exponentList\[Transpose]);
g[zz_,ww_]:=Simplify[ifunc[zz,ww]/(zz^mz ww^mw)];
Return[g]];

f1[z_,w_]:=D[Expand[func[z,w]],z];
f2[z_,w_]:=D[Expand[func[z,w]],w];

points = Select[{z,w}/.NSolve[RemoveMonomialFactor[f1][z,w]==0 &&RemoveMonomialFactor[f2][z,w]==0,{z,w}],#[[1]]≠0&&#[[2]]≠0&];

Return[ListPlot[Map[Arg,points,{2}],AspectRatio->1,PlotStyle->Red]];
];

CoAmoebaCombPlot[func_]:=Show[CoAmoebaPlot[func],CoAmoebaCriticalPointsPlot[func]];

*)