# source:git/Singular/LIB/intprog.lib@0ae4ce

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