Changeset 27e935 in git


Ignore:
Timestamp:
May 11, 2000, 11:04:05 AM (23 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
54b3568dedfbe0f7886ed616396038a36933c3e3
Parents:
0658d94b02b071f2a019543baea28fe2c95d9226
Message:
*hannes: formatiing


git-svn-id: file:///usr/local/Singular/svn/trunk@4311 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
IntegerProgramming
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • IntegerProgramming/IP.lib

    r0658d9 r27e935  
    1 // $Id: IP.lib,v 1.1.1.1 2000-05-02 12:58:31 Singular Exp $
     1// $Id: IP.lib,v 1.2 2000-05-11 09:04:05 Singular Exp $
    22//
    33// author : Christine Theis
    44//
    55
    6 //version="$Id: IP.lib,v 1.1.1.1 2000-05-02 12:58:31 Singular Exp $";
     6//version="$Id: IP.lib,v 1.2 2000-05-11 09:04:05 Singular Exp $";
    77
    88///////////////////////////////////////////////////////////////////////////////
    99
    1010info="
    11 LIBRARY: IP.lib         PROCEDURES FOR INTEGER PROGRAMMING USING GROEBNER BASIS METHODS   
     11LIBRARY: IP.lib                PROCEDURES FOR INTEGER PROGRAMMING USING GROEBNER BASIS METHODS
    1212
    1313
    1414Let A an integral (mxn)-matrix, b a vector with m integral
    1515coefficients and c a vector with n nonnegative integral coefficients. The
    16 solution of the IP-problem 
    17          
     16solution of the IP-problem
     17
    1818(*)    minimize{ cx | Ax=b, all components of x are nonnegative integers }
    19              
     19
    2020proceeds in two steps:
    21            
     21
    22221. We compute the toric ideal of A and its Groebner basis with respect
    2323   to a term ordering refining the cost function c (such an ordering
    2424   exists because c is nonnegative).
    25            
     25
    26262. We reduce the right hand vector b or an initial solution of the problem
    2727   modulo this ideal.
    28              
     28
    2929For these purposes, we can use seven different algorithms:
    3030The algorithm of Conti/Traverso (ct) can compute an optimal solution
    3131of the IP-problem directly from the right hand vector b. The same is
    3232true 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. 
     33and b have nonnegative entries, but should be faster in that case.
    3434All other algorithms need initial solutions of the IP-problem. Except
    3535from the Conti-Traverso algorithm with elimination (ect),
    … …  
    4040
    4141
    42 solve_IP(intmat A, intvec bx, intvec c, string alg);   
     42solve_IP(intmat A, intvec bx, intvec c, string alg);
    4343solve_IP(intmat A, list bx, intvec c, string alg);
    44 solve_IP(intmat A, intvec bx, intvec c, string alg, intvec prsv);       
     44solve_IP(intmat A, intvec bx, intvec c, string alg, intvec prsv);
    4545solve_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 
     46
     47        procedures for solving the IP-problem (*)
     48    They return the solution(s) of the given problem(s) or the
    4949    message `not solvable'.
    5050    `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
    5252    hand vector b of the system Ax=b. b should then be an intvec of
    5353    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.
     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.
    5959    If `alg' is chosen to be `blr' or `hs', the algorithm needs a
    6060    vector with positive coefficcients in the row space of A. If
    6161    no row of A contains only positive entries, one must use the
    6262    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                          
     63    argument.
     64        solve_IP may also be called with a list bx of intvecs instead
     65        of a single intvec.
     66
     67
    6868";
    6969
    7070///////////////////////////////////////////////////////////////////////////////
    7171
    72 
    73 
    7472static proc solve_IP_1(intmat A, intvec bx, intvec c, string alg)
    75 {       
     73{
    7674  intvec v;
    7775  // to be returned
    7876
    7977  // check arguments as far as necessary
    80   // other inconsistencies are detected by the external program 
     78  // other inconsistencies are detected by the external program
    8179  if(size(c)!=ncols(A))
    8280  {
    8381    "ERROR: number of matrix columns must equal size of cost vector";
    8482    return(v);
    85   }     
     83  }
    8684
    8785  // create first temporary file with that the external program is
    88   // called 
     86  // called
    8987
    9088  int process=system("pid");
    91   string matrixfile="temp_MATRIX"+string(process);     
     89  string matrixfile="temp_MATRIX"+string(process);
    9290  link MATRIX=":w "+matrixfile;
    9391  open(MATRIX);
    … …  
    107105    }
    108106  }
    109  
     107
    110108  // search for positive row space vector, if required by the
    111   // algorithm   
     109  // algorithm
    112110  int found=0;
    113111  if((alg=="blr") || (alg=="hs"))
    114   { 
     112  {
    115113    for(i=1;i<=nrows(A);i++)
    116114    {
    … …  
    130128    if(found==0)
    131129    {
    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";
    133131      close(MATRIX);
    134132      system("sh","rm -f "+matrixfile);
    135133      return(v);
    136     } 
     134    }
    137135    write(MATRIX,"positive row space vector:");
    138136    for(j=1;j<=ncols(A);j++)
    139     { 
     137    {
    140138      write(MATRIX,A[found,j]);
    141139    }
    … …  
    145143  // create second temporary file for the external program
    146144
    147   string problemfile="temp_PROBLEM"+string(process);   
     145  string problemfile="temp_PROBLEM"+string(process);
    148146  link PROBLEM=":w "+problemfile;
    149147  open(PROBLEM);
    150148
    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:");
    152150  for(i=1;i<=size(bx);i++)
    153151  {
    … …  
    165163  string s;
    166164  if(alg=="ct" || alg=="pct")
    167   { 
     165  {
    168166    pos=find(solution,"NO");
    169167    if(pos!=0)
    … …  
    193191  }
    194192  else
    195   { 
     193  {
    196194    pos=find(solution,"optimal");
    197195    pos=find(solution,":",pos);
    … …  
    212210    }
    213211  }
    214      
     212
    215213  // delete all created files
    216214  dummy=system("sh","rm -f "+matrixfile);
    … …  
    222220}
    223221
    224 
    225 
    226222static proc solve_IP_2(intmat A, list bx, intvec c, string alg)
    227 {       
     223{
    228224  list l;;
    229225  // to be returned
    230226
    231227  // check arguments as far as necessary
    232   // other inconsistencies are detected by the external program 
     228  // other inconsistencies are detected by the external program
    233229  if(size(c)!=ncols(A))
    234230  {
    235231    "ERROR: number of matrix columns must equal size of cost vector";
    236232    return(l);
    237   }     
     233  }
    238234
    239235  int k;
    … …  
    248244
    249245  // create first temporary file with that the external program is
    250   // called 
     246  // called
    251247
    252248  int process=system("pid");
    253   string matrixfile="temp_MATRIX"+string(process);     
     249  string matrixfile="temp_MATRIX"+string(process);
    254250  link MATRIX=":w "+matrixfile;
    255251  open(MATRIX);
    … …  
    269265    }
    270266  }
    271  
     267
    272268  // search for positive row space vector, if required by the
    273   // algorithm   
     269  // algorithm
    274270  int found=0;
    275271  if((alg=="blr") || (alg=="hs"))
    276   { 
     272  {
    277273    for(i=1;i<=nrows(A);i++)
    278274    {
    … …  
    292288    if(found==0)
    293289    {
    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";
    295291      close(MATRIX);
    296292      system("sh","rm -f "+matrixfile);
    297293      return(l);
    298     } 
     294    }
    299295    write(MATRIX,"positive row space vector:");
    300296    for(j=1;j<=ncols(A);j++)
    301     { 
     297    {
    302298      write(MATRIX,A[found,j]);
    303299    }
    … …  
    307303  // create second temporary file for the external program
    308304
    309   string problemfile="temp_PROBLEM"+string(process);   
     305  string problemfile="temp_PROBLEM"+string(process);
    310306  link PROBLEM=":w "+problemfile;
    311307  open(PROBLEM);
    312308
    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:");
    314310  for(k=1;k<=size(bx);k++)
    315311  {
    … …  
    331327  string s;
    332328  if(alg=="ct" || alg=="pct")
    333   { 
     329  {
    334330    pos=1;
    335331    for(k=1;k<=size(bx);k++)
    336     {     
     332    {
    337333      pos1=find(solution,"NO",pos);
    338334      pos2=find(solution,"YES",pos);
    … …  
    359355            s=s+solution[pos];
    360356            pos++;
    361           } 
     357          }
    362358          execute("v[j]="+s+";");
    363359        }
    … …  
    367363  }
    368364  else
    369   { 
     365  {
    370366    pos=1;
    371367    for(k=1;k<=size(bx);k++)
    … …  
    385381          s=s+solution[pos];
    386382          pos++;
    387         } 
     383        }
    388384        execute("v[j]="+s+";");
    389385      }
    … …  
    391387    }
    392388  }
    393      
     389
    394390  // delete all created files
    395391  dummy=system("sh","rm -f "+matrixfile);
    … …  
    401397}
    402398
    403 
    404 
    405399static proc solve_IP_3(intmat A, intvec bx, intvec c, string alg, intvec prsv)
    406 {       
     400{
    407401  intvec v;
    408402  // to be returned
    409403
    410404  // check arguments as far as necessary
    411   // other inconsistencies are detected by the external program 
     405  // other inconsistencies are detected by the external program
    412406  if(size(c)!=ncols(A))
    413407  {
    414408    "ERROR: number of matrix columns must equal size of cost vector";
    415409    return(v);
    416   }     
     410  }
    417411
    418412  if(size(prsv)!=ncols(A))
    … …  
    420414    "ERROR: number of matrix columns must equal size of positive row space vector";
    421415    return(v);
    422   }     
     416  }
    423417
    424418  // create first temporary file with that the external program is
    425   // called 
     419  // called
    426420
    427421  int process=system("pid");
    428   string matrixfile="temp_MATRIX"+string(process);     
     422  string matrixfile="temp_MATRIX"+string(process);
    429423  link MATRIX=":w "+matrixfile;
    430424  open(MATRIX);
    … …  
    444438    }
    445439  }
    446  
     440
    447441  // enter positive row space vector, if required by the algorithm
    448442  if((alg=="blr") || (alg=="hs"))
    449   { 
     443  {
    450444    write(MATRIX,"positive row space vector:");
    451445    for(j=1;j<=ncols(A);j++)
    452     { 
     446    {
    453447      write(MATRIX,prsv[j]);
    454448    }
    … …  
    458452  // create second temporary file for the external program
    459453
    460   string problemfile="temp_PROBLEM"+string(process);   
     454  string problemfile="temp_PROBLEM"+string(process);
    461455  link PROBLEM=":w "+problemfile;
    462456  open(PROBLEM);
    463457
    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:");
    465459  for(i=1;i<=size(bx);i++)
    466460  {
    … …  
    478472  string s;
    479473  if(alg=="ct" || alg=="pct")
    480   { 
     474  {
    481475    pos=find(solution,"NO");
    482476    if(pos!=0)
    … …  
    500494          s=s+solution[pos];
    501495          pos++;
    502         } 
     496        }
    503497        execute("v[j]="+s+";");
    504498      }
    … …  
    506500  }
    507501  else
    508   { 
     502  {
    509503    pos=find(solution,"optimal");
    510504    pos=find(solution,":",pos);
    … …  
    521515        s=s+solution[pos];
    522516        pos++;
    523       } 
     517      }
    524518      execute("v[j]="+s+";");
    525519    }
    526520  }
    527      
     521
    528522  // delete all created files
    529523  dummy=system("sh","rm -f "+matrixfile);
    … …  
    535529}
    536530
    537 
    538 
    539531static proc solve_IP_4(intmat A, list bx, intvec c, string alg, intvec prsv)
    540 {       
     532{
    541533  list l;
    542534  // to be returned
    543535
    544536  // check arguments as far as necessary
    545   // other inconsistencies are detected by the external program 
     537  // other inconsistencies are detected by the external program
    546538  if(size(c)!=ncols(A))
    547539  {
    548540    "ERROR: number of matrix columns must equal size of cost vector";
    549541    return(l);
    550   }     
     542  }
    551543
    552544  if(size(prsv)!=ncols(A))
    … …  
    554546    "ERROR: number of matrix columns must equal size of positive row space vector";
    555547    return(v);
    556   }     
     548  }
    557549
    558550  int k;
    … …  
    567559
    568560  // create first temporary file with that the external program is
    569   // called 
     561  // called
    570562
    571563  int process=system("pid");
    572   string matrixfile="temp_MATRIX"+string(process);     
     564  string matrixfile="temp_MATRIX"+string(process);
    573565  link MATRIX=":w "+matrixfile;
    574566  open(MATRIX);
    … …  
    588580    }
    589581  }
    590  
     582
    591583  // enter positive row space vector if required by the algorithm
    592584  if((alg=="blr") || (alg=="hs"))
    593   { 
     585  {
    594586    write(MATRIX,"positive row space vector:");
    595587    for(j=1;j<=ncols(A);j++)
    596     { 
     588    {
    597589      write(MATRIX,prsv[j]);
    598590    }
    … …  
    602594  // create second temporary file for the external program
    603595
    604   string problemfile="temp_PROBLEM"+string(process);   
     596  string problemfile="temp_PROBLEM"+string(process);
    605597  link PROBLEM=":w "+problemfile;
    606598  open(PROBLEM);
    607599
    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:");
    609601  for(k=1;k<=size(bx);k++)
    610602  {
    … …  
    626618  string s;
    627619  if(alg=="ct" || alg=="pct")
    628   { 
     620  {
    629621    pos=1;
    630622    for(k=1;k<=size(bx);k++)
    631     {     
     623    {
    632624      pos1=find(solution,"NO",pos);
    633625      pos2=find(solution,"YES",pos);
    … …  
    654646            s=s+solution[pos];
    655647            pos++;
    656           } 
     648          }
    657649          execute("v[j]="+s+";");
    658650        }
    … …  
    662654  }
    663655  else
    664   { 
     656  {
    665657    pos=1;
    666658    for(k=1;k<=size(bx);k++)
    … …  
    686678    }
    687679  }
    688      
     680
    689681  // delete all created files
    690682  dummy=system("sh","rm -f "+matrixfile);
    … …  
    696688}
    697689
    698 
    699 
    700 
    701 
    702 
    703690proc solve_IP
    704 "USAGE:   
     691"USAGE:
    705692solve_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 
     693solve_IP(A,bx,c,alg);       A intmat, bx list of intvec, c intvec, alg string
    707694solve_IP(A,bx,c,alg,prsv);  A intmat, bx intvec, c intvec, alg string,
    708                             prsv intvec 
     695                            prsv intvec
    709696solve_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
     698RETURN:  solution of the associated integer programming problem as explained
    712699         in IP.lib
    713          return type = type of bx 
     700         return type = type of bx
    714701EXAMPLE: example solve_IP;  shows an example"
    715702{
    … …  
    737724  }
    738725}
    739 
    740 
    741 
    742726example
    743727{
    744 "EXAMPLE"; echo=2;
    745 
    746 // call with single right hand vector
    747 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 vectors
    757 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 vector
    764 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 positive
    770 // row space vector
    771 solution_vector=solve_IP(A,b1,c,"hs");
    772 
    773 // call with single initial solution vector and positive row space vector
    774 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 vector
    780 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;
    785769}
    786 
    787 
    788 

  • IntegerProgramming/toric.lib

    r0658d9 r27e935  
    1 // $Id: toric.lib,v 1.1.1.1 2000-05-02 12:58:31 Singular Exp $
     1// $Id: toric.lib,v 1.2 2000-05-11 09:03:41 Singular Exp $
    22//
    33// author : Christine Theis
    44//
    5 
    6 //version="$Id: toric.lib,v 1.1.1.1 2000-05-02 12:58:31 Singular Exp $";
     5//version="$Id: toric.lib,v 1.2 2000-05-11 09:03:41 Singular Exp $";
    76
    87///////////////////////////////////////////////////////////////////////////////
    98
    109info="
    11 LIBRARY: toric.lib              PROCEDURES FOR COMPUTING TORIC IDEALS   
    12 
    13 
    14 Let A an integral (mxn)-matrix. The toric ideal I_A of A is defined 
     10LIBRARY: toric.lib                PROCEDURES FOR COMPUTING TORIC IDEALS
     11
     12
     13Let A an integral (mxn)-matrix. The toric ideal I_A of A is defined
    1514as the ideal
    1615
    … …  
    1817
    1918in 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.   
     19Toric ideals play an important role in polyhedral geometry and may also be
     20used for integer programming. They are generated by binomials with
     21relatively prime monomials. Buchberger's algorithm can be specialized to
     22these structures in a way that considerably speeds up computation.
    2423
    2524
    2625toric_ideal(intmat A, string alg);
    2726toric_ideal(intmat A, string alg, intvec prsv);
    28  
     27
    2928    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
    3130    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.
    3332    Not all term orderings are supported: The usual global term orderings
    3433    may be used, but no block orderings combining them.
    … …  
    4241
    4342    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.
    4645    If `alg' is chosen to be `blr' or `hs', the algorithm needs a
    4746    vector with positive coefficcients in the row space of A. If
    4847    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
    5251
    5352toric_std(ideal I);
    5453
    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
    5756    of binomials of the form x^u-x^v (although there are other generating
    5857    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
    6059    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
    6261    `between' I and its saturation with respect to all variables.
    63                          
     62
    6463";
    6564
    6665///////////////////////////////////////////////////////////////////////////////
    6766
    68 
    69 
    7067static proc toric_ideal_1(intmat A, string alg)
    71 {       
     68{
    7269  ideal I;
    7370  // to be returned
    … …  
    7875    "ERROR: number of matrix columns must be greater or equal number of ring variables";
    7976    return(I);
    80   }     
     77  }
    8178
    8279  // check suitability of actual term ordering
    … …  
    167164      "Warning: block orderings are not supported; Dp used for computation";
    168165    }
    169   } 
     166  }
    170167
    171168  int pos;
    … …  
    263260    return(I);
    264261  }
    265  
     262
    266263  // check algorithm
    267264  if(alg=="ct" || alg=="pct")
    … …  
    274271  }
    275272
    276   // create temporary file with that the external program is called 
     273  // create temporary file with that the external program is called
    277274
    278275  int dummy;
    279276  int process=system("pid");
    280   string matrixfile="temp_MATRIX"+string(process);     
     277  string matrixfile="temp_MATRIX"+string(process);
    281278  link MATRIX=":w "+matrixfile;
    282279  open(MATRIX);
    … …  
    295292    }
    296293  }
    297  
     294
    298295  // search for positive row space vector, if required by the
    299   // algorithm   
     296  // algorithm
    300297  int found=0;
    301298  if((alg=="blr") || (alg=="hs"))
    302   { 
     299  {
    303300    for(i=1;i<=nrows(A);i++)
    304301    {
    … …  
    318315    if(found==0)
    319316    {
    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";
    321318      close(MATRIX);
    322319      dummy=system("sh","rm -f "+matrixfile);
    323320      return(I);
    324     } 
     321    }
    325322    write(MATRIX,"positive row space vector:");
    326323    for(j=1;j<=ncols(A);j++)
    327     { 
     324    {
    328325      write(MATRIX,A[found,j]);
    329326    }
    … …  
    343340  pos=find(toric_ideal,":",pos);
    344341  pos++;
    345  
     342
    346343  while(toric_ideal[pos]==" " || toric_ideal[pos]==newline)
    347344  {
    … …  
    355352  }
    356353  execute("generators="+number_string+";");
    357  
     354
    358355  intvec v;
    359356  poly head;
    … …  
    388385      if(v[j]>0)
    389386      {
    390         head=head*var(j)^v[j];   
     387        head=head*var(j)^v[j];
    391388      }
    392389    }
    393390    I[i]=head-tail;
    394391  }
    395      
     392
    396393  // delete all created files
    397394  dummy=system("sh","rm -f "+matrixfile);
    … …  
    401398}
    402399
    403 
    404 
    405400static proc toric_ideal_2(intmat A, string alg, intvec prsv)
    406 {       
     401{
    407402  ideal I;
    408403  // to be returned
    … …  
    413408    "ERROR: number of matrix columns must equal size of positive row space vector";
    414409    return(I);
    415   }     
     410  }
    416411
    417412  // check suitability of actual basering
    … …  
    420415    "ERROR: number of matrix columns must be greater or equal number of ring variables";
    421416    return(I);
    422   }     
     417  }
    423418
    424419  // check suitability of actual term ordering
    … …  
    509504      "Warning: block orderings are not supported; Dp used for computation";
    510505    }
    511   } 
     506  }
    512507
    513508  int pos;
    … …  
    605600    return(I);
    606601  }
    607  
     602
    608603  // check algorithm
    609604  if(alg=="ct" || alg=="pct")
    … …  
    616611  }
    617612
    618   // create temporary file with that the external program is called 
     613  // create temporary file with that the external program is called
    619614
    620615  int dummy;
    621616  int process=system("pid");
    622   string matrixfile="temp_MATRIX"+string(process);     
     617  string matrixfile="temp_MATRIX"+string(process);
    623618  link MATRIX=":w "+matrixfile;
    624619  open(MATRIX);
    … …  
    640635  // enter positive row space vector, if required by the algorithm
    641636  if((alg=="blr") || (alg=="hs"))
    642   { 
     637  {
    643638    write(MATRIX,"positive row space vector:");
    644639    for(j=1;j<=ncols(A);j++)
    645     { 
     640    {
    646641      write(MATRIX,prsv[j]);
    647642    }
    … …  
    660655  pos=find(toric_ideal,":",pos);
    661656  pos++;
    662  
     657
    663658  while(toric_ideal[pos]==" " || toric_ideal[pos]==newline)
    664659  {
    … …  
    672667  }
    673668  execute("generators="+number_string+";");
    674  
     669
    675670  intvec v;
    676671  poly head;
    … …  
    705700      if(v[j]>0)
    706701      {
    707         head=head*var(j)^v[j];   
     702        head=head*var(j)^v[j];
    708703      }
    709704    }
    710705    I[i]=head-tail;
    711706  }
    712      
     707
    713708  // delete all created files
    714709  dummy=system("sh","rm -f "+matrixfile);
    … …  
    718713}
    719714
    720 
    721 
    722715proc toric_ideal
    723 "USAGE:   
     716"USAGE:
    724717toric_ideal(A,alg);            A intmat, alg string
    725 toric_ideal(A,alg,prsv);       A intmat, alg string, prsv intvec 
     718toric_ideal(A,alg,prsv);       A intmat, alg string, prsv intvec
    726719RETURN:  toric ideal of A as explained in toric.lib
    727720         return type = ideal
    … …  
    737730  }
    738731}
    739 
    740 
    741 
    742732example
    743733{
    744 "EXAMPLE"; echo=2;
    745 
    746 ring r=0,(x,y,z),dp;
    747 
    748 // call with two arguments
    749 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 arguments
    759 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
    763753}
    764754
    765 
    766 
    767755proc toric_std(ideal I)
    768 "USAGE:   toric_std(I);        I ideal       
     756"USAGE:   toric_std(I);        I ideal
    769757RETURN:  standard basis of I as explained in toric.lib
    770758         return type = ideal
    … …  
    861849      "Warning: block orderings are not supported; Dp used for computation";
    862850    }
    863   } 
     851  }
    864852
    865853  int pos;
    … …  
    887875    }
    888876  }
    889  
     877
    890878  if(external_ord=="" && find(singular_ord,"wp")==3)
    891879  {
    … …  
    958946    return(I);
    959947  }
    960  
     948
    961949  // create first temporary file with which the external program is called
    962  
     950
    963951  int dummy;
    964952  int process=system("pid");
    965   string groebnerfile="temp_GROEBNER"+string(process); 
     953  string groebnerfile="temp_GROEBNER"+string(process);
    966954  link GROEBNER=":w "+groebnerfile;
    967955  open(GROEBNER);
    … …  
    969957  write(GROEBNER,"GROEBNER","computed with algorithm:","pt","term ordering:","elimination block",0,"weighted block",nvars(basering),external_ord);
    970958  // algorithm is totally unimportant, only required by the external program
    971  
     959
    972960  for(i=1;i<=nvars(basering);i++)
    973961  {
    974962    write(GROEBNER,weightvec[i]);
    975963  }
    976  
     964
    977965  write(GROEBNER,"size:",size(I),"Groebner basis:");
    978966  poly head;
    … …  
    983971  for(j=1;j<=size(I);j++)
    984972  {
    985     // test suitability of generator j 
     973    // test suitability of generator j
    986974    rest=I[j];
    987975    head=lead(rest);
    … …  
    989977    tail=lead(rest);
    990978    rest=rest-tail;
    991    
     979
    992980    if(head==0 && tail==0 && rest!=0)
    993981    {
    … …  
    997985      return(J);
    998986    }
    999    
     987
    1000988    if(leadcoef(tail)!=-leadcoef(head))
    1001989      // generator no difference of monomials (or a constant multiple)
    … …  
    1006994      return(J);
    1007995    }
    1008        
     996
    1009997    if(gcd(head,tail)!=1)
    1010998    {
    1011999      "Warning: monomials of generator "+string(j)+" of input ideal are not relatively prime";
    10121000    }
    1013    
     1001
    10141002    // write vector representation of generator j into the file
    10151003    v=leadexp(head)-leadexp(tail);
    … …  
    10201008  }
    10211009  close(GROEBNER);
    1022  
     1010
    10231011  // create second temporary file
    10241012
    1025   string newcostfile="temp_NEW_COST"+string(process);   
     1013  string newcostfile="temp_NEW_COST"+string(process);
    10261014  link NEW_COST=":w "+newcostfile;
    10271015  open(NEW_COST);
    10281016
    1029   write(NEW_COST,"NEW_COST","variables:",nvars(basering),"cost vector:"); 
     1017  write(NEW_COST,"NEW_COST","variables:",nvars(basering),"cost vector:");
    10301018  for(i=1;i<=nvars(basering);i++)
    10311019  {
    10321020    write(NEW_COST,weightvec[i]);
    10331021  }
    1034  
     1022
    10351023  // call external program
    10361024  dummy=system("sh","change_cost "+groebnerfile+" "+newcostfile);
    1037  
     1025
    10381026  // read toric standard basis from created file
    10391027  link TORIC_IDEAL=":r "+newcostfile+".GB.pt";
    … …  
    10441032  pos=find(toric_ideal,":",pos);
    10451033  pos++;
    1046  
     1034
    10471035  while(toric_ideal[pos]==" " || toric_ideal[pos]==newline)
    10481036  {
    … …  
    10851073      if(v[i]>0)
    10861074      {
    1087         head=head*var(i)^v[i];   
     1075        head=head*var(i)^v[i];
    10881076      }
    10891077    }
    10901078    J[j]=head-tail;
    10911079  }
    1092      
     1080
    10931081  // delete all created files
    10941082  dummy=system("sh","rm -f "+groebnerfile);
    10951083  dummy=system("sh","rm -f "+groebnerfile+".GB.pt");
    1096   dummy=system("sh","rm -f "+newcostfile); 
     1084  dummy=system("sh","rm -f "+newcostfile);
    10971085
    10981086  return(J);
    10991087}
    1100 
    1101 
    1102 
    11031088example
    11041089{
    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=x-y,x-z;
    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=x-y,2x-y-z;
    1117 J=toric_std(I);
    1118 // binomials whose monomials are not relatively prime
    1119 I=x-y,xy-yz,y-z;
    1120 J=toric_std(I);
    1121 J;
    1122  
    1123 // call with a non-toric ideal that seems to be toric
    1124 I=x-yz,xy-z;
    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=x-y,x-z;
     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=x-y,2x-y-z;
     1102  J=toric_std(I);
     1103  // binomials whose monomials are not relatively prime
     1104  I=x-y,xy-yz,y-z;
     1105  J=toric_std(I);
     1106  J;
     1107  
     1108  // call with a non-toric ideal that seems to be toric
     1109  I=x-yz,xy-z;
     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);
    11321117}
    11331118
    1134 
    1135 
    11361119//////////////////////////////////////////////////////////////////////////////
    1137 
    1138 
    1139 // toric_std(ideal I); ring term ordering
Note: See TracChangeset for help on using the changeset viewer.