/************************************************************************
* retroLP solves linear programs using the simplex method.             *
* Copyright (C) 2001, 2003, 2005 Gavriel Yarmish & Richard Van Slyke.  *
* For further information, and gpl licensing, see the end of this file.*
************************************************************************/


// New pivot (reversed) for column major format (8/02)
// 9/9/02 Added this new routine which does column choice
// in additon to pivoting. This saves a reading of the
// matrix from memory. If the return value is positive,
// the routine successfully successfully found a good
// pivot column with that index, if negative it didn't find
// one (actually, this routine does not set < 0; but outside
// of routines negative values means no candidate), if 0,
// there is no eligible column
// IMPORTANT! IT IS IMPORTANT TO MAINTAIN CONSISTANCY AMONG
// PIVOT (used in refresh), PIVOT_PLUS (used in regular iteration,
// AND CHOOSE_COLUMN (used if no available column from PIVOT_PLUS)

#include "simplex.h"
#include "params.h"
#include <math.h>
#include <stdlib.h>
#include <float.h>

extern long cntnfx, cntcc;
// To improve locality, we scan matrix alternately left to right
// and right to left, controlled by pdir
extern short pdir; // Used in pivot, and pivot+

long choose_row(LP_state *, long , short,
		double *, iter_result *, double);

long               // Proposed new pivot column
pivot_plus
(
 LP_state *lp,   // Linear programming data structure
 long ip,        // Index of pivot row
 long jp,        // Incoming index of pivot column
 double *val,    // Value of best column coming out
 short *up_flag, // If true we wish to increase the non-basic
 double delta           // EXPAND tolerance; used only for ccr=3
 )
{
  long i, j, jn = 0; // jn is candidate for column proposal out
  double t;
  double piv;
  short e_switch, // 0 if column is inelgible; +1 to ncrease;
    // -1 to decrease
    fix;
  
  double temp;
  double num, den, dx;
  iter_result nonbasicflag;
  piv=1.0/lp->a[jp][ip];
  
  switch(ccr)
    {
    case 1:
      //classical (Dantzig) column choice rule
      *val=epsc;
      break;
    case 2:
      //steepest edge pricing
      *val = -DBL_MAX;
      break;
    case 3:
      *val = -DBL_MAX;
      break;
    }
  //--------------- COLUMN 0 --------------------
  // Handle column 0 separately
  t = (lp->a[0][ip] *= -piv);
  for (i=0;i<=lp->m+1;i++)
    {
      if(i!=ip) lp->a[0][i] += lp->a[jp][i]*t;
    }
  //--------------- GENERAL COLUMN ---------------
  // General case
  // See which direction to go through the matrix
  long jj;
  pdir = -pdir;
  for (jj=1;jj<=lp->n; jj++)
    {
      if(pdir > 0)j=jj;else j= lp->n - jj +1;
      if(j==jp) continue; // Go to next j
                          // Do pivot column at end
      t = (lp->a[j][ip] *= -piv);
      // Check to see if column is eligible
      // Count non-fixed columns
      e_switch = 0;
      fix = (lp->nrange[j] == ATBOTH);
      if(!fix) cntnfx++;
      // Check to see if we may want to increase non-basic
      if((temp=(lp->a[j][lp->mm])+(lp->a[jp][lp->mm])*t) > epsc &&
	 lp->nrange[j]!=ATUPPER) e_switch = +1;
      // Check to see if we may want to decrease non=basic
      else if(temp < -epsc && lp->nrange[j]!=ATLOWER) e_switch = -1;
      if(fix || (e_switch == 0))
	//-------- FIXED OR INELIGIBLE COLUMN JUST PIVOT ---------
	// Only pivot do not evaluate column choice rule; not eligible
	{
	  if(!fix) cntcc++;  // Count ineligible columns
	  
	  for (i = 0; i <= lp->m+1; i++)
	    {
	      lp->a[j][i] += lp->a[jp][i]*t;
	    }
	  lp->a[j][ip] = t;
	  continue; // skip to next j
	}
      // Otherwise pivot AND evaluate for ccr
      switch(ccr)
	{
	case 1:
	case 3:
	  // For classical and greatest change column pivot doesn't change
	  for (i = 0; i <= lp->m+1; i++)
	    {
	      lp->a[j][i] += lp->a[jp][i]*t;
	    }
	  lp->a[j][ip] = t;
	  break;
	  // For steepest edge we accumulate vector length
	case 2:
	  den = 1.0;
	  lp->a[j][0] += lp->a[jp][0]*t; // handle rows 0 and m+1
	  // separately
	  for (i=1;i<=lp->m;i++)
	    {
	      if(i!=ip) lp->a[j][i] += lp->a[jp][i]*t;
	      den+=lp->a[j][i]*lp->a[j][i];
	    }
	  lp->a[j][lp->m+1] += lp->a[jp][lp->m+1]*t;
	  break;
	}
      // after i iteration
      switch(ccr)
	{
	case 1:
	  if(temp > *val && e_switch == +1)
	    {
	      *val=temp;
	      jn=j;
	      *up_flag = TRUE;
	    }
	  if(-temp > *val && e_switch == -1)
	    {
	      *val=-temp;
	      jn=j;
	      *up_flag = FALSE;
	    }
	  break;
	case 2:
	  // Column is now updated look at column choice value
	  {
	    num = temp * temp;
	    if (num/den > *val)
	      {
		*val=num/den;
		jn=j;
		if(e_switch == 1)*up_flag = TRUE;
		else *up_flag = FALSE;
	      }
	  }
	  break;
	case 3:
	  if(e_switch == +1) *up_flag = TRUE;
	  else *up_flag = FALSE;
	  i = choose_row( lp, j, *up_flag, &dx, &nonbasicflag, delta);
	  if(i < 0)
	    {
	      *val = DBL_MAX;
	      return j; // Unbounded
	    }
	  if (lp->a[j][lp->mm]*dx > *val)
	    {
	      *val=lp->a[j][lp->mm]*dx;
	      jn=j;
	    }
	  break;
	}
    } // end of for loop (j)
  
  //--------------- UPDATE PIVOT COLUMN --------------------
  // Ugh! Have to check ccr for pivot column too for Phase 1
  e_switch = 0;
  if(lp->nrange[jp] == ATBOTH)fix = TRUE; else fix = FALSE;
  if(!fix) cntnfx++;
  // Check to see if we may want to increase non-basic
  if(((temp=lp->a[jp][lp->mm]* piv) > epsc) &&
     (lp->nrange[jp]!=ATUPPER)) e_switch = +1;
  // Check to see if we may want to decrease non=basic
  else if((temp < -epsc) && (lp->nrange[jp]!=ATLOWER)) e_switch = -1;
  if(fix || (e_switch == 0))
    //--- PIVOT COL. FIXED OR INELIGIBLE COLUMN JUST PIVOT ---
    {
      if(!fix) cntcc++;  // Count ineligible columns
      
      for (i=0;i<=lp->m+1;i++)
	{
	  lp->a[jp][i] *= piv;
	}
      lp->a[jp][ip]=piv;
    }
  else
    {
      assert((e_switch == -1) || (e_switch == +1));
      // Otherwise pivot AND evaluate for ccr
      switch(ccr)
	{
	case 1:
	case 3:
	  // For classical and greatest change column choice pivot
	  // unaffected
	  for (i=0;i<=lp->m+1;i++)
	    {
	      lp->a[jp][i] *= piv;
	    }
	  lp->a[jp][ip]=piv;
	  // For steepest edge we accumulate vector length
	  break;
	case 2:
	  den = 1.0;
	  lp->a[jp][0] *= piv; // handle rows 0, ip, and m+1
	  // separately
	  lp->a[jp][ip] = piv;
	  for (i=1;i<=lp->m;i++)
	    {
	      if(i!=ip) lp->a[jp][i] *= piv;
	      den+=lp->a[jp][i]*lp->a[jp][i];
	    }
	  lp->a[jp][lp->m+1] *= piv;
	  break;
	}
      // after i iteration
      switch(ccr)
	{
	case 1:
	  if(temp > *val && e_switch == +1)
	    {
	      *val=temp;
	      jn=jp;
	      *up_flag = TRUE;
	    }
	  if(-temp > *val && e_switch == -1)
	    {
	      *val=-temp;
	      jn=jp;
	      *up_flag = FALSE;
	    }
	  break;
	case 2:
	  // Column is now updated look at column choice value
	  num = temp * temp;
	  if (num/den > *val)
	    {
	      *val=num/den;
	      jn=jp;
	      if(e_switch == 1)*up_flag = TRUE;
	      else *up_flag = FALSE;
	    }
	  break;
	case 3:
	  if(e_switch == +1) *up_flag = TRUE;
	  else *up_flag = FALSE;
	  i = choose_row( lp, jp, *up_flag, &dx, &nonbasicflag, delta);
	  if(i < 0)
	    {
	      *val = DBL_MAX;
	      return jp; // Unbounded
	    }
	  if (lp->a[jp][lp->mm]*dx > *val)
	    {
	      *val=lp->a[jp][lp->mm]*dx;
	      jn=jp;
	    }
	  break;
	}
    } // end of check for jp
  
  
  //--------------- CHOOSE COLUMN -------------------------
  // Check to see if any column selected
  if(!jn)*val=0.0;
  return jn;
  // This should not be reached
  assert(0);
  return -2;
  
}

/************************************************************************
* retroLP solves linear programs using the simplex method.             *
* Copyright (C) 2001, 2003, 2005 Gavriel Yarmish & Richard Van Slyke.  *
*                                                                      *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or    *
* (at your option) any later version.                                  *
*                                                                      *
* This program is distributed in the hope that it will be useful,      *
* but WITHOUT ANY WARRANTY; without even the implied warranty of       *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
* GNU General Public License for more details.                         *
*                                                                      *
* You should have received a copy of the GNU General Public License    *  
* along with this program; if not, write to the Free Software          *
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,           *
*                                                                      *
************************************************************************/
