///////////////////////////////////////////////////////////////////////////////
version="$Id$";
category="Commutative Algebra";
info="
LIBRARY: intprog.lib      Integer Programming with Groebner Basis Methods
AUTHOR:  Christine Theis, email: ctheis@math.uni-sb.de

PROCEDURES:
 solve_IP(..);   procedures for solving Integer Programming  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: same type as 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 (ct),
        @*- the positive variant of the algorithm of Conti/Traverso (pct),
        @*- 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: ct, pct, ect, pt, blr, hs or du.

        `ct' allows computation of 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.
SEE ALSO: intprog_lib, toric_lib, Integer programming
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;
  // 1. call with single right-hand vector
  intmat A[2][3]=1,1,0,0,1,1;
  intvec b1=1,1;
  intvec c=2,2,1;
  intvec solution_vector=solve_IP(A,b1,c,"pct");
  solution_vector;"";

  // 2. 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;"";

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

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

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

  // 6. 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;
}