/************************************************************************
* 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.*
************************************************************************/

#include "simplex.h"
#include "params.h"
#include <math.h> //for fabs()
#include <float.h>
#include <stdio.h>
#include <stdlib.h>

// This row choice rule uses the EXPAND methodology (Gill et al, 1989)
long choose_row
(
  LP_state *lp,              // LP data structure
  long jp,                   // Pivot column
  short up_flag,             // True if non-basic variable increasing
  double *dx,                // Change in non-basic variable
  iter_result *nonbasicflag, // Type of state change
  double d                   // Relaxation delta from EXPAND
)
{
    double  tempmin, t, alphamin, pivot = 0.0;
    long ip, i;
    iter_result tempflag;
// If ip = -1 is returned we have unbounded problem; ip = 0 if non-basic
// goes to one of its bounds; otherwise ip = pivot row.
    ip= -1;
    double delta=DBL_MAX; // Change in non-basic

/////////////////////////////////////////////////////////////////////////////////
// PASS I (RELAXATION)
/////////////////////////////////////////////////////////////////////////////////

// In Pass I we don't need to calculate nonbasicflag, tempflag and/or
// ip. Leave in for now for debugging. Take out later.

    *nonbasicflag=UNBOUNDED;

// HOW MUCH CAN NON-BASIC VARIABLE MOVE CONSTRAINED ONLY BY ITS BOUNDS?

    if(   up_flag
       &&( (lp->ntype[jp]==UBOUNDED)
          ||(lp->ntype[jp]==BOTHBOUNDS)))
// Trying to increase non-basic variable with upper bound
    {
       switch(lp->nrange[jp])
       {
       case BELOW:
          assert(lp->x[jp] < lp->nl[jp] + d);
          delta=lp->nl[jp]-lp->x[jp] + d;
          assert(delta > 0);
          *nonbasicflag=NONBASICTOFEASIBLEL;
          break;
       case ATLOWER:
          assert(fabs(lp->x[jp] - lp->nl[jp]) < d);
       // Fall through
       case BETWEEN:
          assert((lp->x[jp] > lp->nl[jp] - d)&&(lp->x[jp] < lp->nu[jp] + d));
          delta=lp->nu[jp]-lp->x[jp] + d;
          assert(delta > 0);
          *nonbasicflag=NOPIVOTUBOUND;
          break;
       case ABOVE:
          assert(lp->x[jp] > lp->nu[jp] - d);
          delta=DBL_MAX;
          *nonbasicflag=UNBOUNDED;
          break;
       // These cases shouldn't happen
       default: // Including ATUPPER and ATBOTH
          assert((lp->nrange[jp]!=ATUPPER) && (lp->nrange[jp]!=ATBOTH));
          printf("Incorrect call to row_choice\n");
          exit(501);
       }
       ip= 0;
    }
    else if (    !up_flag
              && ((   lp->ntype[jp]==LBOUNDED)
                   ||(lp->ntype[jp]==BOTHBOUNDS)))
// Trying to decrease non-basic variable with lower bound
    {
       switch(lp->nrange[jp])
       {
       case ABOVE:
          assert(lp->x[jp] > lp->nu[jp] - d);
          delta=lp->x[jp]-lp->nu[jp] + d;
          assert(delta > 0);
          *nonbasicflag=NONBASICTOFEASIBLEU;
          break;
       case ATUPPER:
          assert(fabs(lp->x[jp] - lp->nu[jp]) < d);
       // Fall through
       case BETWEEN:
          assert((lp->x[jp] > lp->nl[jp] - d)&&(lp->x[jp] < lp->nu[jp] + d));
          delta=lp->x[jp]-lp->nl[jp] + d;
          assert(delta > 0);
          *nonbasicflag=NOPIVOTLBOUND;
          break;
       case BELOW:
          assert(lp->x[jp] < lp->nl[jp] + d);
          delta=DBL_MAX;
          *nonbasicflag=UNBOUNDED;
          break;
       // These cases shouldn't happen
       default: // Including ATLOWER and ATBOTH
          assert((lp->nrange[jp]!=ATLOWER) && (lp->nrange[jp]!=ATBOTH));
          printf("Incorrect call to row_choice\n");
          exit(502);
       }
       ip=0;
    }
// WHAT LIMITS DO THE BASIC VARIABLES IMPOSE?
// There are three cases:
//    Basic value is increasing
//    Basic value is decreasing
//    Matrix element too small to pivot on.
        for (i=1;i<=lp->m;i++)  // Go through rows
        {
// INCREASING BASIC VALUE
        // Note: Basic variables should not be ATLOWER or ATUPPER
            assert((lp->brange[i]!=ATLOWER) && (lp->brange[i]!=ATUPPER));
            if(
                  (!up_flag && lp->a[jp][i]< - epsp)
               ||  ( up_flag && lp->a[jp][i]>   epsp)
              )
           {
              switch(lp->brange[i])
              {
              case BELOW:
                 assert(lp->b[i] < lp->bl[i] + d);
                 tempmin=(lp->bl[i]-lp->b[i] + d)/fabs(lp->a[jp][i]);
                 assert(tempmin > 0);
                 tempflag=BASICTOFEASIBLEL;
                 break;
              case BETWEEN:
                 assert((lp->b[i] > lp->bl[i] - d)&&(lp->b[i] < lp->bu[i] + d));
                 if(lp->btype[i] == LBOUNDED||lp->btype[i] == FREE)
                 {
                    tempmin=DBL_MAX;
                    tempflag=UNBOUNDED;
                    break;
                 }
                 else
                 {
                    tempmin=(lp->bu[i]-lp->b[i] + d)/fabs(lp->a[jp][i]);
                    assert(tempmin > 0);
                    tempflag=PIVOTUBOUND;
                    break;
                 }
              case ATBOTH:
                assert(fabs(lp->bl[i] - lp->bu[i])<d);
                assert(fabs(lp->b[i] - lp->bl[i]) < d);
                tempflag=PIVOTUBOUND;
                tempmin = (lp->bu[i]-lp->b[i] + d)/fabs(lp->a[jp][i]);
                break;
              case ABOVE:
                assert(lp->b[i] > lp->bu[i] - d);
                tempmin=DBL_MAX;
                tempflag=UNBOUNDED;
                break;
              default:
                printf("\nImproper result code");
                exit (302);
              }
              if(tempmin < delta)
                {
                   delta = tempmin;
                   ip = i;
                   *nonbasicflag = tempflag;
                }
          }

// DECREASING BASIC VALUE
        // Note: Basic variables should not be ATLOWER or ATUPPER
           else if(
                      ( up_flag && lp->a[jp][i] < -epsp)
                   ||  (!up_flag && lp->a[jp][i] >  epsp)
                  )
           {
              assert((lp->brange[i]!=ATLOWER) && (lp->brange[i]!=ATUPPER));
              switch(lp->brange[i])
              {
              case ABOVE:
                 assert(lp->b[i] > lp->bu[i] - d);
                 tempmin=(lp->b[i] - lp->bu[i] + d)/fabs(lp->a[jp][i]);
                 assert(tempmin > 0);
                 tempflag=BASICTOFEASIBLEU;
                 break;
              case BETWEEN:
                 assert((lp->b[i] > lp->bl[i] - d)&&(lp->b[i] < lp->bu[i] + d));
                 if(lp->btype[i] == BOTHBOUNDS||lp->btype[i] == LBOUNDED)
                 {
                    tempmin=(lp->b[i] - lp->bl[i] + d)/fabs(lp->a[jp][i]);
                    assert(tempmin > 0);
                    tempflag=PIVOTLBOUND;
                    break;
                 }
            
                 else
                 {
                    tempmin=DBL_MAX;
                    tempflag=UNBOUNDED;
                    break;
                 }
              case ATBOTH:
                 assert(fabs(lp->bl[i] - lp->bu[i])<d);
                 assert(fabs(lp->b[i] - lp->bl[i]) < d);
                 tempflag=PIVOTLBOUND;
                 tempmin=(lp->b[i] - lp->bl[i] + d)/fabs(lp->a[jp][i]);
                 break;
              case BELOW:
                 assert(lp->b[i] < lp->bl[i] + d);
                 tempflag=UNBOUNDED;
                 tempmin=DBL_MAX;
                 break;
              default:
                 printf("\nImproper result code");
                 exit (303);

              }
              if(tempmin < delta)
                 {
                    delta = tempmin;
                    ip = i;
                    *nonbasicflag = tempflag;
                 }
          } 
//  If a[i][j] is smaller than epsp zero it out
          else lp->a[jp][i]=0.0;
        }

/////////////////////////////////////////////////////////////////////////////////
// PASS II
/////////////////////////////////////////////////////////////////////////////////

    pivot = 0.0;
    *nonbasicflag=UNBOUNDED;
    tempmin = DBL_MAX;
// delta is the (non-negative) maximum amount that the non-basic
// variable is allowed to change. We now look for largest pivot element..
    assert(delta > 0);

// CONSTRAINED BY BOUNDS?

    if(   up_flag
       &&( (lp->ntype[jp]==UBOUNDED)
          ||(lp->ntype[jp]==BOTHBOUNDS)))
// Trying to increase non-basic variable with upper bound
    {
       switch(lp->nrange[jp])
       {
       case BELOW:
          if((t=(lp->nl[jp]-lp->x[jp])) <= delta)
          {
            pivot=DBL_MAX;
            tempmin = t;
            *nonbasicflag=NONBASICTOFEASIBLEL;
          }
          break;
       case ATLOWER:
       // Fall through
       case BETWEEN:
          if((t=(lp->nu[jp]-lp->x[jp])) <= delta)
          { 
            pivot=DBL_MAX;
            tempmin = t;
            *nonbasicflag=NOPIVOTUBOUND;
          }
          break;
       case ABOVE:
          pivot=DBL_MAX;
          *nonbasicflag=UNBOUNDED;
          break;
       // These cases shouldn't happen
       default: // Including ATUPPER and ATBOTH
          assert((lp->nrange[jp]!=ATUPPER) && (lp->nrange[jp]!=ATBOTH));
          printf("Incorrect call to row_choice\n");
          exit(501);
       }
    }
    else if (    !up_flag
              && ((   lp->ntype[jp]==LBOUNDED)
                   ||(lp->ntype[jp]==BOTHBOUNDS)))
// Trying to decrease non-basic variable with lower bound
    {
       switch(lp->nrange[jp])
       {
       case ABOVE:
          if((t=(lp->x[jp]-lp->nu[jp])) <= delta)
          {
            pivot=DBL_MAX;
            tempmin = t;
            *nonbasicflag=NONBASICTOFEASIBLEU;
          }
          break;
       case ATUPPER:
       // Fall through
       case BETWEEN:
          if((t=(lp->x[jp]-lp->nl[jp])) <= delta)
          {
            pivot=DBL_MAX;
            tempmin = t;
            *nonbasicflag=NOPIVOTLBOUND;
          }
          break;
       case BELOW:
          pivot = DBL_MAX;
          *nonbasicflag = UNBOUNDED;
          break;
       // These cases shouldn't happen
       default: // Including ATLOWER and ATBOTH
          assert((lp->nrange[jp]!=ATLOWER) && (lp->nrange[jp]!=ATBOTH));
          printf("Incorrect call to row_choice\n");
          exit(502);
       }
    }
// WHAT LIMITS DO THE BASIC VARIABLES IMPOSE?
// There are three cases:
//    Basic value is increasing
//    Basic value is decreasing
//    Matrix element too small to pivot on.
        ip = 0;
        for (i=1;i<=lp->m;i++)  // Go through rows
        {
// INCREASING BASIC VALUE
        // Note: Basic variables should not be ATLOWER or ATUPPER
            if(
                  (!up_flag && lp->a[jp][i]< - epsp)
               ||  ( up_flag && lp->a[jp][i]>   epsp)
              )
              switch(lp->brange[i])
              {
              case BELOW:
                 if((t=((lp->bl[i]-lp->b[i])/fabs(lp->a[jp][i]))) <= delta)
                 {
                   if(fabs(lp->a[jp][i]) > pivot)
                   {
                     pivot = fabs(lp->a[jp][i]);
                     tempmin = t;
                     ip = i;
                     *nonbasicflag=BASICTOFEASIBLEL;
                   }
                 }
                 break;
              case BETWEEN:
                 if(lp->btype[i] == LBOUNDED||lp->btype[i] == FREE)
                    break;
                 else
                 {
                    if((t=((lp->bu[i]-lp->b[i])/fabs(lp->a[jp][i]))) <= delta)
                    {
                      if(fabs(lp->a[jp][i])>pivot)
                      {
                        pivot = fabs(lp->a[jp][i]);
                        tempmin = t;
                        ip = i;  
                        *nonbasicflag=PIVOTUBOUND;
                      }
                    }
                    break;
                 }
              case ATBOTH:
                if((t=((lp->bu[i]-lp->b[i])/fabs(lp->a[jp][i]))) <= delta)
                {
                  if(fabs(lp->a[jp][i]) > pivot)
                  {
                    pivot = fabs(lp->a[jp][i]);
                    tempmin = t;
                    ip = i;
                    *nonbasicflag=PIVOTUBOUND;
                  }
                }
                break;
              case ABOVE:
                break;
              default:
                printf("\nImproper result code");
                exit (502);
              }

// DECREASING BASIC VALUE
        // Note: Basic variables should not be ATLOWER or ATUPPER
           else if(
                      ( up_flag && lp->a[jp][i] < -epsp)
                   ||  (!up_flag && lp->a[jp][i] >  epsp)
                  )
           {
              switch(lp->brange[i])
              {
              case ABOVE:
                 if((t=((lp->b[i] - lp->bu[i])/fabs(lp->a[jp][i]))) <= delta)
                 {
                   if(fabs(lp->a[jp][i]) > pivot)
                   {
                     pivot = fabs(lp->a[jp][i]);
                     tempmin = t;
                     ip = i;
                     *nonbasicflag = BASICTOFEASIBLEU;
                   }
                 }
                 break;
              case BETWEEN:
                 if(lp->btype[i] == BOTHBOUNDS||lp->btype[i] == LBOUNDED)
                 {
                    if((t=((lp->b[i] - lp->bl[i])/fabs(lp->a[jp][i]))) <= delta)
                    {
                      if(fabs(lp->a[jp][i]) > pivot)
                      {
                        pivot = fabs(lp->a[jp][i]);
                        tempmin = t;
                        ip = i;
                        *nonbasicflag = PIVOTLBOUND;
                      }
                    }
                 }         
                 break;
              case ATBOTH:
                if((t=((lp->b[i] - lp->bl[i])/fabs(lp->a[jp][i]))) <= delta)
                {                
                  if(fabs(lp->a[jp][i]) > pivot)
                  {
                    pivot = fabs(lp->a[jp][i]);
                    tempmin = t;
                    ip = i;
                    *nonbasicflag=PIVOTUBOUND;
                  }
                }
                break;
              case BELOW:
                 break;
              default:
                 printf("\nImproper result code");
                 exit (303);

              }
          } 
        }
// Move non-basic to bound
        if(pivot == DBL_MAX)
        {
          alphamin = tau;
          delta = (tempmin > alphamin)?tempmin:alphamin;         
          *dx = (up_flag)?delta:-delta;
          return 0;
        }
// Unbounded solution
        if(*nonbasicflag == UNBOUNDED)return -1;
// Pivot
        alphamin = tau/pivot;
        delta = (tempmin > alphamin)?tempmin:alphamin;
        *dx = (up_flag)?delta:-delta;
        return ip;
}

/************************************************************************
* 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,           *
*                                                                      *
************************************************************************/
