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


// This file contains general utilities.
// 4/11/05 Added a C++ flavor file opener get_file
// 3/26/05 Changed so that readmps could be integrated 
//         into retroLP
#include <stdio.h>
#include <math.h>
#include "simplex.h"
#include "params.h"
#include <malloc.h>
#include <stdlib.h>
#include <string>
#include <fstream>
#include <iostream>

void tableau_out(LP_state *);
LP_state *allocate_LP_state();
double **allocate_a(long , long);
short *allocate_svector(long);
int *allocate_ivector(long);
long *allocate_lvector(long);
double *allocate_dvector(long);
void free_a(double ** , long , long );
void sswitch(short *, short *);
void lswitch(long *, long *);
void dswitch(double *, double *);
char *out_type(var_type);
char *out_range(var_range);
char *out_result(iter_result);
void out_istate(long,LP_state *);
void out_jstate(long,LP_state *);
void init_brange(LP_state *);
int  get_file(std::ifstream&, const std::string&, int);
void empty_file(const char * , FILE *);

//===========================TABLEAUX_OUT==============
void tableau_out(LP_state *s)
{
   long i, j;
   printf("\nVariable indices     ");
   for(j=1;j<=s->n;j++)
   {
      printf(" %5d ", s->jnonbasic[j]);
      if(!(j%8))printf("\n");
   }
   printf("\n");

   // Phase 2 objective
   printf("Ph2 Obj: %10.2e |\n           ", -maxfactor*s->b[0] );
   for(j=0;j<=s->n;j++)
   {
      printf("%10.2e",s->a[j][0]);
      if(!(j%8))printf("\n");
   }
   printf("\n");
   printf("Basic indices/values\n");
   for(i=1;i<=s->m;i++)
   {
      printf("%3d %10.2e |\n           ", s->ibasic[i],s->b[i] );
      for(j=0;j<=s->n;j++)
      {
         printf("%10.2e",s->a[j][i]);
         if(!(j%8))printf("\n");
      }
      printf("\n");
   }
   // Phase 1 objective
   printf("Ph1 Obj: %10.2e |\n           ", s->b[s->m + 1] );
   for(j=0;j<=s->n;j++)
   {
      printf("%10.2e",s->a[j][s->m + 1]);
      if(!(j%8))printf("\n");
   }
   printf("\n");
   printf("\nNon-basics values\n");
   for(j=1;j<=s->n;j++)printf("%10.2e", s->x[j]);
   printf("\n\n");
}
//=======================ALLOCATE LP_STATE====================

LP_state *allocate_LP_state()
{
   size_t size;
   size=sizeof(LP_state);
   LP_state *state;
   /* allocate memory for LP_state structure, state */
   if ((state = (LP_state *)malloc(size)) == NULL)
   {
       printf("Not enough memory to allocate LP_state\n");
       exit(1);  /* terminate program if out of memory */
   }
    return state;
 }
//<?> This could be a good spot for a template!?
//==========================ALLOCATE_A========================

double **allocate_a(long m, long n)
{   int i, j;
    double **a;
    a=(double **) malloc((n+1)*sizeof(double *)); // allocate column
    // pointers
    if(a==NULL) {printf("malloc failed allocating matrix a"); exit(10);}
    for(j=0;j<=n;j++)
    { a[j]=(double *) malloc((m+2)*sizeof(double));
      if(a[j]==NULL) {printf("malloc failed allocating matrix a"); exit(11);}
      for(i=0;i<=m+1;i++) a[j][i]=0.0;
    }
    return a;
}
//===========================ALLOCATE_LVECTOR==========================
long *allocate_lvector(long nn)
/* allocate a long vector with subscript range v[0..nn] */
{
    long *v;
    v=(long *)malloc((size_t) ((nn+1)*sizeof(long)));
    if (!v) {printf("malloc failed allocating long vector"); exit(11);}
    v[0]=0; // not always used
    return v;
}

//===========================ALLOCATE_DVECTOR==========================
double *allocate_dvector(long nn)
/* allocate a double vector with subscript range v[0..nn] */
{
    double *v;
    v=(double *)malloc((size_t)((nn+1)*sizeof(double)));
    if (!v) {printf("malloc failed allocating double vector"); exit(12);}
    v[0]=0.0; // not always used
    return v;
}

//===========================ALLOCATE_IVECTOR==========================
int *allocate_ivector(long nn)
/* allocate a int vector with subscript range v[0..nn] */
{
    int *v;
    v=(int *)malloc((size_t) ((nn+1)*sizeof(int)));
    if (!v) {printf("malloc failed allocating int vector"); exit(13);}
    v[0]=0; // not always used
    return v;
}

//===========================ALLOCATE_SVECTOR==========================
short *allocate_svector(long nn)
/* allocate a short vector with subscript range v[0..nn] */
{
    short *v;
    v=(short *)malloc((size_t) ((nn+1)*sizeof(short)));
    if (!v) {printf("\nmalloc failed allocating short vector\n"); exit(14);}
    v[0]=0; // not always used
    return v;
}
//===========================ALLOCATE_EVECTOR==========================
// Allocation for enumerated types
void *allocate_evector(long nn, size_t size)
/* allocate vector of enumerated type with subscript range v[0..nn] */
{
    void *v;
    v=(void *)malloc((size_t) ((nn+1)*size));
    if (!v) {printf("\nmalloc failed allocating enumerated vector\n");
       exit(14);}
    return v;
}

//===========================ALLOCATE_CVECTOR==========================
char *allocate_cvector(long nn)
/* allocate a char buffer of length nn + 1
   allows for string terminator */
{
    char *v;
    v=(char *)malloc((size_t) ((nn+1)*sizeof(char)));
    if (!v) {printf("\nmalloc failed allocating char buffer\n"); exit(14);}
    v[0]=' '; // in case v[0] not be used
    return v;
}
//==========================ALLOCATE_STRVECTOR=========================
char *allocate_strvector(int nn,long n)
/* allocate enough space for a avector of n+1 (char) buffers each
   of length nn + 1 allows for string terminator, and for n+1
   entries
*/
{
    char *v;
    v=(char *)malloc((size_t) ((n+1)*(nn+1)*sizeof(char)));
    if (!v) {printf("\nmalloc failed allocating string vector buffer\n"); exit(15);}
    return v;
}

//===========================FREE_A====================================
void free_a(double ** a, long m, long n)
{
    long j;

    for(j=0;j<n+1;j++)
    {
       free((void *) a[j]);
       a[j] = NULL;
    }
    free((void *) a);
}
//<!?> Use template?!
//=============================== TSWITCH =================================
void tswitch(var_type *x, var_type *y)
{
   var_type temp;
   temp = *x;
   *x = *y;
   *y = temp;
   return;
}
//=============================== LSWITCH =================================
void lswitch(long *x, long *y)
{
   long temp;
   temp = *x;
   *x = *y;
   *y = temp;
   return;
}
//=============================== DSWITCH =================================
void dswitch(double *x, double *y)
{
   double temp;
   temp = *x;
   *x = *y;
   *y = temp;
   return;
}
//================================ OUT_TYPE ==========================
char *out_type(var_type type)
{
   char *tp;
   switch(type){case FREE: tp = "   FREE   \0";break;
         case LBOUNDED: tp = " LBOUNDED \0";break; case UBOUNDED:
         tp = " UBOUNDED \0";break; case BOTHBOUNDS: tp = "BOTHBOUNDS\0";
         break; case FIXED: tp = "  FIXED   \0";break;default: tp = "?????\0";}
         return tp;
}
//=============================== OUT_RANGE ========================
char *out_range(var_range range)
{
  char *rng;
  switch(range){case BELOW: rng = "  BELOW   \0";break;
         case ATLOWER: rng = " ATLOWER  \0";break; case BETWEEN:
         rng = " BETWEEN  \0";break; case ATUPPER: rng = " ATUPPER  \0";
         break; case ABOVE: rng = "  ABOVE   \0";break;
         case ATBOTH: rng = "  ATBOTH  \0";break; default: rng = "?????\0";}
         return rng;
}
//============================== OUT_RESULT =======================
char *out_result(iter_result result)
{
   char *rslt;
   switch(result)
   {
     case UNBOUNDED:             rslt = "    UNBOUNDED    \0";break;
     case PIVOTLBOUND:           rslt = "PIVOTLBOUND  \0";break;
     case PIVOTUBOUND:           rslt = "PIVOTUBOUND  \0";break;
     case NOPIVOTLBOUND:         rslt = "NOPIVOTLBOUND\0";break;
     case NOPIVOTUBOUND:         rslt = "NOPIVOTUBOUND\0";break;
     case BASICTOFEASIBLEL:      rslt = "BASICTOFEASIBLEL\0";break;
     case BASICTOFEASIBLEU:      rslt = "BASICTOFEASIBLEU\0";break;
     case NONBASICTOFEASIBLEL:   rslt = "NONBASICTOFEASIBLEL\0";break;
     case NONBASICTOFEASIBLEU:   rslt = "NONBASICTOFEASIBLEU\0";break;
     case NONDEGENERATEPIVOT:    rslt = "NONDEGENERATEPIVOT\0";break;
     default: rslt = "?????\0";}
     return rslt;
};
//============================= OUT_ISTATE =====================
void out_istate(long i,LP_state *lp)
{
   printf("State of basic variable %ld\n", i);
   if(i!=0)
   {
      printf("i, index, type, upper, lower, range, value:\n %ld, %ld,\
      %s, %12.2e, %12.2e, %s, %12.4e\n", i, lp->ibasic[i], out_type(lp->btype[i]),
      lp->bu[i], lp->bl[i], out_range(lp->brange[i]), lp->b[i]);
   }
   else // Variable going to bound; no pivot row
   {
      printf("To bound; no pivot row\n");
   }

}
//============================= OUT_JSTATE =====================
void out_jstate(long j,LP_state *lp)
{
   printf("State of non-basic variable %ld\n", j);
   printf("j, index, type, upper, lower, range, value:\n %ld, %ld,\
   %s, %12.2e, %12.2e, %s, %12.4e\n", j, lp->jnonbasic[j], out_type(lp->ntype[j]),
   lp->nu[j], lp->nl[j], out_range(lp->nrange[j]), lp->x[j]);
}
//======================== INIT_BRANGE ==============
// CALCULATE INITIAL STATE FOR BASIC VARIABLES
void init_brange(LP_state *lp)
{
   long row, col;
// CALCULATE INITIAL VALUE FOR BASIC VARIABLES
   for(row=0;row<=lp->m + 1; row++)
   {
       lp->b[row]=lp->a[0][row];
       for(col=1; col<=lp->n; col++)
          lp->b[row]+=lp->x[col]*lp->a[col][row];
   }
// Fill in brange based on b and bounds
// We don't allow BASIC variables to be AT bounds
// except ATBOTH(11/28/99)
   for (row=1; row<=lp->m; row++)
   {
      switch(lp->btype[row])
      {
         case FREE:
         lp->brange[row]= BETWEEN;
         break;
         case LBOUNDED:
         if(lp->b[row] < lp->bl[row]-epsb)
         {
            lp->brange[row]= BELOW;
         }
         else
         {
            lp->brange[row]= BETWEEN;
         }
         break;
         case UBOUNDED:
         if(lp->b[row] > lp->bu[row]+epsb)
         {
            lp->brange[row]= ABOVE;
         }
         else
         {
            lp->brange[row]= BETWEEN;
         }
         break;
         case BOTHBOUNDS:
         if(lp->b[row] < lp->bl[row]-epsb)
         {
            lp->brange[row]= BELOW;
         }
         else
         {
            if (lp->b[row] > lp->bu[row]+epsb)
            {
               lp->brange[row]= ABOVE;
            }
            else lp->brange[row]= BETWEEN;
         }
         break;
         case FIXED:
         if(lp->b[row] < lp->bl[row]-epsb)
         {
            lp->brange[row]= BELOW;
         }
         else
         {
            if((lp->b[row] >= lp->bl[row]-epsb)&&(lp->b[row]<=lp->bl[row]+epsb))
              lp->brange[row]= ATBOTH;
            else
            {
               lp->brange[row]= ABOVE;
            }
         }
         break;
         default:
         printf("\n Error in brange: row = %ld\n",row);
         exit(304);
      } //endswitch
   } // endfor
}

//======================== GET_FILE ==============

// These routines open a file and attaches it to handle
// get_file does this for C++ streams
// cget_file is legacy for C type treatment of mps files
// cbget_file is legacy for C type treatment of binary data file

#include <string>
#include <fstream>
#include <iostream>

int get_file( std::ifstream& file_handle, const std::string &file_path, int errorno)
{
  int e = 0;
  file_handle.close(); // close if already open
  file_handle.clear(); // clear existing errors, if any
  // Give it a try!
  file_handle.open(file_path.c_str());
  if(!file_handle)
  {
    std::cerr << "error: cannot open file: "
	      << file_path << std::endl;
    e = errorno ? errorno:1;
  }
  return e;
}

//======================== CGET_FILE ==============

// Opens C type text file and attaches it to handle

#include <stdio.h>

FILE * cget_file( const char * file_path, int errorno)
{
  FILE * file_handle;
  file_handle = fopen(file_path, "rt");
  if(file_handle == NULL)
  {
    printf("error: cannot open file: %s\n", file_path);
  }
  return file_handle;
}

// Opens C type binary file and attaches it to handle
// Used for bfile.lpb, binary problem file

//======================== CBGET_FILE ==============

#include <stdio.h>

FILE * cbget_file(const char * file_path, int errorno)
{
  FILE * file_handle;
  file_handle = fopen(file_path, "w+b");
  if(file_handle == NULL)
  {
    printf("error: cannot open binary output file: %s\n", file_path);
  }
  return file_handle;
}
// This is used to empty out bfile.lpb at end of run
void empty_file(const char * file_path, FILE * file_handle)
{
  freopen(file_path, "w", file_handle);
  fclose(file_handle);
}

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