# source:git/IntegerProgramming/IP.lib@6ba162

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