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


// 3/26/05 Changed so that readmps could be integrated 
//         into retroLP
// 3/9/05  Corrected bug; wasn't counting pivot time in statistics.

#define DI 1 // Debugging printout for refresh

//<!> Put in checks for degenerate cases

#include<math.h>
#include<stdio.h>
#include <iostream>
#include<stdlib.h>
#include<malloc.h>
#include<time.h>
#include"simplex.h"
#include"params.h"
#include "long_clock.h"

// External subroutines used
void reset(LP_state* );
void pivot(LP_state *s, long, long );
double **allocate_a(long , long);
short *allocate_svector(long);
long *allocate_lvector(long);
double *allocate_dvector(long);
void tswitch(var_type *, var_type *);
void lswitch(long *, long *);
void dswitch(double *, double *);
int to_state(FILE *, short, LP_state **);
short check_ph1(LP_state * );
void report(long, LP_state *);


// This routine eventually will be used to "refresh" the tableau
// after a given number of iterations; to recover from "check points"
// and to "crash" a basis and/or a solution vector. We also
// "reset" in the sense of Gill et al.

long refresh
(
   FILE * lpb,   // Handle for binary file 
   LP_state* lp,    // Data structure for LP (see simplex.h)
   short* p         // Phase
)
{
   long *row_list, *col_list, *row_link, *col_link;
   long *save_basic, *save_nonbasic;
   long row_count=0, col_count=0, col_count0, iter;
   long incount=0, m, n, h, i, j, k, hh, save_h, save_k;
   double *save_x, *save_b;
   double size, temp_size;
   int result;

// Save information before overwriting lp
   m = lp->m;
   n = lp->n;
   save_b = lp->b;
   lp->b = allocate_dvector(m+1);
   save_basic = lp->ibasic;
   lp->ibasic = allocate_lvector(m+1);
   save_x = lp->x;
   lp->x = allocate_dvector(n);
   save_nonbasic = lp->jnonbasic;
   lp->jnonbasic = allocate_lvector(n);

// Read in original structure
// Second argument 0 -> don't allocate memory for arrays,
   result = to_state(lpb, 0, &lp);
   if(result)
     {
       std::cerr << "Not able to read binary file (code: " << result <<")"<< std::endl;
       exit(result);
     }
// Fill in desired values
   for(i=1; i<=m; i++)
   {
      h = save_basic[i];
// Auxiallary variables
      if(h>n) lp->b[h-n] = save_b[i];
// Model variables
      else lp->x[h] = save_b[i];
   }
   for(j=1; j<=n; j++)
   {
// Auxiallary variables
      k = save_nonbasic[j];
      if(k>n) lp->b[k-n] = save_x[j];
// Model variables
      else lp->x[k] = save_x[j];
   }


// Now create linked lists of rows to remove basic variables from,
// and columns to bring into the basis (col_list)

// The net number of basic and non-basic variables interchanged
// can be no greater than the minimum of m and n.
   long minmn;
   minmn = (m<n)?m:n;

   row_link = allocate_lvector(minmn+1);
   row_list = allocate_lvector(minmn+1);
   col_link = allocate_lvector(minmn+1);
   col_list = allocate_lvector(minmn+1);
   row_link[1]=0;
   col_link[1]=0;
// Recognize end of list by 0 in link field.
   for(i=1; i<=m; i++)
   {
// These model variables should be in the basis
      if(save_basic[i]<=n)
      {
// You want to pivot save_basic[i] into the basis
         col_count++;
         col_link[col_count] = col_count+1;
         col_link[col_count+1] = 0;
         col_list[col_count] = save_basic[i];
      }
   }
   for(j=1; j<=n; j++)
   {
      if(save_nonbasic[j]>n)
// You are allowed to pivot in rows corresponding to auxilliary
// variables which leave the basis.
      {
         row_count++;
         row_link[row_count] = row_count+1;
         row_link[row_count+1] = 0;
         row_list[row_count] = save_nonbasic[j] - n;
      }
   }
#if DI
// For reinversion row_count and column_count should be the same
// May want to change latter for other applications
   if(row_count!=col_count)
   {
     printf("Row count: %ld not equal column count: %ld \n",
       row_count, col_count);
     exit(44);
   }
   printf("REFRESH with %ld rows and %ld columns\n",row_count, col_count);
   col_count0=col_count;
#endif
// Now run through columns in list to find the largest pivot
// element and repeat until all variables on list pivoted in
//<?> The way this is implemented you look for the largest
//<?> pivot over all columns and repeat; this is quadratic
//<?> in the number of columns; an alternative is to pivot
//(?) on the largest element in the NEXT column,
//<?> which is linear in the number of columns.
//(?) Better yet is to find largest pivot element as you pivot!

// Main iteration starts here:
assert(col_count0<=minmn);
for(iter=1; iter<=col_count0; iter++)
{
   k=1;
   size = 0.0;
   save_k = 0;
// Go through columns
   while(col_link[k]!=0)
   {
      j = col_list[k];
// Run through elgible rows in column j for largest pivot
      temp_size = epsi;
      hh = 0;
      h = 1;
// Look through the rows of column j
      while(row_link[h]!=0)
      {
          i = row_list[h];
          if(fabs(lp->a[j][i]) > temp_size)
          {
             temp_size = fabs(lp->a[j][i]);
             hh = h;
          }
          h = row_link[h];
      }

// If hh = 0, column j cannot enter basis without removing another
// good basis variable. Close up list and continue
      if(hh == 0)
      {
#if DI
          printf(" Desired variable could not be brought into basis in refresh\n");
          printf(" because column %ld; (variable %ld) is all 0's in refresh rows\n",
           j, lp->jnonbasic[j]);
          printf(" Index: %ld; There have been %d such columns\n\n",
           j, col_count0-incount-col_count+1);
#endif
          col_count--;
          col_list[k] = col_list[col_link[k]];
          col_link[k] = col_link[col_link[k]];
          continue;
// Continue returns to while link not 0 with same
// k as before which now points to new, closed in value
      }
      if(temp_size > size)
      {
         size = temp_size;
         save_k = k;
         save_h = hh;
      }
      k =  col_link[k];
   }
// Return if no possible columns
// Else check to see if pivot large enough
   if(save_k == 0)
   {
#if DI
          printf(" Variables which were supposed to be made basic could not be (in refresh)\n");
          printf(" There are %ld such variables\n", col_count0-incount);
          printf(" incount: %ld; col_count0: %ld; col_count: %ld; row_count: %ld\n",
            incount, col_count0, col_count, row_count);
#endif

     break;
   }
   if(size < epsi)
   {
      printf("\nPivot too small: %12.4e \n", size);
      printf(" %ld pivots successfully made\n",
        incount);
#if DI
      printf(" incount: %ld; col_count0: %ld; col_count: %ld; row_count: %ld\n",
        incount, col_count0, col_count, row_count);
#endif

      break;
   }
   incount++;
   i = row_list[save_h];
   j = col_list[save_k];
   riter++;
   pstart = clock();
   pivot(lp, i, j );
   pstop = clock();
   pclock += pstop - pstart;
   pclock = tlong.reinitLongClock(pclock);

// Remove row i from rowlist, and j from column list
   row_count--;
   row_list[save_h] = row_list[row_link[save_h]];
   row_link[save_h] = row_link[row_link[save_h]];
   col_count--;
   col_list[save_k] = col_list[col_link[save_k]];
   col_link[save_k] = col_link[col_link[save_k]];
// Exchange row and column properties
   tswitch(&(lp->btype[i]), &(lp->ntype[j]));
   lswitch(&(lp->ibasic[i]), &(lp->jnonbasic[j]));
   dswitch(&(lp->bu[i]),&(lp->nu[j]));
   dswitch(&(lp->bl[i]),&(lp->nl[j]));
   dswitch(&(lp->b[i]),&(lp->x[j]));
   if((iter%k_iter == 0)&&l_iter)
   {
      printf("======================================\n");
      printf("======================================\n");
      printf("REFRESH: Iteration: %d\n", iter);
#if DI
      printf(" incount: %ld; col_count0: %ld; col_count: %ld; row_count: %ld\n",
        incount, col_count0, col_count, row_count);
#endif
      report(l_iter,lp);
   }
}  // endfor
// "reset" for EXPAND; sets non-basics
   reset(lp);
// Calculate correct b
   for(i=0;i<=m; i++)
   {
      lp->b[i]=lp->a[0][i];
      for(j=1; j<=n; j++)
         lp->b[i]+=lp->x[j]*lp->a[j][i];
   }
// Assign ranges and check to see if reverted to Ph1
   *p = check_ph1(lp);
   if(l_start > 0)printf("REFRESH\n");
// Release allocates
   free(save_b);
   save_b = NULL;
   free(save_basic);
   save_basic = NULL;
   free(save_x);
   save_x = NULL;
   free(save_nonbasic);
   save_nonbasic = NULL;
   free(row_link);
   row_link = NULL;
   free(row_list);
   row_list = NULL;
   free(col_link);
   col_link = NULL;
   free(col_list);
   col_list = NULL;
#if DI
   printf(" %ld pivots made\n", incount);
#endif
   return incount;
}

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