BeginPackage["PolygonalSG`", {"ComputationalGeometry`"}]
(* Authors: 
I. Garc\'{i}a-Garc\'{i}a <ignacio.garcia@uca.es>, 
A. S\'{a}nchez-R.-Navarro <alfredo.sanchez@uca.es>, 
A. Vigneron-Tenorio <alberto.vigneron@uca.es>
*)

PolygonalSG::usage = 
  "PolyginalSG[p], given p, a polygon rational 2D vertex set, if \
exists, it computes the minimal system of generators of the affine \
semigroup associated to the convex hull of the polygon. \
An external program is needed. It must be specified assigning to the variable Solucionator any value similar to this example: \n \
Windows operating system: c:\\path\\dickson.exe";
BelongToSG::usage = 
  "BelongToSG[p, v], given an integer pair p, and the set of vertex \
of a 2D rational v, it returns True if p belongs to the affine \
semigroup associated to the polygon, and False otherwise.";
Solucionator = "Dropbox/DPSolve/build/DPSolve";
Begin["Private`"]
PolygonalSG[p_] := 
  Module[{i, pr = p, clws, lenclws, Ac, aux, pp, u, eq, ti, ToS, 
    fname, s, equation, xmax, ymax, b},
   If[! CheckSolucionator[], 
    Print[StringJoin["External program ", Solucionator, " needed not found."]]; 
    Return[];];
   For[i = 1, i <= Length[p], i++,
    pr[[i, 1]] = Rational2Frac[p[[i, 1]]];
    pr[[i, 2]] = Rational2Frac[p[[i, 2]]];
    ];
   Needs["ComputationalGeometry`"];
   clws = 
    ComputationalGeometry`ConvexHull[p, 
     ComputationalGeometry`AllPoints -> False];
   lenclws = Length[clws];
   If[Length[clws]<= 2, 
    Print["The number of points does not correspond to a rational polygon."];
    Return[];];
   Ac = {};
   aux = 0*Range[1, lenclws];
   aux[[1]] = 1;
   pp = pr[[#]] & /@ clws;
   pp = Join[#, {1}] & /@ pp;
   For[i = 1, i <=  Length[clws], i++,
    u = pp[[1]]\[Cross] pp[[2]];
    u = u*(LCM @@ Denominator[u]);
    eq = Join[u, -Sign[u.pp[[3]]]*aux];
    aux = RotateRight[aux];
    pp = RotateRight[pp];
    AppendTo[Ac, eq]; 
    ];
   ti = ConstantArray[0, {lenclws, 1}];(*0*Range[1,lenclws];*)
   
   fname = FileNameJoin[{$TemporaryDirectory, "testfilepolsg"}];
   s = OpenWrite[fname];
   Do[
    equation = "";
    Do[
     equation = equation <> " " <> ToString[Ac[[k, j]]],
     {j, 1, Length[Ac[[k]]]}
     ];
    equation = 
     equation <> " = " <> ToString[ti[[k, 1]]] <> 
      If[k == Length[Ac],(*then*)".\n",(*else*)",\n"];
    WriteString[s, equation];,
    {k, 1, Length[Ac]}
    ];
   Close[fname];
   ToS = Import[StringJoin["! ", Solucionator, " < "] <> fname, 
     "Table"];
   ToS = Union[ToS[[All, 1 ;; 2]]];
   pp=pp[[All, 1 ;; 2]];
   ToS = SistemaMGPol[pp, ToS];
   xmax = Max[ToS[[All, 1]]]; ymax = Max[ToS[[All, 2]]];
   b = Max[IntegerPart[xmax/Max[pp[[All, 1]]]], 
     IntegerPart[ymax/Max[pp[[All, 2]]]]];
   Print[Show[
      Graphics[{Red, Thin, 
        Table[Line[{{0, i}, {xmax+1, i}}], {i, 1, ymax+1}], Thin, 
        Table[Line[{{i, 0}, {i, ymax+1}}], {i, 1, xmax+1}], Blue, 
        Opacity[0.1], 
        Table[Polygon[i*pp[[All, 1 ;; 2]]], {i, 1, b+1}]}], 
      ListPlot[SistemaMGPol[pp[[All, 1 ;; 2]], ToS], 
       PlotStyle -> Directive[Blue, Thick]], Axes -> True, 
      AxesOrigin -> {0, 0}, PlotRange -> {{0, xmax+1}, {0, ymax+1}}, 
      PlotRangeClipping -> True, PlotRangeClipping -> True, 
      AspectRatio -> 3]];
   Return[ToS];
   ];
CheckSolucionator[] := Module[{fname, s, sols},
   fname = FileNameJoin[{$TemporaryDirectory, "testsator"}];
   s = OpenWrite[fname];
   WriteString[s, "-2 1 0 = 0,\n0 3 -2 =0."]; 
   Close[fname];
   sols = 
    Import[StringJoin["! ", Solucionator, " < "] <> fname, "Table"];
   Return[sols == {{1, 2, 3}}];
   ];
Rational2Frac[x_] := Module[{res},
   If[Element[x, Rationals],
     Return[x];,
     res = Rationalize[N[x], 0];
     If[x - res == 0, Return[res];,
      Print["Input is not a rational number."];
      Abort[];
      ];
     ];
   ];
SistemaMGPol[conjvertices_, sistGeneradores_] := 
  Module[{menorIgualProdCartesiano, esElemMinimalPoligono, sistG, 
    ecuaciones, dist},
   menorIgualProdCartesiano[a_, b_] := 
    And @@ Table[a[[i]] <= b[[i]], {i, Length[a]}];
   (*******)
   
   esElemMinimalPoligono[ecuaciones_, dist_, sistG_, elemento_] := 
    Module[{lmenores, i, longitud, seguir = True, cad}, 
     lmenores = 
      Select[sistG, (menorIgualProdCartesiano[#, 
           elemento] \[And] (# != elemento)) &];
     longitud = Length[lmenores];
     i = 1;
     While[i <= longitud \[And] seguir,
      If[perteneceASemigrupo[elemento - lmenores[[i]], ecuaciones, 
        dist], seguir = False];
      i++;
      ];
     (*Print[elemento," >= ",lmenores," ",seguir];*)
     
     Return[seguir];
     ];
   (*******)
   sistG = sistGeneradores;
   ecuaciones = ecuacionesPoligonal[conjvertices];
   dist = distAlOrigen[conjvertices];
   Select[sistG, (esElemMinimalPoligono[ecuaciones, dist, sistG, #]) &]
   ];
BelongToSG[p_, vertices_] := Module[{vertpol, pos = ConvexHull[vertices]},
  vertpol=Table[vertices[[pos[[i]]]], {i, 1, Length[pos]}];
  perteneceASemigrupo[p, ecuacionesPoligonal[vertpol],  
   distAlOrigen[vertpol] ]
   ];
perteneceASemigrupo[{a_, b_}, ecuaciones_, distancia_] := 
  Module[{k = 1, parar = False},
   While[Norm[{a, b}/k] >= distancia \[And] ! parar,
    parar = estaEnInteriorOFrontera[{a, b}/k, ecuaciones];
    k++;
    ];
   parar
   ];
distAlOrigen[listapuntos_] := 
  Min[ Table[Norm[listapuntos[[i]]], {i, Length[listapuntos]}]];
ecuacionesPoligonal[listapuntos_] := 
  Append[Table[ecuacionPar[listapuntos[[i]], listapuntos[[i + 1]]],
    {i, 1, Length[listapuntos] - 1}], 
   ecuacionPar[listapuntos[[Length[listapuntos]]], 
    listapuntos[[1]]]];
ecuacionPar[{a1_, a2_}, {b1_, b2_}] := Det[( {
     {x - a1, y - a2},
     {b1 - a1, b2 - a2}
    } )];
estaEnInteriorOFrontera[{a_, b_}, ecuaciones_] := 
  todoNegativoOCero[ecuaciones /. {x -> a, y -> b}];
todoNegativoOCero[l_] := And @@ (# <= 0 & /@ l);
End[]
EndPackage[];

