//  Implementation of the class linearProgram
//
// This file contains procedures to solve a Linear Programming
// problem of n variables and m constraints, the last p of which
// are equality constraints by the dual simplex method.
//
// The code was originally in Algol 60 and was published in Collected
// Algorithms from CACM (alg #333) by Rudolfo C. Salazar and Subrata
// K. Sen in 1968 under the title MINIT (MINimum ITerations) algorithm
// for Linear Programming.
// It was directly translated into C by Badri Lokanathan, Dept. of EE,
// University of Rochester, with no modification to code structure. 
//
// The problem statement is
// Maximise z = cX
//
// subject to
// AX <= b
// X >=0
//
// c is a (1*n) row vector
// A is a (m*n) matrix
// b is a (m*1) column vector.
// e is a matrix with with (m+1) rows and lcol columns (lcol = m+n-p+1)
//   and forms the initial tableau of the algorithm.
// td is read into the procedure and should be a very small number,
//   say 10**-8
// 
// The condition of optimality is the non-negativity of
// e(1,j) for j = 1 ... lcol-1 and of e(1,lcol) = 2 ... m+1.
// If the e(i,j) values are greater than or equal to -td they are
// considered to be non-negative. The value of td should reflect the 
// magnitude of the co-efficient matrix.


#include "linearProgram.h"



#ifndef LARGENUMBER
  #define LARGENUMBER 1e6
#endif  

#ifndef VLARGENUMBER
  #define VLARGENUMBER 1e8
#endif  

#ifndef VSMALLNUMBER
  #define VSMALLNUMBER 1e-8
#endif  

#define ABS(A) (A > 0 ? A: -A)
#define MAX(A,B) (A > B ? A : B)
#define MIN(A,B) (A < B ? A : B)



Matrix<double> nullD(0);
Matrix<int> nullI(0);



//////////////////
// Constructors //
//////////////////

linearProgram::linearProgram () :
   k(0), m(0), mL(0), mG(0), mE(0), n(0), lcol(0), 
   L(0), im(0), jmin(0), jm(0), imax(0), debug(0), 
   maxim(true),
   td(VSMALLNUMBER), // This is application specific.
   gmin(0.0), phimax(0.0),
   imin(nullI), jmax(nullI), ind(nullI), ind1(nullI), chk(nullI), 
   permRow(nullI), e(nullD), thmin(nullD), delmax(nullD) 
{
}



void linearProgram::newLP (
   int objSense,
   Array<double>& c,
   Array<double>& A,
   Array<char>& R,
   Array<double>& b)
{
   // Check sizes of arrays
   if ( c.n()!=A.n() || R.n()!=A.m() || b.n()!=A.m() )
      throw whereDimMismatch(
	 "void linearProgram::newLP (objSense, Array<double> c,A,R,b)",
	 A.m(),A.n(),b.n(),c.n());

   // Update global sizes
   m = A.m();
   n = A.n();
   mL = mG = mE = 0;
   for (int i=0; i<R.n(); i++)
      if	(R(i) == 'L')	mL++;
      else if	(R(i) == 'G')	mG++;
      else			mE++;
   lcol = n + mL + mG + 1;
   maxim = (objSense==MAXIMIZE);

   // Allocate or reallocate the memory and initialize
   jmax.resize(1,m+1,0);
   ind1.resize(1,m+1,0);
   chk.resize(1,m+1,0);
   delmax.resize(1,m+1,0.0);
   imin.resize(1,lcol,0);
   ind.resize(1,lcol,0);
   thmin.resize(1,lcol,0.0);
   e.resize(m+1,lcol,0.0);
   permRow.resize(1,m,0);

   // Construct the permutation for the rows (<=, >=, ==)
   int imL=0; int imG=mL; int imE=mL+mG;
   for (int row=0; row<R.n(); row++)
      if	( R(row)=='L' )	permRow(imL++)=row;
      else if	( R(row)=='G')	permRow(imG++)=row;
      else			permRow(imE++)=row;   

   // Generate initial tableau.
   e.s(0,1,0,n) = (-objSense) * c;

   for(i = 0; i < mL; i++)
   {   
      e.s(i+1,1,0,n) = A.s(permRow(i));
      e(i+1,lcol-1) = b(permRow(i));
   }
   for(i = mL; i < mL+mG; i++)
   {
      e.s(i+1,1,0,n) = -A.s(permRow(i));
      e(i+1,lcol-1) = -b(permRow(i));
   }
   for(i = mL+mG; i < m; i++)
   {
      e.s(i+1,1,0,n) = A.s(permRow(i));
      e(i+1,lcol-1) = b(permRow(i));
   }

   for(i = 0; i < mL+mG; i++)
      e(1+i,n+i) = 1.0;
}
  


void linearProgram::newLP (
   int objSense,
   Array<int>& c,
   Array<int>& A,
   Array<char>& R,
   Array<int>& b)
{
   register int i,j;
   // Check sizes of arrays
   if ( c.n()!=A.n() || R.n()!=A.m() || b.n()!=A.m() )
      throw whereDimMismatch(
	 "void linearProgram::newLP (objSense, Array<double> c,A,R,b)",
	 A.m(),A.n(),b.n(),c.n());

   // Update global sizes
   m = A.m();
   n = A.n();
   mL = mG = mE = 0;
   for (i=0; i<R.n(); i++)
      if	(R(i) == 'L')	mL++;
      else if	(R(i) == 'G')	mG++;
      else			mE++;
   lcol = n + mL + mG + 1;
   maxim = (objSense==MAXIMIZE);

   // Allocate or reallocate the memory and initialize
   jmax.resize(1,m+1,0);
   ind1.resize(1,m+1,0);
   chk.resize(1,m+1,0);
   delmax.resize(1,m+1,0.0);
   imin.resize(1,lcol,0);
   ind.resize(1,lcol,0);
   thmin.resize(1,lcol,0.0);
   e.resize(m+1,lcol,0.0);
   permRow.resize(1,m,0);

   // Construct the permutation for the rows (<=, >=, ==)
   int imL=0; int imG=mL; int imE=mL+mG;
   for (int row=0; row<R.n(); row++)
      if	( R(row)=='L' )	permRow(imL++)=row;
      else if	( R(row)=='G')	permRow(imG++)=row;
      else			permRow(imE++)=row;   

   // Generate initial tableau.
   for(i=0; i<n; i++)
      e(0,i) = double((-objSense) * c(i));

   for(i=0; i<mL; i++)
   {   
      for (j=0; j<n; j++)
	 e(i+1,j) = double(A(permRow(i),j));
      e(i+1,lcol-1) = double(b(permRow(i)));
   }
   for(i=mL; i<mL+mG; i++)
   {
      for (j=0; j<n; j++)
	 e(i+1,j) = double(-A(permRow(i),j));
      e(i+1,lcol-1) = double(-b(permRow(i)));
   }
   for(i=mL+mG; i<m; i++)
   {
      for (j=0; j<n; j++)
	 e(i+1,j) = double(A(permRow(i),j));
      e(i+1,lcol-1) = double(b(permRow(i)));
   }

   for(i=0; i<mL+mG; i++)
      e(1+i,n+i) = 1.0;
}



/*inline*/ linearProgram::linearProgram (
   int objSense,
   Array<double>& c,
   Array<double>& A,
   Array<char>& R,
   Array<double>& b) :
   k(0), m(0), mL(0), mG(0), mE(0), n(0), lcol(0), 
   L(0), im(0), jmin(0), jm(0), imax(0), debug(0), 
   maxim(true),
   td(VSMALLNUMBER), // This is application specific.
   gmin(0.0), phimax(0.0),
   imin(nullI), jmax(nullI), ind(nullI), ind1(nullI), chk(nullI), 
   permRow(nullI), e(nullD), thmin(nullD), delmax(nullD) 
{
   newLP(objSense,c,A,R,b);
} 



/*inline*/ linearProgram::linearProgram (
   Array<double>& c,
   Array<double>& A,
   Array<double>& b) :
   k(0), m(0), mL(0), mG(0), mE(0), n(0), lcol(0), 
   L(0), im(0), jmin(0), jm(0), imax(0), debug(0), 
   maxim(true),
   td(VSMALLNUMBER), // This is application specific.
   gmin(0.0), phimax(0.0),
   imin(nullI), jmax(nullI), ind(nullI), ind1(nullI), chk(nullI), 
   permRow(nullI), e(nullD), thmin(nullD), delmax(nullD) 
{
   Matrix<char> R(1,A.m(),'L');
   newLP(MAXIMIZE,c,A,R,b);
} 



/*inline*/ linearProgram::linearProgram (
   int objSense,
   Array<int>& c,
   Array<int>& A,
   Array<char>& R,
   Array<int>& b) :
   k(0), m(0), mL(0), mG(0), mE(0), n(0), lcol(0), 
   L(0), im(0), jmin(0), jm(0), imax(0), debug(0), 
   maxim(true),
   td(VSMALLNUMBER), // This is application specific.
   gmin(0.0), phimax(0.0),
   imin(nullI), jmax(nullI), ind(nullI), ind1(nullI), chk(nullI), 
   permRow(nullI), e(nullD), thmin(nullD), delmax(nullD) 
{
   newLP(objSense,c,A,R,b);
} 



/*inline*/ linearProgram::linearProgram (
   Array<int>& c,
   Array<int>& A,
   Array<int>& b) :
   k(0), m(0), mL(0), mG(0), mE(0), n(0), lcol(0), 
   L(0), im(0), jmin(0), jm(0), imax(0), debug(0), 
   maxim(true),
   td(VSMALLNUMBER), // This is application specific.
   gmin(0.0), phimax(0.0),
   imin(nullI), jmax(nullI), ind(nullI), ind1(nullI), chk(nullI), 
   permRow(nullI), e(nullD), thmin(nullD), delmax(nullD) 
{
   Matrix<char> R(1,A.m(),'L');
   newLP(MAXIMIZE,c,A,R,b);
} 



///////////////
// Operators //
///////////////

int linearProgram::solve (
   double& optValue,
   Matrix<double>& primal,
   Matrix<double>& dual)
{
  register int i,j;
   /* Now begins the actual algorithm. */
   for(i = 1; i < m+1; i++)
      chk(i) = -1;	/* Indexing problem; Algol
			 * version treates 0 as out
			 * of bounds, in C we prefer -1.
			 */

   if(!mE)	/* There are no equality constraints */
      goto RCS;
   else
   {
      if(phase1())
	 return(1);
   }
		
  RCS: L = 1; k = 1;

  if(debug)
     tab();	/* Print the tableau. */

  for(j = 0; j < lcol - 1; j++)
  {
     if(e(0,j) < -td)
	ind((L++)-1) = j; /* ind(L) keeps track of the
			   * columns in which e(0,j)
			   * is negative.
			   */
  }

  for(i = 1; i < m + 1; i++)
  {
     if(e(i,lcol-1) < -td)
	ind1((k++)-1) = i; /* ind1(k) keeps track of the
			    * rows in which e(i,lcol)
			    * is negative.
			    */
  }

  if(L == 1)
  {
     if(k == 1) /* Results */
	goto RESULTS;
     else
     {
	if(k == 2)
	{
	   for(j = 0; j < lcol - 1; j++)
	   {
	      if(e(ind1(0),j) < 0)
		 goto R;
	   }
	   /* Primal problem has no feasible solutions. */
	   return(2);
	}
	else
	   goto R;
     }
  }
  else
  {
     if(L == 2)
     {
	if(k == 1)
	{
	   for(i = 1; i < m + 1; i++)
	   {
	      if(e(i,ind(0)) > 0)
		 goto C;
	   }

	   /* Primal problem is unbounded. */
	   return(3);
	}
	else
	   goto S;
     }

     if(k == 1)
	goto C;
     else
	goto S;
  }

  R:
  prophi();
  if(rowtrans(imax,jm))
     return(1);
  goto RCS;

  C:
  progamma();
  if(rowtrans(im,jmin))
     return(1);
  goto RCS;

  S:
  progamma();
  prophi();
  if(gmin == LARGENUMBER)
  {
     if(rowtrans(imax,jm))
	return(1);
     goto RCS;
  }

  if(phimax == - LARGENUMBER)
  {
     if(rowtrans(im,jmin))
	return(1);
     goto RCS;
  }

  if(ABS(phimax) > ABS(gmin))
  {
     if(rowtrans(imax,jm))
	return(1);
  }
  else
  {
     if(rowtrans(im,jmin))
	return(1);
  }
  goto RCS;

  RESULTS:
  /* Output results here. */
  optValue = e(0,lcol-1);
  primal.resize(1,n,0.0);
  dual.resize(1,m,0.0);

  for(i = 1; i < m + 1; i++)
  {
     if(chk(i) >= n)
	chk(i) = -1;

     if(chk(i) > -1)
	primal(chk(i)) = e(i,lcol-1);
  }

  for(j = n; j < lcol - 1; j++)
     dual(permRow(j-n)) = e(0,j);
	
  return(0);
}



/*inline*/ int linearProgram::solve (
   double& optValue,
   Matrix<double>& primal)
{
   Matrix<double> dual;
   return solve(optValue,primal,dual);
}



/*inline*/ int linearProgram::solve (
   double& optValue)
{
   Matrix<double> primal;
   return solve(optValue,primal);
}



ostream& operator << (ostream& s, const linearProgram& lp)
{
   return s << endl << lp.e << endl;
}



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


// Performs the usual tableau transformations in a linear programming
// problem, (p,q) being the pivotal element.
// Returns the following error codes:
// 0: Everything was OK.
// 1: No solution.

int linearProgram::rowtrans(int p,int q)
{
   register int i, j;
   float dummy;

   if(p == -1 || q == -1) /* No solution. */
      return(1);

   dummy = e(p,q);

   if(debug)
      cout<< "--\nPivot element is e(" << p << "," << q << ") = " 
	  << dummy << endl;

   for(j = 0; j < lcol; j++)
      e(p,j) /= dummy; 

   for(i = 0; i < m + 1; i++)
   {
      if(i != p)
      {
	 if(e(i,q) != 0.0)
	 {
	    dummy = e(i,q);
	    for(j = 0; j < lcol; j++)
	       e(i,j) -= e(p,j) * dummy;
	 }
      }
   }

   chk(p) = q;
   return(0);
}



// Performs calculations over columns to determine the pivot element.
void linearProgram::progamma()
{
   float theta, gamma;
   register int i, L1;

   gmin = LARGENUMBER;	/* gmin is set to a large no. for
			 * initialization purposes.
			 */
   jmin = -1;		/* Array indices in C start from 0 */

   for(L1 = 0; L1 < L - 1; L1++)
   {
      imin(ind(L1)) = 0;
      thmin(ind(L1)) = LARGENUMBER;
      for(i = 1; i < m + 1; i++)
      {
	 if(e(i,ind(L1)) > td && e(i,lcol-1) >= -td)
	 {
	    theta = e(i,lcol-1) / e(i,ind(L1));
	    if(theta < thmin(ind(L1)))
	    {
	       thmin(ind(L1)) = theta;
	       imin(ind(L1)) = i;
	    }
	 }
      }

      if(thmin(ind(L1)) == LARGENUMBER)
	 gamma = VLARGENUMBER;
      else
	 gamma = thmin(ind(L1)) * e(0,ind(L1));

      if(gamma < gmin)
      {
	 gmin = gamma;
	 jmin = ind(L1);
      }
   }
   if(jmin > -1)
      im = imin(jmin);
}



// Performs calculations over rows to determine the pivot element.
void linearProgram::prophi()
{
   float delta, phi;
   register int j, k1;

   phimax = - LARGENUMBER;	/* phimax is set to a small no. for
				 * initialization purposes.
				 */
   imax = -1;		/* Array indices in C start from 0 */
	
   for(k1 = 0; k1 < k - 1; k1++)
   {
      jmax(ind1(k1)) = 0;
      delmax(ind1(k1)) = - LARGENUMBER;
      for(j = 0; j < lcol - 1; j++)
      {
	 if(e(ind1(k1),j) < -td && e(0,j) >= -td)
	 {
	    delta = e(0,j) / e(ind1(k1),j);
	    if(delta > delmax(ind1(k1)))
	    {
	       delmax(ind1(k1)) = delta;
	       jmax(ind1(k1)) = j;
	    }
	 }
      }

      if(delmax(ind1(k1)) == - LARGENUMBER)
	 phi = - VLARGENUMBER;
      else
	 phi = delmax(ind1(k1)) * e(ind1(k1),lcol-1);
		
      if(phi > phimax)
      {
	 phimax = phi;
	 imax = ind1(k1);
      }
   }
   if(imax > -1)
      jm = jmax(imax);
}



// Applied only to equality constraints if any.
int linearProgram::phase1()
{
   float theta, gamma;
   register int i, j, r;
   /* Fix suggested by Holmgren, Obradovic, Kolm. */
   register int im1, jmin1, first;

   im1 = jmin1 = -1;
   /* Fix suggested by Messham to allow negative coeffs. in
    * equality constraints.
    */
   for(i = m - mE + 1; i < m + 1; i++)
   {
      if(e(i,lcol - 1) < 0)
      {
	 for(j = 0; j < lcol; j++)
	    e(i,j) = -e(i,j);
      }
   }

   for(r = 0; r < mE; r++)
   {
      gmin = LARGENUMBER;	/* gmin is set to a large no. for
				 * initialization purposes.
				 */
      L = 1;
      jmin = -1;
      first = 1;
      for(j = 0; j < n; j++)
      {
	 thmin(j) = LARGENUMBER;
	 /* Fix suggested by Kolm and Dahlstrand */
	 /* if(e(0,j) < 0) */
	 if(e(0,j) < -td)
	    ind((L++)-1) = j;
      }

     L1:	
     if(L == 1)
     {
	for(j = 0; j < n; j++)
	   ind(j) = j;

	L = n + 1;
     }

     for(k = 0; k < L - 1; k++)
     {
	for(i = m - mE + 1; i < m + 1; i++)
	   if(chk(i) == -1)
	   {
	      /* Fix suggested by Kolm
	       * and Dahlstrand
	       */
	      /* if(e(i,ind(k)) > 0.0) */
	      if(e(i,ind(k)) > td)
	      {
		 if((theta = e(i,lcol-1) /
		     e(i,ind(k))) < thmin(ind(k)))
		 {
		    thmin(ind(k)) = theta;
		    imin(ind(k)) = i;
		 }
	      }
	   }

	/* Fix suggested by Obradovic overrides
	 * fixes suggested by Kolm and Dahstrand
	 * as well as Messham.
	 */
	if(thmin(ind(k)) < LARGENUMBER)
	{
	   gamma = thmin(ind(k)) * e(0,ind(k));
	   if(gamma < gmin)
	   {
	      gmin = gamma;
	      jmin = ind(k);
	   }
	}
     }
     if(jmin == -1)
     {
	if(first)
	{
	   first = 0;
	   L = 1;
	   goto L1;
	}
	else
	   im = -1;
     }
     else
	im = imin(jmin);

     if(im == im1 && jmin == jmin1)
     {
	L = 1;
	goto L1;
     }

     if(debug)
	tab();	/* Print the tableau. */

     if(rowtrans(im,jmin))
	return(1);
		
     im1 = im;
     jmin1 = jmin;
   }
   return(0);
}



// The following procedure is for debugging. It simply prints the
// current tableau.
void linearProgram::tab()
{
   register int i, j;
   cout << endl;
   for(i = 0; i < 35; i++) cout << "-";
   cout << " TABLEAU ";
   for(i = 0; i < 35; i++) cout << "-";
   cout << endl;
   for(i = 0; i < m+1; i++)
   {
      for(j = 0; j < lcol; j++)
	 cout << setw(10) << setprecision(3) << e(i,j);
      cout << endl;
   }
}