source: git/Singular/LIB/intprog.lib @ 93315a2

spielwiese
Last change on this file since 93315a2 was 93315a2, checked in by Hans Schoenemann <hannes@…>, 2 years ago
fix: typo
  • Property mode set to 100644
File size: 19.0 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version intprog.lib 4.2.1.3 Dec_2021 "; // $Id$
3category="Commutative Algebra";
4info="
5LIBRARY: intprog.lib      Integer Programming with Groebner Basis Methods
6AUTHOR:  Christine Theis, email: ctheis@math.uni-sb.de
7
8PROCEDURES:
9 solve_IP(..);   procedures for solving Integer Programming  problems
10";
11
12///////////////////////////////////////////////////////////////////////////////
13static proc solve_IP_1(intmat A, intvec bx, intvec c, string alg)
14{
15  intvec v;
16  // to be returned
17
18  // check arguments as far as necessary
19  // other inconsistencies are detected by the external program
20  if(size(c)!=ncols(A))
21  {
22    "ERROR: The number of matrix columns does not equal the size of the cost vector.";
23    return(v);
24  }
25
26  // create first temporary file with which the external program is
27  // called
28
29  int process=system("pid");
30  string matrixfile="temp_MATRIX"+string(process);
31  link MATRIX=":w "+matrixfile;
32  open(MATRIX);
33
34  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
35  int i,j;
36  for(j=1;j<=ncols(A);j++)
37  {
38    write(MATRIX,c[j]);
39  }
40  write(MATRIX,"rows:",nrows(A),"matrix:");
41  for(i=1;i<=nrows(A);i++)
42  {
43    for(j=1;j<=ncols(A);j++)
44    {
45      write(MATRIX,A[i,j]);
46    }
47  }
48
49  // search for positive row space vector, if required by the
50  // algorithm
51  int found=0;
52  if((alg=="blr") || (alg=="hs"))
53  {
54    for(i=1;i<=nrows(A);i++)
55    {
56      found=i;
57      for(j=1;j<=ncols(A);j++)
58      {
59        if(A[i,j]<=0)
60        {
61          found=0;
62        }
63      }
64      if(found>0)
65      {
66        break;
67      }
68    }
69    if(found==0)
70    {
71      "ERROR: The chosen algorithm needs a positive vector in the row space of the matrix.";
72      close(MATRIX);
73      system("sh","rm -f "+matrixfile);
74      return(v);
75    }
76    write(MATRIX,"positive row space vector:");
77    for(j=1;j<=ncols(A);j++)
78    {
79      write(MATRIX,A[found,j]);
80    }
81  }
82  close(MATRIX);
83
84  // create second temporary file for the external program
85
86  string problemfile="temp_PROBLEM"+string(process);
87  link PROBLEM=":w "+problemfile;
88  open(PROBLEM);
89
90  write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:");
91  for(i=1;i<=size(bx);i++)
92  {
93    write(PROBLEM,bx[i]);
94  }
95  close(PROBLEM);
96
97  // call external program
98  if (status(system("SingularBin")+"solve_IP","exists")=="yes")
99  {
100    int dummy=system("sh",system("SingularBin")+"solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
101  }
102  else
103  {
104    int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
105  }
106
107  // read solution from created file
108  link SOLUTION=":r "+matrixfile+".sol."+alg;
109  string solution=read(SOLUTION);
110  int pos;
111  string s;
112  if(alg=="ct" || alg=="pct")
113  {
114    pos=find(solution,"NO");
115    if(pos!=0)
116    {
117      "not solvable";
118    }
119    else
120    {
121      pos=find(solution,"YES");
122      pos=find(solution,":",pos);
123      pos++;
124      for(j=1;j<=ncols(A);j++)
125      {
126        while(solution[pos]==" " || solution[pos]==newline)
127        {
128          pos++;
129        }
130        s="";
131        while(solution[pos]!=" " && solution[pos]!=newline)
132        {
133          s=s+solution[pos];
134          pos++;
135        }
136        execute("v[j]="+s+";");
137      }
138    }
139  }
140  else
141  {
142    pos=find(solution,"optimal");
143    pos=find(solution,":",pos);
144    pos++;
145    for(j=1;j<=ncols(A);j++)
146    {
147      while(solution[pos]==" " || solution[pos]==newline)
148      {
149        pos++;
150      }
151      s="";
152      while(solution[pos]!=" " && solution[pos]!=newline)
153      {
154        s=s+solution[pos];
155        pos++;
156      }
157      execute("v[j]="+s+";");
158    }
159  }
160
161  // delete all created files
162  dummy=system("sh","rm -f "+matrixfile);
163  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
164  dummy=system("sh","rm -f "+problemfile);
165  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);
166
167  return(v);
168}
169///////////////////////////////////////////////////////////////////////////////
170static proc solve_IP_2(intmat A, list bx, intvec c, string alg)
171{
172  list l;;
173  // to be returned
174
175  // check arguments as far as necessary
176  // other inconsistencies are detected by the external program
177  if(size(c)!=ncols(A))
178  {
179    "ERROR: The number of matrix columns does not equal the size of the cost vector.";
180    return(l);
181  }
182
183  int k;
184  for(k=2;k<=size(bx);k++)
185  {
186    if(size(bx[k])!=size(bx[1]))
187    {
188      "ERROR: The size of all right-hand vectors must be equal.";
189      return(l);
190    }
191  }
192
193  // create first temporary file with which the external program is
194  // called
195
196  int process=system("pid");
197  string matrixfile="temp_MATRIX"+string(process);
198  link MATRIX=":w "+matrixfile;
199  open(MATRIX);
200
201  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
202  int i,j;
203  for(j=1;j<=ncols(A);j++)
204  {
205    write(MATRIX,c[j]);
206  }
207  write(MATRIX,"rows:",nrows(A),"matrix:");
208  for(i=1;i<=nrows(A);i++)
209  {
210    for(j=1;j<=ncols(A);j++)
211    {
212      write(MATRIX,A[i,j]);
213    }
214  }
215
216  // search for positive row space vector, if required by the
217  // algorithm
218  int found=0;
219  if((alg=="blr") || (alg=="hs"))
220  {
221    for(i=1;i<=nrows(A);i++)
222    {
223      found=i;
224      for(j=1;j<=ncols(A);j++)
225      {
226        if(A[i,j]<=0)
227        {
228          found=0;
229        }
230      }
231      if(found>0)
232      {
233        break;
234      }
235    }
236    if(found==0)
237    {
238      "ERROR: The chosen algorithm needs a positive vector in the row space of the matrix.";
239      close(MATRIX);
240      system("sh","rm -f "+matrixfile);
241      return(l);
242    }
243    write(MATRIX,"positive row space vector:");
244    for(j=1;j<=ncols(A);j++)
245    {
246      write(MATRIX,A[found,j]);
247    }
248  }
249  close(MATRIX);
250
251  // create second temporary file for the external program
252
253  string problemfile="temp_PROBLEM"+string(process);
254  link PROBLEM=":w "+problemfile;
255  open(PROBLEM);
256
257  write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:");
258  for(k=1;k<=size(bx);k++)
259  {
260    for(i=1;i<=size(bx[1]);i++)
261    {
262      write(PROBLEM,bx[k][i]);
263    }
264  }
265  close(PROBLEM);
266
267  // call external program
268  if (status(system("SingularBin")+"solve_IP","exists")=="yes")
269  {
270    int dummy=system("sh",system("SingularBin")+"solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
271  }
272  else
273  {
274    int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
275  }
276
277  // read solution from created file
278  link SOLUTION=":r "+matrixfile+".sol."+alg;
279  string solution=read(SOLUTION);
280  intvec v;
281  int pos,pos1,pos2;
282  string s;
283  if(alg=="ct" || alg=="pct")
284  {
285    pos=1;
286    for(k=1;k<=size(bx);k++)
287    {
288      pos1=find(solution,"NO",pos);
289      pos2=find(solution,"YES",pos);
290      if(pos1!=0 && (pos1<pos2 || pos2==0))
291        // first problem not solvable
292      {
293        pos=find(solution,":",pos1);
294        l=insert(l,"not solvable",size(l));
295      }
296      else
297        // first problem solvable
298      {
299        pos=find(solution,":",pos2);
300        pos++;
301        for(j=1;j<=ncols(A);j++)
302        {
303          while(solution[pos]==" " || solution[pos]==newline)
304          {
305            pos++;
306          }
307          s="";
308          while(solution[pos]!=" " && solution[pos]!=newline)
309          {
310            s=s+solution[pos];
311            pos++;
312          }
313          execute("v[j]="+s+";");
314        }
315        l=insert(l,v,size(l));
316      }
317    }
318  }
319  else
320  {
321    pos=1;
322    for(k=1;k<=size(bx);k++)
323    {
324      pos=find(solution,"optimal",pos);
325      pos=find(solution,":",pos);
326      pos++;
327      for(j=1;j<=ncols(A);j++)
328      {
329        while(solution[pos]==" " || solution[pos]==newline)
330        {
331          pos++;
332        }
333        s="";
334        while(solution[pos]!=" " && solution[pos]!=newline)
335        {
336          s=s+solution[pos];
337          pos++;
338        }
339        execute("v[j]="+s+";");
340      }
341      l=insert(l,v,size(l));
342    }
343  }
344
345  // delete all created files
346  dummy=system("sh","rm -f "+matrixfile);
347  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
348  dummy=system("sh","rm -f "+problemfile);
349  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);
350
351  return(l);
352}
353///////////////////////////////////////////////////////////////////////////////
354
355static proc solve_IP_3(intmat A, intvec bx, intvec c, string alg, intvec prsv)
356{
357  intvec v;
358  // to be returned
359
360  // check arguments as far as necessary
361  // other inconsistencies are detected by the external program
362  if(size(c)!=ncols(A))
363  {
364    "ERROR: The number of matrix columns does not equal the size of the cost vector.";
365    return(v);
366  }
367
368  if(size(prsv)!=ncols(A))
369  {
370    "ERROR: The number of matrix columns does not equal the size of the positive row space vector.";
371    return(v);
372  }
373
374  // create first temporary file with which the external program is
375  // called
376
377  int process=system("pid");
378  string matrixfile="temp_MATRIX"+string(process);
379  link MATRIX=":w "+matrixfile;
380  open(MATRIX);
381
382  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
383  int i,j;
384  for(j=1;j<=ncols(A);j++)
385  {
386    write(MATRIX,c[j]);
387  }
388  write(MATRIX,"rows:",nrows(A),"matrix:");
389  for(i=1;i<=nrows(A);i++)
390  {
391    for(j=1;j<=ncols(A);j++)
392    {
393      write(MATRIX,A[i,j]);
394    }
395  }
396
397  // enter positive row space vector, if required by the algorithm
398  if((alg=="blr") || (alg=="hs"))
399  {
400    write(MATRIX,"positive row space vector:");
401    for(j=1;j<=ncols(A);j++)
402    {
403      write(MATRIX,prsv[j]);
404    }
405  }
406  close(MATRIX);
407
408  // create second temporary file for the external program
409
410  string problemfile="temp_PROBLEM"+string(process);
411  link PROBLEM=":w "+problemfile;
412  open(PROBLEM);
413
414  write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:");
415  for(i=1;i<=size(bx);i++)
416  {
417    write(PROBLEM,bx[i]);
418  }
419  close(PROBLEM);
420
421  // call external program
422  if (status(system("SingularBin")+"solve_IP","exists")=="yes")
423  {
424    int dummy=system("sh",system("SingularBin")+"solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
425  }
426  else
427  {
428    int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
429  }
430
431  // read solution from created file
432  link SOLUTION=":r "+matrixfile+".sol."+alg;
433  string solution=read(SOLUTION);
434  int pos;
435  string s;
436  if(alg=="ct" || alg=="pct")
437  {
438    pos=find(solution,"NO");
439    if(pos!=0)
440    {
441      "not solvable";
442    }
443    else
444    {
445      pos=find(solution,"YES");
446      pos=find(solution,":",pos);
447      pos++;
448      for(j=1;j<=ncols(A);j++)
449      {
450        while(solution[pos]==" " || solution[pos]==newline)
451        {
452          pos++;
453        }
454        s="";
455        while(solution[pos]!=" " && solution[pos]!=newline)
456        {
457          s=s+solution[pos];
458          pos++;
459        }
460        execute("v[j]="+s+";");
461      }
462    }
463  }
464  else
465  {
466    pos=find(solution,"optimal");
467    pos=find(solution,":",pos);
468    pos++;
469    for(j=1;j<=ncols(A);j++)
470    {
471      while(solution[pos]==" " || solution[pos]==newline)
472      {
473        pos++;
474      }
475      s="";
476      while(solution[pos]!=" " && solution[pos]!=newline)
477      {
478        s=s+solution[pos];
479        pos++;
480      }
481      execute("v[j]="+s+";");
482    }
483  }
484
485  // delete all created files
486  dummy=system("sh","rm -f "+matrixfile);
487  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
488  dummy=system("sh","rm -f "+problemfile);
489  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);
490
491  return(v);
492}
493///////////////////////////////////////////////////////////////////////////////
494
495static proc solve_IP_4(intmat A, list bx, intvec c, string alg, intvec prsv)
496{
497  list l;
498  // to be returned
499
500  // check arguments as far as necessary
501  // other inconsistencies are detected by the external program
502  if(size(c)!=ncols(A))
503  {
504    "ERROR: The number of matrix columns does not equal the size of the cost vector.";
505    return(l);
506  }
507
508  if(size(prsv)!=ncols(A))
509  {
510    "ERROR: The number of matrix columns does not equal the size of the positive row space vector";
511    return(v);
512  }
513
514  int k;
515  for(k=2;k<=size(bx);k++)
516  {
517    if(size(bx[k])!=size(bx[1]))
518    {
519      "ERROR: The size of all right-hand vectors must be equal.";
520      return(l);
521    }
522  }
523
524  // create first temporary file with which the external program is
525  // called
526
527  int process=system("pid");
528  string matrixfile="temp_MATRIX"+string(process);
529  link MATRIX=":w "+matrixfile;
530  open(MATRIX);
531
532  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
533  int i,j;
534  for(j=1;j<=ncols(A);j++)
535  {
536    write(MATRIX,c[j]);
537  }
538  write(MATRIX,"rows:",nrows(A),"matrix:");
539  for(i=1;i<=nrows(A);i++)
540  {
541    for(j=1;j<=ncols(A);j++)
542    {
543      write(MATRIX,A[i,j]);
544    }
545  }
546
547  // enter positive row space vector if required by the algorithm
548  if((alg=="blr") || (alg=="hs"))
549  {
550    write(MATRIX,"positive row space vector:");
551    for(j=1;j<=ncols(A);j++)
552    {
553      write(MATRIX,prsv[j]);
554    }
555  }
556  close(MATRIX);
557
558  // create second temporary file for the external program
559
560  string problemfile="temp_PROBLEM"+string(process);
561  link PROBLEM=":w "+problemfile;
562  open(PROBLEM);
563
564  write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:");
565  for(k=1;k<=size(bx);k++)
566  {
567    for(i=1;i<=size(bx[1]);i++)
568    {
569      write(PROBLEM,bx[k][i]);
570    }
571  }
572  close(PROBLEM);
573
574  // call external program
575  if (status(system("SingularBin")+"solve_IP","exists")=="yes")
576  {
577    int dummy=system("sh",system("SingularBin")+"solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
578  }
579  else
580  {
581    int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);
582  }
583
584  // read solution from created file
585  link SOLUTION=":r "+matrixfile+".sol."+alg;
586  string solution=read(SOLUTION);
587  intvec v;
588  int pos,pos1,pos2;
589  string s;
590  if(alg=="ct" || alg=="pct")
591  {
592    pos=1;
593    for(k=1;k<=size(bx);k++)
594    {
595      pos1=find(solution,"NO",pos);
596      pos2=find(solution,"YES",pos);
597      if(pos1!=0 && (pos1<pos2 || pos2==0))
598        // first problem not solvable
599      {
600        pos=find(solution,":",pos1);
601        l=insert(l,"not solvable",size(l));
602      }
603      else
604        // first problem solvable
605      {
606        pos=find(solution,":",pos2);
607        pos++;
608        for(j=1;j<=ncols(A);j++)
609        {
610          while(solution[pos]==" " || solution[pos]==newline)
611          {
612            pos++;
613          }
614          s="";
615          while(solution[pos]!=" " && solution[pos]!=newline)
616          {
617            s=s+solution[pos];
618            pos++;
619          }
620          execute("v[j]="+s+";");
621        }
622        l=insert(l,v,size(l));
623      }
624    }
625  }
626  else
627  {
628    pos=1;
629    for(k=1;k<=size(bx);k++)
630    {
631      pos=find(solution,"optimal",pos);
632      pos=find(solution,":",pos);
633      pos++;
634      for(j=1;j<=ncols(A);j++)
635      {
636        while(solution[pos]==" " || solution[pos]==newline)
637        {
638          pos++;
639        }
640        s="";
641        while(solution[pos]!=" " && solution[pos]!=newline)
642        {
643          s=s+solution[pos];
644          pos++;
645        }
646        execute("v[j]="+s+";");
647      }
648      l=insert(l,v,size(l));
649    }
650  }
651
652  // delete all created files
653  dummy=system("sh","rm -f "+matrixfile);
654  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
655  dummy=system("sh","rm -f "+problemfile);
656  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);
657
658  return(l);
659}
660///////////////////////////////////////////////////////////////////////////////
661
662proc solve_IP
663"USAGE:  solve_IP(A,bx,c,alg);  A intmat, bx intvec, c intvec, alg string.@*
664        solve_IP(A,bx,c,alg);  A intmat, bx list of intvec, c intvec,
665                                 alg string.@*
666        solve_IP(A,bx,c,alg,prsv);  A intmat, bx intvec, c intvec,
667                                 alg string, prsv intvec.@*
668        solve_IP(A,bx,c,alg,prsv);  A intmat, bx list of intvec, c intvec,
669                                 alg string, prsv intvec.
670RETURN: same type as bx: solution of the associated integer programming
671        problem(s) as explained in
672   @texinfo
673   @ref{Toric ideals and integer programming}.
674   @end texinfo
675NOTE:   This procedure returns the solution(s) of the given IP-problem(s)
676        or the message `not solvable'.
677        One may call the procedure with several different algorithms:
678        @*- the algorithm of Conti/Traverso (ct),
679        @*- the positive variant of the algorithm of Conti/Traverso (pct),
680        @*- the algorithm of Conti/Traverso using elimination (ect),
681        @*- the algorithm of Pottier (pt),
682        @*- an algorithm of Bigatti/La Scala/Robbiano (blr),
683        @*- the algorithm of Hosten/Sturmfels (hs),
684        @*- the algorithm of DiBiase/Urbanke (du).
685        @*The argument `alg' should be the abbreviation for an algorithm as
686        above: ct, pct, ect, pt, blr, hs or du.
687
688        `ct' allows computation of an optimal solution of the IP-problem
689        directly from the right-hand vector b.
690        The same is true for its `positive' variant `pct' which may only be
691        applied if A and b have nonnegative entries.
692        All other algorithms need initial solutions of the IP-problem.
693
694        If `alg' is chosen to be `ct' or `pct', bx is read as the right hand
695        vector b of the system Ax=b. b should then be an intvec of size m
696        where m is the number of rows of A.
697        Furthermore, bx and A should be nonnegative if `pct' is used.
698        If `alg' is chosen to be `ect',`pt',`blr',`hs' or `du',
699        bx is read as an initial solution x of the system Ax=b.
700        bx should then be a nonnegative intvec of size n where n is the
701        number of columns of A.
702
703        If `alg' is chosen to be `blr' or `hs', the algorithm needs a vector
704        with positive coefficients in the row space of A.
705        If no row of A contains only positive entries, one has to use the
706        versions of solve_IP which take such a vector prsv as an argument.
707
708        solve_IP may also be called with a list bx of intvecs instead of a
709        single intvec.
710SEE ALSO: intprog_lib, toric_lib, Integer programming
711EXAMPLE:  example solve_IP; shows an example
712"
713{
714  if(size(#)==4)
715  {
716    if(typeof(#[2])=="intvec")
717    {
718      return(solve_IP_1(#[1],#[2],#[3],#[4]));
719    }
720    else
721    {
722      return(solve_IP_2(#[1],#[2],#[3],#[4]));
723    }
724  }
725  else
726  {
727    if(typeof(#[2])=="intvec")
728    {
729      return(solve_IP_3(#[1],#[2],#[3],#[4],#[5]));
730    }
731    else
732    {
733      return(solve_IP_4(#[1],#[2],#[3],#[4],#[5]));
734    }
735  }
736}
737
738
739
740example
741{ "EXAMPLE"; echo=2;
742  // 1. call with single right-hand vector
743  intmat A[2][3]=1,1,0,0,1,1;
744  intvec b1=1,1;
745  intvec c=2,2,1;
746  intvec solution_vector=solve_IP(A,b1,c,"pct");
747  solution_vector;"";
748
749  // 2. call with list of right-hand vectors
750  intvec b2=-1,1;
751  list l=b1,b2;
752  l;
753  list solution_list=solve_IP(A,l,c,"ct");
754  solution_list;"";
755
756  // 3. call with single initial solution vector
757  A=2,1,-1,-1,1,2;
758  b1=3,4,5;
759  solve_IP(A,b1,c,"du");"";
760
761  // 4. call with single initial solution vector
762  //    and algorithm needing a positive row space vector
763  solution_vector=solve_IP(A,b1,c,"hs");"";
764
765  // 5. call with single initial solution vector
766  //     and positive row space vector
767  intvec prsv=1,2,1;
768  solution_vector=solve_IP(A,b1,c,"hs",prsv);
769  solution_vector;"";
770
771  // 6. call with list of initial solution vectors
772  //    and positive row space vector
773  b2=7,8,0;
774  l=b1,b2;
775  l;
776  solution_list=solve_IP(A,l,c,"blr",prsv);
777  solution_list;
778}
Note: See TracBrowser for help on using the repository browser.