Changeset 27e935 in git for IntegerProgramming
 Timestamp:
 May 11, 2000, 11:04:05 AM (23 years ago)
 Branches:
 (u'jengelhdatetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
 Children:
 54b3568dedfbe0f7886ed616396038a36933c3e3
 Parents:
 0658d94b02b071f2a019543baea28fe2c95d9226
 Location:
 IntegerProgramming
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

IntegerProgramming/IP.lib
r0658d9 r27e935 1 // $Id: IP.lib,v 1. 1.1.1 20000502 12:58:31Singular Exp $1 // $Id: IP.lib,v 1.2 20000511 09:04:05 Singular Exp $ 2 2 // 3 3 // author : Christine Theis 4 4 // 5 5 6 //version="$Id: IP.lib,v 1. 1.1.1 20000502 12:58:31Singular Exp $";6 //version="$Id: IP.lib,v 1.2 20000511 09:04:05 Singular Exp $"; 7 7 8 8 /////////////////////////////////////////////////////////////////////////////// 9 9 10 10 info=" 11 LIBRARY: IP.lib PROCEDURES FOR INTEGER PROGRAMMING USING GROEBNER BASIS METHODS11 LIBRARY: IP.lib PROCEDURES FOR INTEGER PROGRAMMING USING GROEBNER BASIS METHODS 12 12 13 13 14 14 Let A an integral (mxn)matrix, b a vector with m integral 15 15 coefficients and c a vector with n nonnegative integral coefficients. The 16 solution of the IPproblem 17 16 solution of the IPproblem 17 18 18 (*) minimize{ cx  Ax=b, all components of x are nonnegative integers } 19 19 20 20 proceeds in two steps: 21 21 22 22 1. We compute the toric ideal of A and its Groebner basis with respect 23 23 to a term ordering refining the cost function c (such an ordering 24 24 exists because c is nonnegative). 25 25 26 26 2. We reduce the right hand vector b or an initial solution of the problem 27 27 modulo this ideal. 28 28 29 29 For these purposes, we can use seven different algorithms: 30 30 The algorithm of Conti/Traverso (ct) can compute an optimal solution 31 31 of the IPproblem directly from the right hand vector b. The same is 32 32 true for its `positive' variant (pct) which can only be applied if A 33 and b have nonnegative entries, but should be faster in that case. 33 and b have nonnegative entries, but should be faster in that case. 34 34 All other algorithms need initial solutions of the IPproblem. Except 35 35 from the ContiTraverso algorithm with elimination (ect), … … 40 40 41 41 42 solve_IP(intmat A, intvec bx, intvec c, string alg); 42 solve_IP(intmat A, intvec bx, intvec c, string alg); 43 43 solve_IP(intmat A, list bx, intvec c, string alg); 44 solve_IP(intmat A, intvec bx, intvec c, string alg, intvec prsv); 44 solve_IP(intmat A, intvec bx, intvec c, string alg, intvec prsv); 45 45 solve_IP(intmat A, list bx, intvec c, string alg, intvec prsv); 46 47 procedures for solving the IPproblem (*) 48 They return the solution(s) of the given problem(s) or the 46 47 procedures for solving the IPproblem (*) 48 They return the solution(s) of the given problem(s) or the 49 49 message `not solvable'. 50 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 51 If `alg' is chosen to be `ct' or `pct', bx is read as the right 52 52 hand vector b of the system Ax=b. b should then be an intvec of 53 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 56 57 58 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 59 If `alg' is chosen to be `blr' or `hs', the algorithm needs a 60 60 vector with positive coefficcients in the row space of A. If 61 61 no row of A contains only positive entries, one must use the 62 62 versions of solve_IP which take such a vector prsv as 63 argument. 64 65 of a single intvec. 66 67 63 argument. 64 solve_IP may also be called with a list bx of intvecs instead 65 of a single intvec. 66 67 68 68 "; 69 69 70 70 /////////////////////////////////////////////////////////////////////////////// 71 71 72 73 74 72 static proc solve_IP_1(intmat A, intvec bx, intvec c, string alg) 75 { 73 { 76 74 intvec v; 77 75 // to be returned 78 76 79 77 // check arguments as far as necessary 80 // other inconsistencies are detected by the external program 78 // other inconsistencies are detected by the external program 81 79 if(size(c)!=ncols(A)) 82 80 { 83 81 "ERROR: number of matrix columns must equal size of cost vector"; 84 82 return(v); 85 } 83 } 86 84 87 85 // create first temporary file with that the external program is 88 // called 86 // called 89 87 90 88 int process=system("pid"); 91 string matrixfile="temp_MATRIX"+string(process); 89 string matrixfile="temp_MATRIX"+string(process); 92 90 link MATRIX=":w "+matrixfile; 93 91 open(MATRIX); … … 107 105 } 108 106 } 109 107 110 108 // search for positive row space vector, if required by the 111 // algorithm 109 // algorithm 112 110 int found=0; 113 111 if((alg=="blr")  (alg=="hs")) 114 { 112 { 115 113 for(i=1;i<=nrows(A);i++) 116 114 { … … 130 128 if(found==0) 131 129 { 132 "ERROR: algorithm needs positive vector in the row space of the matrix"; 130 "ERROR: algorithm needs positive vector in the row space of the matrix"; 133 131 close(MATRIX); 134 132 system("sh","rm f "+matrixfile); 135 133 return(v); 136 } 134 } 137 135 write(MATRIX,"positive row space vector:"); 138 136 for(j=1;j<=ncols(A);j++) 139 { 137 { 140 138 write(MATRIX,A[found,j]); 141 139 } … … 145 143 // create second temporary file for the external program 146 144 147 string problemfile="temp_PROBLEM"+string(process); 145 string problemfile="temp_PROBLEM"+string(process); 148 146 link PROBLEM=":w "+problemfile; 149 147 open(PROBLEM); 150 148 151 write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:"); 149 write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:"); 152 150 for(i=1;i<=size(bx);i++) 153 151 { … … 165 163 string s; 166 164 if(alg=="ct"  alg=="pct") 167 { 165 { 168 166 pos=find(solution,"NO"); 169 167 if(pos!=0) … … 193 191 } 194 192 else 195 { 193 { 196 194 pos=find(solution,"optimal"); 197 195 pos=find(solution,":",pos); … … 212 210 } 213 211 } 214 212 215 213 // delete all created files 216 214 dummy=system("sh","rm f "+matrixfile); … … 222 220 } 223 221 224 225 226 222 static proc solve_IP_2(intmat A, list bx, intvec c, string alg) 227 { 223 { 228 224 list l;; 229 225 // to be returned 230 226 231 227 // check arguments as far as necessary 232 // other inconsistencies are detected by the external program 228 // other inconsistencies are detected by the external program 233 229 if(size(c)!=ncols(A)) 234 230 { 235 231 "ERROR: number of matrix columns must equal size of cost vector"; 236 232 return(l); 237 } 233 } 238 234 239 235 int k; … … 248 244 249 245 // create first temporary file with that the external program is 250 // called 246 // called 251 247 252 248 int process=system("pid"); 253 string matrixfile="temp_MATRIX"+string(process); 249 string matrixfile="temp_MATRIX"+string(process); 254 250 link MATRIX=":w "+matrixfile; 255 251 open(MATRIX); … … 269 265 } 270 266 } 271 267 272 268 // search for positive row space vector, if required by the 273 // algorithm 269 // algorithm 274 270 int found=0; 275 271 if((alg=="blr")  (alg=="hs")) 276 { 272 { 277 273 for(i=1;i<=nrows(A);i++) 278 274 { … … 292 288 if(found==0) 293 289 { 294 "ERROR: algorithm needs positive vector in the row space of the matrix"; 290 "ERROR: algorithm needs positive vector in the row space of the matrix"; 295 291 close(MATRIX); 296 292 system("sh","rm f "+matrixfile); 297 293 return(l); 298 } 294 } 299 295 write(MATRIX,"positive row space vector:"); 300 296 for(j=1;j<=ncols(A);j++) 301 { 297 { 302 298 write(MATRIX,A[found,j]); 303 299 } … … 307 303 // create second temporary file for the external program 308 304 309 string problemfile="temp_PROBLEM"+string(process); 305 string problemfile="temp_PROBLEM"+string(process); 310 306 link PROBLEM=":w "+problemfile; 311 307 open(PROBLEM); 312 308 313 write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:"); 309 write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:"); 314 310 for(k=1;k<=size(bx);k++) 315 311 { … … 331 327 string s; 332 328 if(alg=="ct"  alg=="pct") 333 { 329 { 334 330 pos=1; 335 331 for(k=1;k<=size(bx);k++) 336 { 332 { 337 333 pos1=find(solution,"NO",pos); 338 334 pos2=find(solution,"YES",pos); … … 359 355 s=s+solution[pos]; 360 356 pos++; 361 } 357 } 362 358 execute("v[j]="+s+";"); 363 359 } … … 367 363 } 368 364 else 369 { 365 { 370 366 pos=1; 371 367 for(k=1;k<=size(bx);k++) … … 385 381 s=s+solution[pos]; 386 382 pos++; 387 } 383 } 388 384 execute("v[j]="+s+";"); 389 385 } … … 391 387 } 392 388 } 393 389 394 390 // delete all created files 395 391 dummy=system("sh","rm f "+matrixfile); … … 401 397 } 402 398 403 404 405 399 static proc solve_IP_3(intmat A, intvec bx, intvec c, string alg, intvec prsv) 406 { 400 { 407 401 intvec v; 408 402 // to be returned 409 403 410 404 // check arguments as far as necessary 411 // other inconsistencies are detected by the external program 405 // other inconsistencies are detected by the external program 412 406 if(size(c)!=ncols(A)) 413 407 { 414 408 "ERROR: number of matrix columns must equal size of cost vector"; 415 409 return(v); 416 } 410 } 417 411 418 412 if(size(prsv)!=ncols(A)) … … 420 414 "ERROR: number of matrix columns must equal size of positive row space vector"; 421 415 return(v); 422 } 416 } 423 417 424 418 // create first temporary file with that the external program is 425 // called 419 // called 426 420 427 421 int process=system("pid"); 428 string matrixfile="temp_MATRIX"+string(process); 422 string matrixfile="temp_MATRIX"+string(process); 429 423 link MATRIX=":w "+matrixfile; 430 424 open(MATRIX); … … 444 438 } 445 439 } 446 440 447 441 // enter positive row space vector, if required by the algorithm 448 442 if((alg=="blr")  (alg=="hs")) 449 { 443 { 450 444 write(MATRIX,"positive row space vector:"); 451 445 for(j=1;j<=ncols(A);j++) 452 { 446 { 453 447 write(MATRIX,prsv[j]); 454 448 } … … 458 452 // create second temporary file for the external program 459 453 460 string problemfile="temp_PROBLEM"+string(process); 454 string problemfile="temp_PROBLEM"+string(process); 461 455 link PROBLEM=":w "+problemfile; 462 456 open(PROBLEM); 463 457 464 write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:"); 458 write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:"); 465 459 for(i=1;i<=size(bx);i++) 466 460 { … … 478 472 string s; 479 473 if(alg=="ct"  alg=="pct") 480 { 474 { 481 475 pos=find(solution,"NO"); 482 476 if(pos!=0) … … 500 494 s=s+solution[pos]; 501 495 pos++; 502 } 496 } 503 497 execute("v[j]="+s+";"); 504 498 } … … 506 500 } 507 501 else 508 { 502 { 509 503 pos=find(solution,"optimal"); 510 504 pos=find(solution,":",pos); … … 521 515 s=s+solution[pos]; 522 516 pos++; 523 } 517 } 524 518 execute("v[j]="+s+";"); 525 519 } 526 520 } 527 521 528 522 // delete all created files 529 523 dummy=system("sh","rm f "+matrixfile); … … 535 529 } 536 530 537 538 539 531 static proc solve_IP_4(intmat A, list bx, intvec c, string alg, intvec prsv) 540 { 532 { 541 533 list l; 542 534 // to be returned 543 535 544 536 // check arguments as far as necessary 545 // other inconsistencies are detected by the external program 537 // other inconsistencies are detected by the external program 546 538 if(size(c)!=ncols(A)) 547 539 { 548 540 "ERROR: number of matrix columns must equal size of cost vector"; 549 541 return(l); 550 } 542 } 551 543 552 544 if(size(prsv)!=ncols(A)) … … 554 546 "ERROR: number of matrix columns must equal size of positive row space vector"; 555 547 return(v); 556 } 548 } 557 549 558 550 int k; … … 567 559 568 560 // create first temporary file with that the external program is 569 // called 561 // called 570 562 571 563 int process=system("pid"); 572 string matrixfile="temp_MATRIX"+string(process); 564 string matrixfile="temp_MATRIX"+string(process); 573 565 link MATRIX=":w "+matrixfile; 574 566 open(MATRIX); … … 588 580 } 589 581 } 590 582 591 583 // enter positive row space vector if required by the algorithm 592 584 if((alg=="blr")  (alg=="hs")) 593 { 585 { 594 586 write(MATRIX,"positive row space vector:"); 595 587 for(j=1;j<=ncols(A);j++) 596 { 588 { 597 589 write(MATRIX,prsv[j]); 598 590 } … … 602 594 // create second temporary file for the external program 603 595 604 string problemfile="temp_PROBLEM"+string(process); 596 string problemfile="temp_PROBLEM"+string(process); 605 597 link PROBLEM=":w "+problemfile; 606 598 open(PROBLEM); 607 599 608 write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:"); 600 write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:"); 609 601 for(k=1;k<=size(bx);k++) 610 602 { … … 626 618 string s; 627 619 if(alg=="ct"  alg=="pct") 628 { 620 { 629 621 pos=1; 630 622 for(k=1;k<=size(bx);k++) 631 { 623 { 632 624 pos1=find(solution,"NO",pos); 633 625 pos2=find(solution,"YES",pos); … … 654 646 s=s+solution[pos]; 655 647 pos++; 656 } 648 } 657 649 execute("v[j]="+s+";"); 658 650 } … … 662 654 } 663 655 else 664 { 656 { 665 657 pos=1; 666 658 for(k=1;k<=size(bx);k++) … … 686 678 } 687 679 } 688 680 689 681 // delete all created files 690 682 dummy=system("sh","rm f "+matrixfile); … … 696 688 } 697 689 698 699 700 701 702 703 690 proc solve_IP 704 "USAGE: 691 "USAGE: 705 692 solve_IP(A,bx,c,alg); A intmat, bx intvec, c intvec, alg string 706 solve_IP(A,bx,c,alg); A intmat, bx list of intvec, c intvec, alg string 693 solve_IP(A,bx,c,alg); A intmat, bx list of intvec, c intvec, alg string 707 694 solve_IP(A,bx,c,alg,prsv); A intmat, bx intvec, c intvec, alg string, 708 prsv intvec 695 prsv intvec 709 696 solve_IP(A,bx,c,alg,prsv); A intmat, bx list of intvec, c intvec, alg string, 710 prsv intvec 711 RETURN: solution of the associated integer programming problem as explained 697 prsv intvec 698 RETURN: solution of the associated integer programming problem as explained 712 699 in IP.lib 713 return type = type of bx 700 return type = type of bx 714 701 EXAMPLE: example solve_IP; shows an example" 715 702 { … … 737 724 } 738 725 } 739 740 741 742 726 example 743 727 { 744 "EXAMPLE"; echo=2;745 746 // call with single right hand vector747 intmat A[2][3]=1,1,0,0,1,1;748 A;749 intvec b1=1,1;750 b1;751 intvec c=2,2,1;752 c;753 intvec solution_vector=solve_IP(A,b1,c,"pct");754 solution_vector;755 756 // call with list of right hand vectors757 intvec b2=1,1;758 list l=b1,b2;759 l;760 list solution_list=solve_IP(A,l,c,"ct");761 solution_list; 762 763 // call with single initial solution vector764 A=2,1,1,1,1,2;765 A;766 b1=3,4,5; 767 solution_vector=solve_IP(A,b1,c,"du");768 769 // call with single initial solution vector and algorithm needing a positive770 // row space vector771 solution_vector=solve_IP(A,b1,c,"hs");772 773 // call with single initial solution vector and positive row space vector774 intvec prsv=1,2,1;775 prsv;776 solution_vector=solve_IP(A,b1,c,"hs",prsv);777 solution_vector;778 779 // call with list of initial solution vectors and positive row space vector780 b2=7,8,0;781 l=b1,b2;782 l;783 solution_list=solve_IP(A,l,c,"blr",prsv);784 solution_list;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; 785 769 } 786 787 788 
IntegerProgramming/toric.lib
r0658d9 r27e935 1 // $Id: toric.lib,v 1. 1.1.1 20000502 12:58:31 Singular Exp $1 // $Id: toric.lib,v 1.2 20000511 09:03:41 Singular Exp $ 2 2 // 3 3 // author : Christine Theis 4 4 // 5 6 //version="$Id: toric.lib,v 1.1.1.1 20000502 12:58:31 Singular Exp $"; 5 //version="$Id: toric.lib,v 1.2 20000511 09:03:41 Singular Exp $"; 7 6 8 7 /////////////////////////////////////////////////////////////////////////////// 9 8 10 9 info=" 11 LIBRARY: toric.lib PROCEDURES FOR COMPUTING TORIC IDEALS12 13 14 Let A an integral (mxn)matrix. The toric ideal I_A of A is defined 10 LIBRARY: toric.lib PROCEDURES FOR COMPUTING TORIC IDEALS 11 12 13 Let A an integral (mxn)matrix. The toric ideal I_A of A is defined 15 14 as the ideal 16 15 … … 18 17 19 18 in the ring of n variables x:=x1,...,xn. 20 Toric ideals play an important role in polyhedral geometry and may also be 21 used for integer programming. They are generated by binomials with 22 relatively prime monomials. Buchberger's algorithm can be specialized to 23 these structures in a way that considerably speeds up computation. 19 Toric ideals play an important role in polyhedral geometry and may also be 20 used for integer programming. They are generated by binomials with 21 relatively prime monomials. Buchberger's algorithm can be specialized to 22 these structures in a way that considerably speeds up computation. 24 23 25 24 26 25 toric_ideal(intmat A, string alg); 27 26 toric_ideal(intmat A, string alg, intvec prsv); 28 27 29 28 procedures for computing the toric ideal of A 30 They return the standard basis of the toric ideal of A with respect 29 They return the standard basis of the toric ideal of A with respect 31 30 to the term ordering in the actual basering. 32 When calling this procedure, a ring with n variables should be active. 31 When calling this procedure, a ring with n variables should be active. 33 32 Not all term orderings are supported: The usual global term orderings 34 33 may be used, but no block orderings combining them. … … 42 41 43 42 The last two seem to be the fastest in the actual implementation. 44 `alg' should be the abbreviation (in brackets) for an algorithm 45 as above. 43 `alg' should be the abbreviation (in brackets) for an algorithm 44 as above. 46 45 If `alg' is chosen to be `blr' or `hs', the algorithm needs a 47 46 vector with positive coefficcients in the row space of A. If 48 47 no row of A contains only positive entries, one must use the 49 second version of toric_ideal which takes such a vector in the 50 third argument. 51 48 second version of toric_ideal which takes such a vector in the 49 third argument. 50 52 51 53 52 toric_std(ideal I); 54 53 55 computes the standard basis of I using the specialized Buchberger 56 algorithm. The generating system by which I is given has to consist 54 computes the standard basis of I using the specialized Buchberger 55 algorithm. The generating system by which I is given has to consist 57 56 of binomials of the form x^ux^v (although there are other generating 58 57 systems of toric ideals). There is no real check if I is toric. 59 If the generator list of I contains a binomial whose monomials are not 58 If the generator list of I contains a binomial whose monomials are not 60 59 relatively prime, the procedure outputs a warning. If I is generated by 61 binomials of the above form, but not toric, toric_std computes an ideal 60 binomials of the above form, but not toric, toric_std computes an ideal 62 61 `between' I and its saturation with respect to all variables. 63 62 64 63 "; 65 64 66 65 /////////////////////////////////////////////////////////////////////////////// 67 66 68 69 70 67 static proc toric_ideal_1(intmat A, string alg) 71 { 68 { 72 69 ideal I; 73 70 // to be returned … … 78 75 "ERROR: number of matrix columns must be greater or equal number of ring variables"; 79 76 return(I); 80 } 77 } 81 78 82 79 // check suitability of actual term ordering … … 167 164 "Warning: block orderings are not supported; Dp used for computation"; 168 165 } 169 } 166 } 170 167 171 168 int pos; … … 263 260 return(I); 264 261 } 265 262 266 263 // check algorithm 267 264 if(alg=="ct"  alg=="pct") … … 274 271 } 275 272 276 // create temporary file with that the external program is called 273 // create temporary file with that the external program is called 277 274 278 275 int dummy; 279 276 int process=system("pid"); 280 string matrixfile="temp_MATRIX"+string(process); 277 string matrixfile="temp_MATRIX"+string(process); 281 278 link MATRIX=":w "+matrixfile; 282 279 open(MATRIX); … … 295 292 } 296 293 } 297 294 298 295 // search for positive row space vector, if required by the 299 // algorithm 296 // algorithm 300 297 int found=0; 301 298 if((alg=="blr")  (alg=="hs")) 302 { 299 { 303 300 for(i=1;i<=nrows(A);i++) 304 301 { … … 318 315 if(found==0) 319 316 { 320 "ERROR: algorithm needs positive vector in the row space of the matrix"; 317 "ERROR: algorithm needs positive vector in the row space of the matrix"; 321 318 close(MATRIX); 322 319 dummy=system("sh","rm f "+matrixfile); 323 320 return(I); 324 } 321 } 325 322 write(MATRIX,"positive row space vector:"); 326 323 for(j=1;j<=ncols(A);j++) 327 { 324 { 328 325 write(MATRIX,A[found,j]); 329 326 } … … 343 340 pos=find(toric_ideal,":",pos); 344 341 pos++; 345 342 346 343 while(toric_ideal[pos]==" "  toric_ideal[pos]==newline) 347 344 { … … 355 352 } 356 353 execute("generators="+number_string+";"); 357 354 358 355 intvec v; 359 356 poly head; … … 388 385 if(v[j]>0) 389 386 { 390 head=head*var(j)^v[j]; 387 head=head*var(j)^v[j]; 391 388 } 392 389 } 393 390 I[i]=headtail; 394 391 } 395 392 396 393 // delete all created files 397 394 dummy=system("sh","rm f "+matrixfile); … … 401 398 } 402 399 403 404 405 400 static proc toric_ideal_2(intmat A, string alg, intvec prsv) 406 { 401 { 407 402 ideal I; 408 403 // to be returned … … 413 408 "ERROR: number of matrix columns must equal size of positive row space vector"; 414 409 return(I); 415 } 410 } 416 411 417 412 // check suitability of actual basering … … 420 415 "ERROR: number of matrix columns must be greater or equal number of ring variables"; 421 416 return(I); 422 } 417 } 423 418 424 419 // check suitability of actual term ordering … … 509 504 "Warning: block orderings are not supported; Dp used for computation"; 510 505 } 511 } 506 } 512 507 513 508 int pos; … … 605 600 return(I); 606 601 } 607 602 608 603 // check algorithm 609 604 if(alg=="ct"  alg=="pct") … … 616 611 } 617 612 618 // create temporary file with that the external program is called 613 // create temporary file with that the external program is called 619 614 620 615 int dummy; 621 616 int process=system("pid"); 622 string matrixfile="temp_MATRIX"+string(process); 617 string matrixfile="temp_MATRIX"+string(process); 623 618 link MATRIX=":w "+matrixfile; 624 619 open(MATRIX); … … 640 635 // enter positive row space vector, if required by the algorithm 641 636 if((alg=="blr")  (alg=="hs")) 642 { 637 { 643 638 write(MATRIX,"positive row space vector:"); 644 639 for(j=1;j<=ncols(A);j++) 645 { 640 { 646 641 write(MATRIX,prsv[j]); 647 642 } … … 660 655 pos=find(toric_ideal,":",pos); 661 656 pos++; 662 657 663 658 while(toric_ideal[pos]==" "  toric_ideal[pos]==newline) 664 659 { … … 672 667 } 673 668 execute("generators="+number_string+";"); 674 669 675 670 intvec v; 676 671 poly head; … … 705 700 if(v[j]>0) 706 701 { 707 head=head*var(j)^v[j]; 702 head=head*var(j)^v[j]; 708 703 } 709 704 } 710 705 I[i]=headtail; 711 706 } 712 707 713 708 // delete all created files 714 709 dummy=system("sh","rm f "+matrixfile); … … 718 713 } 719 714 720 721 722 715 proc toric_ideal 723 "USAGE: 716 "USAGE: 724 717 toric_ideal(A,alg); A intmat, alg string 725 toric_ideal(A,alg,prsv); A intmat, alg string, prsv intvec 718 toric_ideal(A,alg,prsv); A intmat, alg string, prsv intvec 726 719 RETURN: toric ideal of A as explained in toric.lib 727 720 return type = ideal … … 737 730 } 738 731 } 739 740 741 742 732 example 743 733 { 744 "EXAMPLE"; echo=2;745 746 ring r=0,(x,y,z),dp; 747 748 // call with two arguments749 intmat A[2][3]=1,1,0,0,1,1;750 A;751 752 ideal I=toric_ideal(A,"du");753 I;754 755 I=toric_ideal(A,"blr");756 I;757 758 // call with three arguments759 intvec prsv=1,2,1;760 I=toric_ideal(A,"blr",prsv);761 I;762 734 "EXAMPLE"; echo=2; 735 736 ring r=0,(x,y,z),dp; 737 738 // call with two arguments 739 intmat A[2][3]=1,1,0,0,1,1; 740 A; 741 742 ideal I=toric_ideal(A,"du"); 743 I; 744 745 I=toric_ideal(A,"blr"); 746 I; 747 748 // call with three arguments 749 intvec prsv=1,2,1; 750 I=toric_ideal(A,"blr",prsv); 751 I; 752 763 753 } 764 754 765 766 767 755 proc toric_std(ideal I) 768 "USAGE: toric_std(I); I ideal 756 "USAGE: toric_std(I); I ideal 769 757 RETURN: standard basis of I as explained in toric.lib 770 758 return type = ideal … … 861 849 "Warning: block orderings are not supported; Dp used for computation"; 862 850 } 863 } 851 } 864 852 865 853 int pos; … … 887 875 } 888 876 } 889 877 890 878 if(external_ord=="" && find(singular_ord,"wp")==3) 891 879 { … … 958 946 return(I); 959 947 } 960 948 961 949 // create first temporary file with which the external program is called 962 950 963 951 int dummy; 964 952 int process=system("pid"); 965 string groebnerfile="temp_GROEBNER"+string(process); 953 string groebnerfile="temp_GROEBNER"+string(process); 966 954 link GROEBNER=":w "+groebnerfile; 967 955 open(GROEBNER); … … 969 957 write(GROEBNER,"GROEBNER","computed with algorithm:","pt","term ordering:","elimination block",0,"weighted block",nvars(basering),external_ord); 970 958 // algorithm is totally unimportant, only required by the external program 971 959 972 960 for(i=1;i<=nvars(basering);i++) 973 961 { 974 962 write(GROEBNER,weightvec[i]); 975 963 } 976 964 977 965 write(GROEBNER,"size:",size(I),"Groebner basis:"); 978 966 poly head; … … 983 971 for(j=1;j<=size(I);j++) 984 972 { 985 // test suitability of generator j 973 // test suitability of generator j 986 974 rest=I[j]; 987 975 head=lead(rest); … … 989 977 tail=lead(rest); 990 978 rest=resttail; 991 979 992 980 if(head==0 && tail==0 && rest!=0) 993 981 { … … 997 985 return(J); 998 986 } 999 987 1000 988 if(leadcoef(tail)!=leadcoef(head)) 1001 989 // generator no difference of monomials (or a constant multiple) … … 1006 994 return(J); 1007 995 } 1008 996 1009 997 if(gcd(head,tail)!=1) 1010 998 { 1011 999 "Warning: monomials of generator "+string(j)+" of input ideal are not relatively prime"; 1012 1000 } 1013 1001 1014 1002 // write vector representation of generator j into the file 1015 1003 v=leadexp(head)leadexp(tail); … … 1020 1008 } 1021 1009 close(GROEBNER); 1022 1010 1023 1011 // create second temporary file 1024 1012 1025 string newcostfile="temp_NEW_COST"+string(process); 1013 string newcostfile="temp_NEW_COST"+string(process); 1026 1014 link NEW_COST=":w "+newcostfile; 1027 1015 open(NEW_COST); 1028 1016 1029 write(NEW_COST,"NEW_COST","variables:",nvars(basering),"cost vector:"); 1017 write(NEW_COST,"NEW_COST","variables:",nvars(basering),"cost vector:"); 1030 1018 for(i=1;i<=nvars(basering);i++) 1031 1019 { 1032 1020 write(NEW_COST,weightvec[i]); 1033 1021 } 1034 1022 1035 1023 // call external program 1036 1024 dummy=system("sh","change_cost "+groebnerfile+" "+newcostfile); 1037 1025 1038 1026 // read toric standard basis from created file 1039 1027 link TORIC_IDEAL=":r "+newcostfile+".GB.pt"; … … 1044 1032 pos=find(toric_ideal,":",pos); 1045 1033 pos++; 1046 1034 1047 1035 while(toric_ideal[pos]==" "  toric_ideal[pos]==newline) 1048 1036 { … … 1085 1073 if(v[i]>0) 1086 1074 { 1087 head=head*var(i)^v[i]; 1075 head=head*var(i)^v[i]; 1088 1076 } 1089 1077 } 1090 1078 J[j]=headtail; 1091 1079 } 1092 1080 1093 1081 // delete all created files 1094 1082 dummy=system("sh","rm f "+groebnerfile); 1095 1083 dummy=system("sh","rm f "+groebnerfile+".GB.pt"); 1096 dummy=system("sh","rm f "+newcostfile); 1084 dummy=system("sh","rm f "+newcostfile); 1097 1085 1098 1086 return(J); 1099 1087 } 1100 1101 1102 1103 1088 example 1104 1089 { 1105 "EXAMPLE"; echo=2;1106 1107 ring r=0,(x,y,z),wp(3,2,1);1108 1109 // call with toric ideal (of the matrix A=(1,1,1))1110 ideal I=xy,xz;1111 ideal J=toric_std(I);1112 J;1113 1114 // call with the same ideal, but badly chosen generators:1115 // not only binomials 1116 I=xy,2xyz;1117 J=toric_std(I);1118 // binomials whose monomials are not relatively prime1119 I=xy,xyyz,yz;1120 J=toric_std(I);1121 J;1122 1123 // call with a nontoric ideal that seems to be toric1124 I=xyz,xyz;1125 J=toric_std(I);1126 J;1127 // comparison with real standard basis and saturation 1128 ideal H=std(I);1129 H;1130 LIB "elim.lib";1131 sat(H,xyz);1090 "EXAMPLE"; echo=2; 1091 1092 ring r=0,(x,y,z),wp(3,2,1); 1093 1094 // call with toric ideal (of the matrix A=(1,1,1)) 1095 ideal I=xy,xz; 1096 ideal J=toric_std(I); 1097 J; 1098 1099 // call with the same ideal, but badly chosen generators: 1100 // not only binomials 1101 I=xy,2xyz; 1102 J=toric_std(I); 1103 // binomials whose monomials are not relatively prime 1104 I=xy,xyyz,yz; 1105 J=toric_std(I); 1106 J; 1107 1108 // call with a nontoric ideal that seems to be toric 1109 I=xyz,xyz; 1110 J=toric_std(I); 1111 J; 1112 // comparison with real standard basis and saturation 1113 ideal H=std(I); 1114 H; 1115 LIB "elim.lib"; 1116 sat(H,xyz); 1132 1117 } 1133 1118 1134 1135 1136 1119 ////////////////////////////////////////////////////////////////////////////// 1137 1138 1139 // toric_std(ideal I); ring term ordering
Note: See TracChangeset
for help on using the changeset viewer.