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

// $Id: IP.lib,v 1.2 2000-05-11 09:04:05 Singular Exp $
//
// author : Christine Theis
//

//version="$Id: IP.lib,v 1.2 2000-05-11 09:04:05 Singular Exp $";

///////////////////////////////////////////////////////////////////////////////

info="
LIBRARY: IP.lib                PROCEDURES FOR INTEGER PROGRAMMING USING GROEBNER BASIS METHODS


Let A an integral (mxn)-matrix, b a vector with m integral
coefficients and c a vector with n nonnegative integral coefficients. The
solution of the IP-problem

(*)    minimize{ cx | Ax=b, all components of x are nonnegative integers }

proceeds in two steps:

1. We compute the toric ideal of A and its Groebner basis with respect
   to a term ordering refining the cost function c (such an ordering
   exists because c is nonnegative).

2. We reduce the right hand vector b or an initial solution of the problem
   modulo this ideal.

For these purposes, we can use seven different algorithms:
The algorithm of Conti/Traverso (ct) can compute an optimal solution
of the IP-problem directly from the right hand vector b. The same is
true for its `positive' variant (pct) which can only be applied if A
and b have nonnegative entries, but should be faster in that case.
All other algorithms need initial solutions of the IP-problem. Except
from the Conti-Traverso algorithm with elimination (ect),
they should be more efficient than the algorithm mentionned before.
These are the algorithms of Pottier (pt), Bigatti/La Scala/Robbiano
(blr), Hosten/Sturmfels (hs) and Di Biase/Urbanke (du). The last two
seem to be the fastest in the actual implementation.


solve_IP(intmat A, intvec bx, intvec c, string alg);
solve_IP(intmat A, list bx, intvec c, string alg);
solve_IP(intmat A, intvec bx, intvec c, string alg, intvec prsv);
solve_IP(intmat A, list bx, intvec c, string alg, intvec prsv);

        procedures for solving the IP-problem (*)
    They return the solution(s) of the given problem(s) or the
    message `not solvable'.
    `alg' may be one of the algorithm abbreviations as above.
    If `alg' is chosen to be `ct' or `pct', bx is read as the right
    hand vector b of the system Ax=b. b should then be an intvec of
    size m where m is the number of rows of A. Furthermore, bx and
    A should be nonnegative if `pct' is used.
        If `alg' is chosen to be `ect',`pt',`blr',`hs' or `du', bx is
        read as an initial solution x of the system Ax=b. bx should
        then be a nonnegative intvec of size n where n is the number
        of columns of A.
    If `alg' is chosen to be `blr' or `hs', the algorithm needs a
    vector with positive coefficcients in the row space of A. If
    no row of A contains only positive entries, one must use the
    versions of solve_IP which take such a vector prsv as
    argument.
        solve_IP may also be called with a list bx of intvecs instead
        of a single intvec.


";

///////////////////////////////////////////////////////////////////////////////

static proc solve_IP_1(intmat A, intvec bx, intvec c, string alg)
{
  intvec v;
  // to be returned

  // check arguments as far as necessary
  // other inconsistencies are detected by the external program
  if(size(c)!=ncols(A))
  {
    "ERROR: number of matrix columns must equal size of cost vector";
    return(v);
  }

  // create first temporary file with that the external program is
  // called

  int process=system("pid");
  string matrixfile="temp_MATRIX"+string(process);
  link MATRIX=":w "+matrixfile;
  open(MATRIX);

  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
  int i,j;
  for(j=1;j<=ncols(A);j++)
  {
    write(MATRIX,c[j]);
  }
  write(MATRIX,"rows:",nrows(A),"matrix:");
  for(i=1;i<=nrows(A);i++)
  {
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,A[i,j]);
    }
  }

  // search for positive row space vector, if required by the
  // algorithm
  int found=0;
  if((alg=="blr") || (alg=="hs"))
  {
    for(i=1;i<=nrows(A);i++)
    {
      found=i;
      for(j=1;j<=ncols(A);j++)
      {
        if(A[i,j]<=0)
        {
          found=0;
        }
      }
      if(found>0)
      {
        break;
      }
    }
    if(found==0)
    {
      "ERROR: algorithm needs positive vector in the row space of the matrix";
      close(MATRIX);
      system("sh","rm -f "+matrixfile);
      return(v);
    }
    write(MATRIX,"positive row space vector:");
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,A[found,j]);
    }
  }
  close(MATRIX);

  // create second temporary file for the external program

  string problemfile="temp_PROBLEM"+string(process);
  link PROBLEM=":w "+problemfile;
  open(PROBLEM);

  write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:");
  for(i=1;i<=size(bx);i++)
  {
    write(PROBLEM,bx[i]);
  }
  close(PROBLEM);

  // call external program
  int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);

  // read solution from created file
  link SOLUTION=":r "+matrixfile+".sol."+alg;
  string solution=read(SOLUTION);
  int pos;
  string s;
  if(alg=="ct" || alg=="pct")
  {
    pos=find(solution,"NO");
    if(pos!=0)
    {
      "not solvable";
    }
    else
    {
      pos=find(solution,"YES");
      pos=find(solution,":",pos);
      pos++;
      for(j=1;j<=ncols(A);j++)
      {
        while(solution[pos]==" " || solution[pos]==newline)
        {
          pos++;
        }
        s="";
        while(solution[pos]!=" " && solution[pos]!=newline)
        {
          s=s+solution[pos];
          pos++;
        }
        execute("v[j]="+s+";");
      }
    }
  }
  else
  {
    pos=find(solution,"optimal");
    pos=find(solution,":",pos);
    pos++;
    for(j=1;j<=ncols(A);j++)
    {
      while(solution[pos]==" " || solution[pos]==newline)
      {
        pos++;
      }
      s="";
      while(solution[pos]!=" " && solution[pos]!=newline)
      {
        s=s+solution[pos];
        pos++;
      }
      execute("v[j]="+s+";");
    }
  }

  // delete all created files
  dummy=system("sh","rm -f "+matrixfile);
  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
  dummy=system("sh","rm -f "+problemfile);
  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);

  return(v);
}

static proc solve_IP_2(intmat A, list bx, intvec c, string alg)
{
  list l;;
  // to be returned

  // check arguments as far as necessary
  // other inconsistencies are detected by the external program
  if(size(c)!=ncols(A))
  {
    "ERROR: number of matrix columns must equal size of cost vector";
    return(l);
  }

  int k;
  for(k=2;k<=size(bx);k++)
  {
    if(size(bx[k])!=size(bx[1]))
    {
      "ERROR: size of all right hand vectors must be equal";
      return(l);
    }
  }

  // create first temporary file with that the external program is
  // called

  int process=system("pid");
  string matrixfile="temp_MATRIX"+string(process);
  link MATRIX=":w "+matrixfile;
  open(MATRIX);

  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
  int i,j;
  for(j=1;j<=ncols(A);j++)
  {
    write(MATRIX,c[j]);
  }
  write(MATRIX,"rows:",nrows(A),"matrix:");
  for(i=1;i<=nrows(A);i++)
  {
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,A[i,j]);
    }
  }

  // search for positive row space vector, if required by the
  // algorithm
  int found=0;
  if((alg=="blr") || (alg=="hs"))
  {
    for(i=1;i<=nrows(A);i++)
    {
      found=i;
      for(j=1;j<=ncols(A);j++)
      {
        if(A[i,j]<=0)
        {
          found=0;
        }
      }
      if(found>0)
      {
        break;
      }
    }
    if(found==0)
    {
      "ERROR: algorithm needs positive vector in the row space of the matrix";
      close(MATRIX);
      system("sh","rm -f "+matrixfile);
      return(l);
    }
    write(MATRIX,"positive row space vector:");
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,A[found,j]);
    }
  }
  close(MATRIX);

  // create second temporary file for the external program

  string problemfile="temp_PROBLEM"+string(process);
  link PROBLEM=":w "+problemfile;
  open(PROBLEM);

  write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:");
  for(k=1;k<=size(bx);k++)
  {
    for(i=1;i<=size(bx[1]);i++)
    {
      write(PROBLEM,bx[k][i]);
    }
  }
  close(PROBLEM);

  // call external program
  int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);

  // read solution from created file
  link SOLUTION=":r "+matrixfile+".sol."+alg;
  string solution=read(SOLUTION);
  intvec v;
  int pos,pos1,pos2;
  string s;
  if(alg=="ct" || alg=="pct")
  {
    pos=1;
    for(k=1;k<=size(bx);k++)
    {
      pos1=find(solution,"NO",pos);
      pos2=find(solution,"YES",pos);
      if(pos1!=0 && (pos1<pos2 || pos2==0))
        // first problem not solvable
      {
        pos=find(solution,":",pos1);
        l=insert(l,"not solvable",size(l));
      }
      else
        // first problem solvable
      {
        pos=find(solution,":",pos2);
        pos++;
        for(j=1;j<=ncols(A);j++)
        {
          while(solution[pos]==" " || solution[pos]==newline)
          {
            pos++;
          }
          s="";
          while(solution[pos]!=" " && solution[pos]!=newline)
          {
            s=s+solution[pos];
            pos++;
          }
          execute("v[j]="+s+";");
        }
        l=insert(l,v,size(l));
      }
    }
  }
  else
  {
    pos=1;
    for(k=1;k<=size(bx);k++)
    {
      pos=find(solution,"optimal",pos);
      pos=find(solution,":",pos);
      pos++;
      for(j=1;j<=ncols(A);j++)
      {
        while(solution[pos]==" " || solution[pos]==newline)
        {
          pos++;
        }
        s="";
        while(solution[pos]!=" " && solution[pos]!=newline)
        {
          s=s+solution[pos];
          pos++;
        }
        execute("v[j]="+s+";");
      }
      l=insert(l,v,size(l));
    }
  }

  // delete all created files
  dummy=system("sh","rm -f "+matrixfile);
  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
  dummy=system("sh","rm -f "+problemfile);
  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);

  return(l);
}

static proc solve_IP_3(intmat A, intvec bx, intvec c, string alg, intvec prsv)
{
  intvec v;
  // to be returned

  // check arguments as far as necessary
  // other inconsistencies are detected by the external program
  if(size(c)!=ncols(A))
  {
    "ERROR: number of matrix columns must equal size of cost vector";
    return(v);
  }

  if(size(prsv)!=ncols(A))
  {
    "ERROR: number of matrix columns must equal size of positive row space vector";
    return(v);
  }

  // create first temporary file with that the external program is
  // called

  int process=system("pid");
  string matrixfile="temp_MATRIX"+string(process);
  link MATRIX=":w "+matrixfile;
  open(MATRIX);

  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
  int i,j;
  for(j=1;j<=ncols(A);j++)
  {
    write(MATRIX,c[j]);
  }
  write(MATRIX,"rows:",nrows(A),"matrix:");
  for(i=1;i<=nrows(A);i++)
  {
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,A[i,j]);
    }
  }

  // enter positive row space vector, if required by the algorithm
  if((alg=="blr") || (alg=="hs"))
  {
    write(MATRIX,"positive row space vector:");
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,prsv[j]);
    }
  }
  close(MATRIX);

  // create second temporary file for the external program

  string problemfile="temp_PROBLEM"+string(process);
  link PROBLEM=":w "+problemfile;
  open(PROBLEM);

  write(PROBLEM,"PROBLEM","vector size:",size(bx),"number of instances:",1,"right hand or initial solution vectors:");
  for(i=1;i<=size(bx);i++)
  {
    write(PROBLEM,bx[i]);
  }
  close(PROBLEM);

  // call external program
  int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);

  // read solution from created file
  link SOLUTION=":r "+matrixfile+".sol."+alg;
  string solution=read(SOLUTION);
  int pos;
  string s;
  if(alg=="ct" || alg=="pct")
  {
    pos=find(solution,"NO");
    if(pos!=0)
    {
      "not solvable";
    }
    else
    {
      pos=find(solution,"YES");
      pos=find(solution,":",pos);
      pos++;
      for(j=1;j<=ncols(A);j++)
      {
        while(solution[pos]==" " || solution[pos]==newline)
        {
          pos++;
        }
        s="";
        while(solution[pos]!=" " && solution[pos]!=newline)
        {
          s=s+solution[pos];
          pos++;
        }
        execute("v[j]="+s+";");
      }
    }
  }
  else
  {
    pos=find(solution,"optimal");
    pos=find(solution,":",pos);
    pos++;
    for(j=1;j<=ncols(A);j++)
    {
      while(solution[pos]==" " || solution[pos]==newline)
      {
        pos++;
      }
      s="";
      while(solution[pos]!=" " && solution[pos]!=newline)
      {
        s=s+solution[pos];
        pos++;
      }
      execute("v[j]="+s+";");
    }
  }

  // delete all created files
  dummy=system("sh","rm -f "+matrixfile);
  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
  dummy=system("sh","rm -f "+problemfile);
  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);

  return(v);
}

static proc solve_IP_4(intmat A, list bx, intvec c, string alg, intvec prsv)
{
  list l;
  // to be returned

  // check arguments as far as necessary
  // other inconsistencies are detected by the external program
  if(size(c)!=ncols(A))
  {
    "ERROR: number of matrix columns must equal size of cost vector";
    return(l);
  }

  if(size(prsv)!=ncols(A))
  {
    "ERROR: number of matrix columns must equal size of positive row space vector";
    return(v);
  }

  int k;
  for(k=2;k<=size(bx);k++)
  {
    if(size(bx[k])!=size(bx[1]))
    {
      "ERROR: size of all right hand vectors must be equal";
      return(l);
    }
  }

  // create first temporary file with that the external program is
  // called

  int process=system("pid");
  string matrixfile="temp_MATRIX"+string(process);
  link MATRIX=":w "+matrixfile;
  open(MATRIX);

  write(MATRIX,"MATRIX","columns:",ncols(A),"cost vector:");
  int i,j;
  for(j=1;j<=ncols(A);j++)
  {
    write(MATRIX,c[j]);
  }
  write(MATRIX,"rows:",nrows(A),"matrix:");
  for(i=1;i<=nrows(A);i++)
  {
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,A[i,j]);
    }
  }

  // enter positive row space vector if required by the algorithm
  if((alg=="blr") || (alg=="hs"))
  {
    write(MATRIX,"positive row space vector:");
    for(j=1;j<=ncols(A);j++)
    {
      write(MATRIX,prsv[j]);
    }
  }
  close(MATRIX);

  // create second temporary file for the external program

  string problemfile="temp_PROBLEM"+string(process);
  link PROBLEM=":w "+problemfile;
  open(PROBLEM);

  write(PROBLEM,"PROBLEM","vector size:",size(bx[1]),"number of instances:",size(bx),"right hand or initial solution vectors:");
  for(k=1;k<=size(bx);k++)
  {
    for(i=1;i<=size(bx[1]);i++)
    {
      write(PROBLEM,bx[k][i]);
    }
  }
  close(PROBLEM);

  // call external program
  int dummy=system("sh","solve_IP -alg "+alg+" "+matrixfile+" "+problemfile);

  // read solution from created file
  link SOLUTION=":r "+matrixfile+".sol."+alg;
  string solution=read(SOLUTION);
  intvec v;
  int pos,pos1,pos2;
  string s;
  if(alg=="ct" || alg=="pct")
  {
    pos=1;
    for(k=1;k<=size(bx);k++)
    {
      pos1=find(solution,"NO",pos);
      pos2=find(solution,"YES",pos);
      if(pos1!=0 && (pos1<pos2 || pos2==0))
        // first problem not solvable
      {
        pos=find(solution,":",pos1);
        l=insert(l,"not solvable",size(l));
      }
      else
        // first problem solvable
      {
        pos=find(solution,":",pos2);
        pos++;
        for(j=1;j<=ncols(A);j++)
        {
          while(solution[pos]==" " || solution[pos]==newline)
          {
            pos++;
          }
          s="";
          while(solution[pos]!=" " && solution[pos]!=newline)
          {
            s=s+solution[pos];
            pos++;
          }
          execute("v[j]="+s+";");
        }
        l=insert(l,v,size(l));
      }
    }
  }
  else
  {
    pos=1;
    for(k=1;k<=size(bx);k++)
    {
      pos=find(solution,"optimal",pos);
      pos=find(solution,":",pos);
      pos++;
      for(j=1;j<=ncols(A);j++)
      {
        while(solution[pos]==" " || solution[pos]==newline)
        {
          pos++;
        }
        s="";
        while(solution[pos]!=" " && solution[pos]!=newline)
        {
          s=s+solution[pos];
          pos++;
        }
        execute("v[j]="+s+";");
      }
      l=insert(l,v,size(l));
    }
  }

  // delete all created files
  dummy=system("sh","rm -f "+matrixfile);
  dummy=system("sh","rm -f "+matrixfile+".GB."+alg);
  dummy=system("sh","rm -f "+problemfile);
  dummy=system("sh","rm -f "+matrixfile+".sol."+alg);

  return(l);
}

proc solve_IP
"USAGE:
solve_IP(A,bx,c,alg);       A intmat, bx intvec, c intvec, alg string
solve_IP(A,bx,c,alg);       A intmat, bx list of intvec, c intvec, alg string
solve_IP(A,bx,c,alg,prsv);  A intmat, bx intvec, c intvec, alg string,
                            prsv intvec
solve_IP(A,bx,c,alg,prsv);  A intmat, bx list of intvec, c intvec, alg string,
                            prsv intvec
RETURN:  solution of the associated integer programming problem as explained
         in IP.lib
         return type = type of bx
EXAMPLE: example solve_IP;  shows an example"
{
  if(size(#)==4)
  {
    if(typeof(#[2])=="intvec")
    {
      return(solve_IP_1(#[1],#[2],#[3],#[4]));
    }
    else
    {
      return(solve_IP_2(#[1],#[2],#[3],#[4]));
    }
  }
  else
  {
    if(typeof(#[2])=="intvec")
    {
      return(solve_IP_3(#[1],#[2],#[3],#[4],#[5]));
    }
    else
    {
      return(solve_IP_4(#[1],#[2],#[3],#[4],#[5]));
    }
  }
}
example
{
  "EXAMPLE"; echo=2;

  // call with single right hand vector
  intmat A[2][3]=1,1,0,0,1,1;
  A;
  intvec b1=1,1;
  b1;
  intvec c=2,2,1;
  c;
  intvec solution_vector=solve_IP(A,b1,c,"pct");
  solution_vector;

  // call with list of right hand vectors
  intvec b2=-1,1;
  list l=b1,b2;
  l;
  list solution_list=solve_IP(A,l,c,"ct");
  solution_list;

  // call with single initial solution vector
  A=2,1,-1,-1,1,2;
  A;
  b1=3,4,5;
  solution_vector=solve_IP(A,b1,c,"du");

  // call with single initial solution vector and algorithm needing a positive
  // row space vector
  solution_vector=solve_IP(A,b1,c,"hs");

  // call with single initial solution vector and positive row space vector
  intvec prsv=1,2,1;
  prsv;
  solution_vector=solve_IP(A,b1,c,"hs",prsv);
  solution_vector;

  // call with list of initial solution vectors and positive row space vector
  b2=7,8,0;
  l=b1,b2;
  l;
  solution_list=solve_IP(A,l,c,"blr",prsv);
  solution_list;
}