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

version="$Id: intprog.lib,v 1.1 2000-09-14 12:37:09 Singular Exp $";

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

info="
LIBRARY: intprog.lib      INTEGER PROGRAMMING WITH GROEBNER BASIS METHODS

AUTHOR:  Christine Theis, email: ctheis@math.uni-sb.de

PROCEDURES:

solve_IP(intmat A, [intvec|list] bx, intvec c, string alg [, intvec prsv]);
    procedures for solving IP-problems
";

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



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: The number of matrix columns does not equal the size of the cost vector.";
    return(v);
  }

  // create first temporary file with which 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: The chosen algorithm needs a 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: The number of matrix columns does not equal the size of the cost vector.";
    return(l);
  }

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

  // create first temporary file with which 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: The chosen algorithm needs a 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: The number of matrix columns does not equal the size of the cost vector.";
    return(v);
  }

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

  // create first temporary file with which 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: The number of matrix columns does not equal the size of the cost vector.";
    return(l);
  }

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

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

  // create first temporary file with which 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:   type of bx: solution of the associated integer programming problem(s) as explained in
@texinfo
@ref{Toric ideals and integer programming}.
@end texinfo
NOTE:     This procedure returns the solution(s) of the given IP-problem(s) or the message `not solvable'.

One may call the procedure with several different algorithms:

- the algorithm of Conti/Traverso using elimination (ect),

- the algorithm of Pottier (pt),

- an algorithm of Bigatti/La Scala/Robbiano (blr),

- the algorithm of Hosten/Sturmfels (hs),

- the algorithm of DiBiase/Urbanke (du).

The argument `alg' should be the abbreviation for an algorithm as above: ect, pt, blr, hs or du.

`ct' allows to 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 may only be applied if A and b have nonnegative entries. All other algorithms need initial solutions of the IP-problem.

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 coefficients in the row space of A. If no row of A contains only positive entries, one has to use the versions of solve_IP which take such a vector prsv as an argument.

solve_IP may also be called with a list bx of intvecs instead of a single intvec.
EXAMPLE:  example solve_IP; shows an example
SEE ALSO: intprog_lib, toric_lib, Integer programming
"
{
  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;
}