1#############################################################################
2##
3#W  affine-extra-4ti2.gi
4#W                          Manuel Delgado <mdelgado@fc.up.pt>
5#W                          Pedro Garcia-Sanchez <pedro@ugr.es>
6##
7#Y  Copyright 2015-- Centro de Matemática da Universidade do Porto, Portugal and Universidad de Granada, Spain
8#############################################################################
9InstallOtherMethod(DegreesOfPrimitiveElementsOfAffineSemigroup,
10        "Computes the set of primitive elements of an affine semigroup",
11        [IsAffineSemigroup],4,
12        function(a)
13    local dir, filename, exec, filestream, matrix,
14				 facs, mat, trunc, ls;
15
16    ls:=MinimalGenerators(a);
17
18    dir := DirectoryTemporary();
19    filename := Filename( dir, "gap_4ti2_temp_matrix" );
20
21	mat:=TransposedMat(ls);
22    4ti2Interface_Write_Matrix_To_File( mat, Concatenation( filename, ".mat" ) );
23    exec := IO_FindExecutable( "graver" );
24    filestream := IO_Popen2( exec, [ filename ]);
25    while IO_ReadLine( filestream.stdout ) <> "" do od;
26    matrix := 4ti2Interface_Read_Matrix_From_File( Concatenation( filename, ".gra" ) );
27
28    trunc:=function(ls)
29		return List(ls, y->Maximum(y,0));
30	end;
31
32	matrix:=Set(matrix,trunc);
33    return Union(Set(matrix, x->x*ls),ls);
34end);
35
36
37InstallOtherMethod(HilbertBasisOfSystemOfHomogeneousEquations,
38        "Computes a Hilbert basiss of a system of linear Diophantine equations, some eventually in congruences.",
39        [IsRectangularTable,IsHomogeneousList],4,
40        function(ls,md)
41    local  homogeneous, withCongruences;
42
43    homogeneous:= function(l)
44        local  dir, filename, exec, filestream, matrix,mat,sign;
45
46        Info(InfoNumSgps,2,"Using 4ti2 for Hilbert.");
47
48        #if not(IsRectangularTable(l)) then
49        #    Error("The argument must be a matrix.");
50        #fi;
51        if not(IsInt(l[1][1])) then
52            Error("The matrix must be of integers.");
53        fi;
54
55        dir := DirectoryTemporary();
56        filename := Filename( dir, "gap_4ti2_temp_matrix" );
57
58	mat:=l;
59        sign:=[List(l[1],_->1)];
60        #Print(mat,"\n");
61        4ti2Interface_Write_Matrix_To_File( mat, Concatenation( filename, ".mat" ) );
62        4ti2Interface_Write_Matrix_To_File( sign, Concatenation( filename, ".sign" ) );
63        exec := IO_FindExecutable( "zsolve" );
64        filestream := IO_Popen2( exec, [ filename ]);
65        while IO_ReadLine( filestream.stdout ) <> "" do od;
66        matrix := 4ti2Interface_Read_Matrix_From_File( Concatenation( filename, ".zhom" ) );
67        return matrix;
68
69    end;
70
71  withCongruences:=function(ls,md)
72      local l,n,m,diag,dim,d, hil, zero, leq;
73
74      leq:= function(v1,v2)
75          local v;
76          v:=v2-v1;
77          return (First(v,n->n<0)=fail);
78      end;
79
80      #if not(IsRectangularTable(ls)) then
81      #    Error("The first argument must be a matrix.");
82      #fi;
83
84      if not(IsListOfIntegersNS(md)) or ForAny(md, x->not(IsPosInt(x))) then
85          Error("The second argument must be a list of positive integers.");
86      fi;
87
88      n:=Length(ls);
89      dim:=Length(ls[1]);
90      m:=Length(md);
91      if m>n then
92          Error("There are more modulus than equations.");
93      fi;
94
95      diag:=Concatenation(md,List([1..n-m],_->0));
96      d:=DiagonalMat(diag);
97      l:=TransposedMat(Concatenation(TransposedMat(ls),d,-d));
98      zero:=List([1..dim],_->0);
99
100      hil:=Difference(List(homogeneous(l), x->x{[1..dim]}),[zero]);
101      return hil;
102
103      return Filtered(hil, y->Filtered(hil,x->leq(x,y))=[y]);
104  end;
105  ## end of local functions ...
106
107  #ls := arg[1][1];
108  #md := arg[1][2];
109  if md = [] then
110      return homogeneous(ls);
111  else
112      return withCongruences(ls,md);
113
114  fi;
115
116end);
117
118InstallOtherMethod(HilbertBasisOfSystemOfHomogeneousInequalities,
119        "Computes a Hilbert basis of l*x>=0, x>=0",
120        [IsRectangularTable],4,
121        function(l)
122    local  dir, filename, exec, filestream, matrix,mat,sign,rel;
123
124    Info(InfoNumSgps,2,"Using 4ti2 for Hilbert.");
125
126    #if not(IsRectangularTable(l)) then
127    #    Error("The argument must be a matrix.");
128    #fi;
129    if not(IsInt(l[1][1])) then
130        Error("The matrix must be of integers.");
131    fi;
132
133    dir := DirectoryTemporary();
134    filename := Filename( dir, "gap_4ti2_temp_matrix" );
135
136    mat:=l;
137    sign:=[List(l[1],_->1)];
138    rel:=[List(l[1],_->">")];
139    #Print(mat,"\n");
140    4ti2Interface_Write_Matrix_To_File( mat, Concatenation( filename, ".mat" ) );
141    4ti2Interface_Write_Matrix_To_File( sign, Concatenation( filename, ".sign" ) );
142    4ti2Interface_Write_Matrix_To_File( rel, Concatenation( filename, ".rel" ) );
143    exec := IO_FindExecutable( "zsolve" );
144    filestream := IO_Popen2( exec, [ filename ]);
145    while IO_ReadLine( filestream.stdout ) <> "" do od;
146    matrix := 4ti2Interface_Read_Matrix_From_File( Concatenation( filename, ".zhom" ) );
147    return matrix;
148
149end);
150
151
152InstallOtherMethod(FactorizationsVectorWRTList,
153        "Computes the factorizations of v in terms of the elments in ls",
154        [IsHomogeneousList,IsMatrix],4,
155        function(v,l)
156    local  dir, filename, exec, filestream, matrix,mat,rhs,sign;
157
158    Info(InfoNumSgps,2,"Using 4ti2 for factorizations.");
159
160    if not(IsListOfIntegersNS(v)) then
161        Error("The first argument must be a list of integers.");
162    fi;
163
164    if not(IsInt(l[1][1])) then
165        Error("The matrix must be of integers.");
166    fi;
167
168    dir := DirectoryTemporary();
169    filename := Filename( dir, "gap_4ti2_temp_matrix" );
170
171    mat:=TransposedMat(l);
172    sign:=[List(mat[1],_->1)];
173    rhs:=[v];
174    #Print(mat,"\n");
175    4ti2Interface_Write_Matrix_To_File( mat, Concatenation( filename, ".mat" ) );
176    4ti2Interface_Write_Matrix_To_File( sign, Concatenation( filename, ".sign" ) );
177    4ti2Interface_Write_Matrix_To_File( rhs, Concatenation( filename, ".rhs" ) );
178    exec := IO_FindExecutable( "zsolve" );
179    filestream := IO_Popen2( exec, [ filename ]);
180    while IO_ReadLine( filestream.stdout ) <> "" do od;
181    matrix := 4ti2Interface_Read_Matrix_From_File( Concatenation( filename, ".zinhom" ) );
182    return matrix;
183
184end);
185
186
187InstallOtherMethod(GeneratorsOfKernelCongruence,
188        "Computes a set of generators of the kernel congruence of the monoid morphism associated to a matrix",
189        [IsRectangularTable],7,
190        function(m)
191    local positivenegative, gr;
192
193    positivenegative:=function(p)
194        local d1, d2;
195        d1:=List(p, i->Maximum(i,0));
196        d2:=List(p, i->-Minimum(0,i));
197        return [d1,d2];
198    end;
199
200    if not(ForAll(m, l->ForAll(l, x->(x=0) or IsPosInt(x)))) then
201        Error("The argument must be a matrix of nonnegative integer.");
202    fi;
203
204    gr:=4ti2Interface_groebner_matrix(m);
205    Info(InfoNumSgps,2,"4ti output:",gr);
206
207    return List(gr, x->positivenegative(x));
208end);
209
210
211############################################################
212# computes a canonical basis of the kernel congruence
213# of the monoid morphism associated to the matrix m with
214# nonnegative integer coefficients wrt the term ordering
215# the kernel is the pairs (x,y) such that xm=ym
216############################################################
217InstallMethod(CanonicalBasisOfKernelCongruence,
218"Computes a canonical basis for the congruence of of the monoid morphism associated to the matrix",
219	[IsRectangularTable, IsMonomialOrdering],7,
220  function(m,ord)
221    local positivenegative, gr, nord, to,dim,ones;
222
223  	positivenegative:=function(p)
224  		local d1, d2;
225  		d1:=List(p, i->Maximum(i,0));
226  		d2:=List(p, i->-Minimum(0,i));
227  		return [d1,d2];
228  	end;
229
230  	if not(ForAll(m, l->ForAll(l, x->(x=0) or IsPosInt(x)))) then
231  		Error("The argument must be a matrix of nonnegative integer.");
232  	fi;
233
234    dim:= Length(m);
235    ones:=List([1..dim],_->1);
236  	# trick taken from the package Singular
237  	nord := Name( ord );
238  	nord := nord{[ 1 .. Position( nord, '(' ) - 1 ]};
239  	if nord = "MonomialLexOrdering"  then
240  		#to := List([1..dim],_->0);
241        #to[1]:=1;
242        #Info(InfoNumSgps,1,"Warning using block ordering that discriminates the first variable wrt to the rest (4ti2Interface current release).");
243        to:=IdentityMat(dim);
244  	elif nord = "MonomialGrevlexOrdering"  then
245      to :=Concatenation([ones],Reversed(-IdentityMat(dim)){[1..dim-1]});
246    elif nord = "MonomialGrlexOrdering" then
247      to :=Concatenation([ones],IdentityMat(dim){[1..dim-1]});
248  	else
249  			Error( "the ordering ", ord, " is not yet supported in 4ti2Interface." );
250  	fi;
251
252  	gr:=4ti2Interface_groebner_matrix(m,to);
253  	Info(InfoNumSgps,2,"4ti output:",gr);
254
255  	return Set(gr, x->positivenegative(x));
256  end);
257
258
259  ############################################################
260  # computes the Graver basis of matrix with integer entries
261  ############################################################
262  InstallMethod(GraverBasis,
263          "Computes the Graver basis of the matrix",
264          [IsRectangularTable],7,
265function(a)
266  #4ti2Interface implementation
267  local gr;
268
269
270  if not(IsRectangularTable(a)) then
271    Error("The argument must be a matrix.");
272  fi;
273
274  if not(IsInt(a[1][1])) then
275    Error("The entries of the matrix must be integers.");
276  fi;
277
278  Info(InfoNumSgps,2,"Using 4ti2Interface for Graver.");
279
280  gr:=4ti2Interface_graver_equalities_in_positive_orthant(a);
281  return Union(gr,-gr);
282end);
283
284
285InstallOtherMethod(MinimalPresentationOfAffineSemigroup,
286        "Computes a minimimal presentation of the affine semigroup",
287        [IsAffineSemigroup],3,
288        function(a)
289    local gens, positive, gr, candidates, pres, rclass,exps, c;
290
291    positive:=function(x)
292        local p,i;
293
294        p:=[];
295        for i in [1..Length(x)] do
296            p[i]:=Maximum(x[i],0);
297        od;
298
299        return p;
300    end;
301    if not(IsAffineSemigroup(a)) then
302        Error("The argument must be an affine semigroup.");
303    fi;
304
305    gens:=MinimalGenerators(a);
306
307    gr:=4ti2Interface_groebner_matrix(gens);
308    Info(InfoNumSgps,2,"4ti output:",gr);
309
310    candidates:=Set(gr,q->positive(q));
311    candidates:=Set(candidates,c->c*gens);
312    Info(InfoNumSgps,2, "Candidates to Betti elements",candidates);
313    pres:=[];
314    for c in candidates do
315        exps:=FactorizationsVectorWRTList(c,gens);
316        rclass:=RClassesOfSetOfFactorizations(exps);
317        if Length(rclass)>1 then
318            pres:=Concatenation(pres,List([2..Length(rclass)],
319                          i->[rclass[1][1],rclass[i][1]]));
320        fi;
321    od;
322    return pres;
323
324
325end);
326
327
328
329#####################################################################
330# Computes the omega-primality of v in the affine semigroup a
331#####################################################################
332InstallOtherMethod(OmegaPrimalityOfElementInAffineSemigroup,
333        "Computes the omega-primality of v in the affine semigroup a",
334        [IsHomogeneousList,IsAffineSemigroup],4,
335        function(v,a)
336    local  ls, n, mat,extfact,par,tot,le;
337
338    le:=function(a,b)  #ordinary partial order
339    	return ForAll(b-a,x-> x>=0);
340    end;
341
342    if not(IsAffineSemigroup(a)) then
343        Error("The second argument must be an affine semigroup");
344    fi;
345
346    if not(IsListOfIntegersNS(v)) then
347        Error("The first argument must be a list of integers.");
348    fi;
349
350    if not(ForAll(v, x-> x>=0)) then
351        Error("The first argument must be a list of on nonnegative integers.");
352    fi;
353
354    ls:=MinimalGenerators(a);
355    n:=Length(ls);
356    mat:=TransposedMat(Concatenation(ls,-ls,[-v]));
357
358    if not(IsRectangularTable(mat)) then
359        Error("The first argument has not the dimension of the second.");
360    fi;
361
362    extfact:=FactorizationsVectorWRTList(v,Concatenation(ls,-ls));
363
364    par:=Set(extfact, f->f{[1..n]});
365    tot:=Filtered(par, f-> Filtered(par, g-> le(g,f))=[f]);
366    Info(InfoNumSgps,2,"Minimals of v+ls =",tot);
367    if tot=[] then
368        return 0;
369    fi;
370
371    return Maximum(Set(tot, Sum));
372
373
374end);
375