Changeset 27e935 in git for IntegerProgramming/toric.lib


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


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

Legend:

Unmodified
Added
Removed
  • 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.