# source:git/IntegerProgramming/IP.lib@27e935

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