source: git/Singular/LIB/intprog.lib @ fa8122

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