
#include <stdio.h>
#include <sys/types.h>
#include <sys/times.h>
#include <sys/param.h>
#include "cpxdefs.inc"
#include "cplex_vars.h"

double precision = 1e-5; /* arithmetic precision */

int delta;   /* for memory allocation */

int nonnegative; /* = 1 if the system [s -S] should be completed by
                        adding the non negativity of the variables
                    = 0 otherwise */

int extrays;     /* = 1 if a set of extremal rays is given
                    = 0 otherwise */
                                    
double **ray_mat; /* matrix storing the rays */
double **B_mat;   /* matrix storing the system [B b] */
double **S_mat;   /* matrix storing the dominating system [S s] */

char indic_nonneg[256];  /* for output */

void
read_matrix(type, f_name, nb_lines, nb_cols)
char type;     /* 'R', 'B' or 'S' (depending on the matrix to be read) */
char *f_name;  /* reading from the file f_name */
int *nb_lines; /* #lines of the read matrix */ 
int *nb_cols;  /* #columns of the read matrix */
{
int i, j, extreme_point_saw;
double ray_type;
char buff[256];
FILE *f;

  f = fopen(f_name, "r");

  fscanf(f,"%s", buff);

  while (strcmp(buff, "begin") != 0) {
    fscanf(f,"%s", buff);
    
#ifdef TRACE
    printf("(%s)\n", buff);
#endif
  };


  if (type == 'R') {                               /* reading the rays matrix */
    fscanf(f, "%d %d", nb_lines, nb_cols);
    printf("#rows : %6d   #columns : %6d\n", *nb_lines, *nb_cols);
    
    *nb_cols = *nb_cols-1; 
                  /* dropping the first entry of the rays */
    fgets(buff, sizeof(buff),f);
    
#ifdef TRACE
    printf("(%s)  ", buff);
#endif
    
    ray_mat = (double**) malloc (sizeof(double*)*(*nb_lines));

    for (i=0; i<*nb_lines; i++)
      ray_mat[i] = (double*) malloc (sizeof(double)*(*nb_cols));
    
    extreme_point_saw = 0;
      
    for (i=0; i<*nb_lines; i++) {
    
#ifdef TRACE
      printf("%d : ", (i+1));   
#endif 
      fscanf (f, "%lf", &ray_type);
      if (ray_type == 1.0) {     /* extreme point of the projection cone */
        extreme_point_saw = 1;
      }
      else {
        for (j=0; j<*nb_cols; j++) {
          fscanf(f, "%lf", &(ray_mat [i-extreme_point_saw] [j]));

#ifdef TRACE
          printf("%lf  ", ray_mat [i-extreme_point_saw] [j]);    
#endif 
        
        }
      }
      fgets(buff, sizeof(buff),f);
       
#ifdef TRACE
      printf("(%s)  \n", buff); 
#endif 
    }
    *nb_lines = *nb_lines - extreme_point_saw;
    close(f);
  }
  else if (type == 'B') {      /* reading the [b -B] system */
  
    fscanf(f, "%d %d", nb_lines, nb_cols);
    printf("#rows : %6d   #columns : %6d\n",*nb_lines, *nb_cols);
    fgets(buff, sizeof(buff),f);

    B_mat = (double**) malloc (sizeof(double*)*(*nb_lines));

    for (i=0; i<*nb_lines; i++)
      B_mat[i] = (double*) malloc (sizeof(double)*(*nb_cols));
        
    for (i=0; i<*nb_lines; i++) {
    
#ifdef TRACE
      printf("%d : ", (i+1));
#endif

      fscanf(f, "%lf", &(B_mat [i] [*nb_cols-1]));  /* reading the rhs */

#ifdef TRACE
      printf("%lf  ", B_mat [i] [*nb_cols-1]);
#endif
      for (j=0; j<(*nb_cols)-1; j++) {
        fscanf(f, "%lf", &(B_mat [i] [j])); 
        B_mat [i] [j] = -B_mat [i] [j];     /* stored in the form B x <= b */
        
#ifdef TRACE
        printf("%lf  ", B_mat [i] [j]);
#endif

      }

      fgets(buff, sizeof(buff),f);

#ifdef TRACE
      printf("\n");
#endif

    }
    close(f);
  }
  else if (type == 'S') {
    fscanf(f, "%d %d", nb_lines, nb_cols);
    printf("#rows : %6d   #columns : %6d\n\n", *nb_lines, *nb_cols);
    fgets(buff, sizeof(buff),f);

    S_mat = (double**) malloc (sizeof(double*)*(*nb_lines));

    for (i=0; i<*nb_lines; i++)
      S_mat[i] = (double*) malloc (sizeof(double)*(*nb_cols));
        
    for (i=0; i<*nb_lines; i++) {

#ifdef TRACE
      printf("%d : ", (i+1));
#endif

      fscanf(f, "%lf", &(S_mat [i] [*nb_cols-1]));  /* reading the rhs */

#ifdef TRACE
        printf("%lf  ", S_mat [i] [*nb_cols-1]);
#endif

      for (j=0; j<(*nb_cols)-1; j++) {
        fscanf(f, "%lf", &(S_mat [i] [j])); 
        S_mat [i] [j] = -S_mat [i] [j];         /* stored in the form S x <= s */

#ifdef TRACE
        printf("%lf  ", S_mat [i] [j]);
#endif

      }

      fgets(buff, sizeof(buff),f);

#ifdef TRACE
      printf("\n");
#endif
    }
   
    nonnegative = 0;
    sprintf (indic_nonneg, "");    
    if (fscanf(f, "%s", buff) != EOF) {
      if (strcmp(buff, "end") != 0) {
        printf("## READ ERROR : unable to locate 'end' in the file for [s S]\n");
        exit(0);
      }
      else {
        if (fscanf(f, "%s", buff) != EOF) {
          if (strcmp(buff, "nonnegative") == 0) {
            nonnegative = 1;
            sprintf(indic_nonneg, "(+ non negativity inequalities)");
          }
          else {
            printf("## READ ERROR : unknow keyword '%s' in the file for [s S]\n",
                                                                           buff);
            exit(0);
          }
        }
      }
    }
    close(f);
  }
} /* read_matrix */

void
prod(mat1, nb_lines1, nb_cols1,
     mat2, nb_lines2, nb_cols2,
     mat_res)
double **mat1;    /* matrix 1 */
int nb_lines1;    /* #lines of matrix 1 */
int nb_cols1;     /* #columns of matrix 1 */
double **mat2;    /* matrix 2 */
int nb_lines2;    /* #lines of matrix 2 */
int nb_cols2;     /* #columns of matrix 1 */
double **mat_res; /* matrix for the result of the product matrix 1 * matrix 2 */
{
 int i, j, k;
 
  if (nb_cols1 != nb_lines2) {
    printf("## ERROR: Incompatible dimensions !\n");
    exit(0);
  }
  else {
    for (i=0; i<nb_lines1; i++)
      for (j=0; j<nb_cols2; j++) {
        mat_res [i] [j] = 0.0;
        for (k=0; k<nb_cols1; k++)
        mat_res [i] [j] = mat_res [i] [j] +mat1 [i] [k] * mat2 [k] [j];
    };
  };
  
#ifdef TRACE
  printf ("prod :\n");
  for (i=0; i<nb_lines1; i++) {
    printf("%d : ", i);
    for (j=0; j<nb_cols2; j++) {
      printf("%10lf  ", mat_res [i] [j]);
    };
    printf("\n");
  };
#endif
  
} /* prod */


int
find_index (to_find)
int to_find;

{  int index;
   char current_name[256];

   sprintf(current_name, "x%d", to_find);

   index = to_find;
    
   if (strcmp(cname[index],current_name) != 0) {
     index ++;
     while ((strcmp(cname[index],current_name) != 0) && (index != (to_find))) {
       if (index < mac-1) {
         index++;
       }
       else {
         index = 0;
       }
     }
     if (index == to_find) {
       printf("## ERROR : unable to find the %dth inequality\n", to_find);
       exit(0);      
     }
   }
           
   if (index != to_find) {
     printf("## WARNING : columns permutation in the LP-solver (no = %d)\n",
                                                                  to_find);
   }
   
   return(index);
} /* find_index */

main ()
{
int i, j, index, test_index, status_coef, status_bound;
int v_index[1];
double big[1], zero[1];
char upper[1];
char current_name[256];
char *probname = {"example"};
FILE *f_res;
int non_dom;            /* #non dominated inequalities */
int S_dom;              /* = 1 if S is a dominating system, = 0 otherwise */
int nz_coeff;           /* #of non zero in the LP matrix */
int nz_coeff_var;       /* #of non zero in a column of the LP matrix */
int nb_max_col_row;     /* maximum between # columns and #rows of the LP matrix */
int nb_cara;            /* maximum #characters for any name (rows or columns) */
int nb_rays;            /* #rays = #lines in the system [R] */
int dim_rays;           /* rays dimension = #columns in the system [R] */
int nb_ineqs;           /* #ineqalities in the system [b (-B)] */
int nb_vars;            /* #columns in the system [b (-B)] */
int nb_ineqs_S;         /* #ineqalities in the system [s (-S)] */
int nb_vars_S;          /* #columns in the system [s (-S)] */
double **gen_ineqs;     /* matrix [RB Rb]
                                  [S   s] */
char f_name_rays[256];  /* file name for the rays */
char f_name_B[256];     /* file name for the system [b (-B)] */
char f_name_S[256];     /* file name for the system [s (-S)] */
char f_name_res[256];   /* file name for output */

struct tms cpu;        /* for computing cpu time */

printf("For a linear system describing a polyhedron P of the form\n\n");
  
printf("                                b - A y - B x >= 0\n\n");

printf("and the list of the extremal rays u of the cone u A = 0, u >= 0, the program\n");
printf("domcheck generates the inequalities\n\n");
printf("                               (u b) - (u B) x >= 0\n\n");

printf("and outputs an equivalent minimal system of inequalities (i.e. a minimal system\n");
printf("describing the projection of the polyhedron onto the space of the variables x).\n\n");
  
printf("If a set of inequalities\n");

printf("                                    s - S x >= 0\n\n");

printf("is known to be valid for this projection, there is the option of specifying\n");
printf("this set. The program will then report if this set of inequalities is\n");
printf("sufficient or not to describe the projection (useful in particular for\n");
printf("non full dimensional polyhedra).\n\n");

printf("This program may also be used simply to remove redundancy in a given linear\n");
printf("system of the form\n");
printf("                                b - B x >= 0.\n\n");

printf("In that case, just put 'none' for the file containing the extremal rays.\n\n\n\n\n");

  printf("\nFile name for the extremal rays of the cone u A = 0, u >= 0 :\n");
  scanf("%s",f_name_rays);

  if (strcmp(f_name_rays, "none") == 0) { 
    extrays = 0;
  }
  else { 
    extrays = 1; 
    read_matrix ('R', f_name_rays, &nb_rays, &dim_rays);
  }

  printf("\nFile name for system b - B x <= 0 :\n");
  scanf("%s",f_name_B);

  read_matrix ('B', f_name_B, &nb_ineqs, &nb_vars);

  printf("\nFile name for the dominating system s - S x >= 0 : ('none' if no system)\n");
  scanf("%s", f_name_S);

  if (strcmp(f_name_S, "none") == 0) {
    nb_ineqs_S = 0;
    nb_vars_S = nb_vars;
  }
  else {
    read_matrix ('S', f_name_S, &nb_ineqs_S, &nb_vars_S);  
    if (nb_vars != nb_vars_S) {
      printf("## ERROR : [b -B] and [s -S] of incompatible dimensions\n");
      exit(0);
    };  
  };

  if (extrays == 1) {
    gen_ineqs = (double**) malloc (sizeof(double*)*(nb_rays+nb_ineqs_S));

    for (i=0; i<nb_rays+nb_ineqs_S; i++)
      gen_ineqs [i] = (double*) malloc (sizeof(double)*(nb_vars));
    
    prod(ray_mat, nb_rays, dim_rays, B_mat, nb_ineqs, nb_vars, gen_ineqs);
        /* gen ineqs is a matrix of dimensions (nb_rays+nb_ineqs_S) * nb_vars */
      
    for (i=nb_rays; i<nb_rays+nb_ineqs_S; i++) {
      for (j=0; j<nb_vars; j++)
        gen_ineqs[i] [j] = S_mat [i-nb_rays] [j];
    }
  }
  else {                       /* no rays matrix to be read */
    nb_rays = nb_ineqs;
    gen_ineqs = (double**) malloc (sizeof(double*)*(nb_ineqs+nb_ineqs_S));
/*
    for (i=0; i<nb_ineqs+nb_ineqs_S; i++)
      gen_ineqs [i] = (double*) malloc (sizeof(double)*(nb_vars));
    
    for (i=0; i<nb_ineqs; i++) {
      for (j=0; j<nb_vars; j++) {
        gen_ineqs [i] [j] = B_mat [i] [j];
      }
    }
*/ 
    for (i=nb_ineqs; i<nb_ineqs+nb_ineqs_S; i++)
      gen_ineqs [i] = (double*) malloc (sizeof(double)*(nb_vars));
    
    for (i=0; i<nb_ineqs; i++) {
      gen_ineqs [i] = B_mat [i];
    }
   
    for (i=nb_ineqs; i<nb_ineqs+nb_ineqs_S; i++) {
      for (j=0; j<nb_vars; j++)
        gen_ineqs[i] [j] = S_mat [i-nb_ineqs] [j];
    }

  }
      
 /* Initializing variables for cplex */
  

  mac = nb_rays+nb_ineqs_S;
  mar = nb_vars-1;
  delta = mar;

  if (mac >= mar) {
    nb_max_col_row = mac;
  } else {
      nb_max_col_row = mar;
    }

  nb_cara = 1;
    
  while (nb_max_col_row > 0) {
    nb_cara++;
    nb_max_col_row = (int)(nb_max_col_row/10);
  }

  if (nb_cara > name_length) {
    printf("## Number of rows/columns too big for the defined name length\n");
    exit(0);
  }

  objsen = 1; /* minimization */

  objx = (double*) malloc (sizeof(double)*(mac+delta));
  cstore = (char*) malloc (sizeof(char)*(mac+delta)*(name_length+2));
  cname = (char**) malloc (sizeof(char*)*(mac+delta));

  for (i=0; i<mac; i++) {
    sprintf(&cstore[(name_length+2)*i], "x%d", i);
    cname[i] = &cstore[(name_length+2)*i];
  }
  
/*
  rstore = (char*) malloc (sizeof(char)*(mac+delta)*(name_length+2));
  rname = (char**) malloc (sizeof(char*)*(mac+delta));

  for (i=0; i<mac; i++) {
    strcpy(&rstore[(name_length+2)*i], "C%d", i);
    rname[i] = &rstore[(name_length+2)*i];
  }
*/

  rhsx = (double*) malloc (sizeof(double)*(mar+delta));

  nz_coeff = 0;
  for (i=0; i<mac; i++) {
    for (j=0; j<mar; j++) {
      if (gen_ineqs[i] [j] != 0.0) {
        nz_coeff++;
      }
    }
  }
        
  matval = (double*) malloc (sizeof(double)*(nz_coeff+delta));
  bdl = (double*) malloc (sizeof(double)*(mac+delta));
  bdu = (double*) malloc (sizeof(double)*(mac+delta));
  x = (double*) malloc (sizeof(double)*(mac+delta));
  obj = (double*) malloc (sizeof(double));

  matbeg = (int*) malloc (sizeof(int)*(mac+delta));
  matcnt = (int*) malloc (sizeof(int)*(mac+delta));
  matind = (int*) malloc (sizeof(int)*(nz_coeff+delta));
  status_sol = (int*) malloc (sizeof(int));
  
  senx = (char*) malloc (sizeof(char)*(mar+delta));
  rhsx = (double*) malloc (sizeof(double)*(mar+delta));

  if (nonnegative) {
    for (i=0; i<mar; i++) {
      senx[i] = 'G';
    }
  }
  else {
    for (i=0; i<mar; i++) {
      senx[i] = 'E';
    }
  }

  nz_coeff = 0;

  for (i=0; i<mac; i++) {
    nz_coeff_var = 0;
    matbeg [i] = nz_coeff;
    bdl [i] = 0.0;
    bdu [i] = INFBOUND;
    for (j=0; j<mar; j++) {
      if ((gen_ineqs [i] [j] > precision) || (gen_ineqs [i] [j] < -precision)) {
        matval [nz_coeff] = gen_ineqs [i] [j];
        matind[nz_coeff] = j;
        nz_coeff++;
        nz_coeff_var++;
      }
    }
    matcnt [i] = nz_coeff_var;
  };
    
  matbeg [mac] = nz_coeff;


  for (i=0; i<mac; i++) {
    objx [i] = gen_ineqs [i] [nb_vars-1];
  }

  non_dom = 0;
  zero[0] = 0.0;
  big[0] = INFBOUND;
  upper[0] = 'U';
  status_bound = 1;

  times(&cpu);
  printf("\n\nInitialisation completed       cpu: %8.2lf [sec]\n",
              (double)cpu.tms_utime/(double)HZ);
  printf("\n\nStarting the domination check\n\n");
  

  for (test_index=0; test_index<mac; test_index++) {
  
    index = find_index(test_index);
  
      /* test if the (test_index)th generated inequation is dominated by the others */
     
    if (test_index == 0) {
          /* first tested inequation */
          
      bdu [index] = 0.0;
      for (i=0; i<mar; i++) {
        rhsx [i] = gen_ineqs [index] [i];
      }
      
#ifdef TRACEALL
      for (i = 0; i<mac; i++) {
        printf("objx[%d] = %lf  matbeg[%d] = %d  matcnt[%d] = %d  bdl[%d] = %lf  bdu[%d] = %lf\n",
                i, objx[i], i, matbeg[i], i, matcnt[i], i, bdl[i], i, bdu[i]);};
  
      for (i = 0; i<mar; i++) { 
        printf("rhsx [%d] = %lf  senx [%d] = %c\n", i, rhsx[i], i, senx[i]);};
  
      for (i = 0; i<nz_coeff; i++) { 
        printf("matind [%d] = %d  matval [%d] = %lf\n", i, matind[i], i, matval[i]);};
#endif

      macsz = mac+delta;
      marsz = mar+delta;
      matsz = nz_coeff+delta;
 
      lp =  loadprob (probname, mac, mar, mae, objsen, objx, rhsx,
                      senx, matbeg, matcnt, matind, matval, bdl,
                      bdu, rngval, nrowind, etype, enzbeg, enzcnt,
                      enzind, enzval, dataname, objname, rhsname, rngname,
                      bndname, cname, cstore, rname, rstore, ename,
                      estore, macsz, marsz, matsz, maesz, enzsz, cstorsz, rstorsz,
                      estorsz);

      if (lp == NULL) { 
        printf("## CPLEX-ERROR: can't create the LP-variable \n");
        exit(0);
      }
          
    }
    else { 
       /* test_index > 0 */
      bdu [index] = 0.0;
      v_index[0] = index; 
      status_bound = chgbds(lp, 1, v_index, upper, zero);
    
      if (status_bound == 1) {
        printf("##  CPLEX-ERROR: unable to change bound to zero for variable %d\n", index);
        exit(0);
      }

      for (i = 0; i<mar; i++) {
        status_coef = chgcoef(lp, i, -1, gen_ineqs [index] [i]);
        
        if (status_coef == 0) {
#ifdef TRACEALL
          printf("coeff %d change en %lf\n", i, gen_ineqs [index] [i]);
#endif
        }
        else {
          printf("##  CPLEX-ERROR: unable to change rhs for row %d\n", i);
          exit(0);
        }
      }    
    }

    status = optimize (lp);
  
    status = solution (lp, status_sol, obj, x, pi, slack, dj);

    if (((*status_sol == 1)&&(*obj > gen_ineqs [index] [nb_vars-1])) ||
         (*status_sol == 2)) {
         
      bdu [index] = INFBOUND;   
      status_bound = chgbds(lp, 1, v_index, upper, big);
    
      if (status_bound == 1) {
        printf("##  CPLEX-ERROR: error in changing bound to INFBOUND for variable %d\n", index);
        exit(0);
      }
      non_dom++;
    }
        /* the inequation is not dominated and kept in the pool */

    if (*status_sol > 3) {
      printf("LP not solved to optimality for index = %d\n", index);
    }
    /*    
    unloadprob(&lp);
    */ 
    if ((((test_index+1) % 50) == 0) || (test_index == mac-1)) {
      times(&cpu);
      printf("#tested inequalities   : %6d    #non dominated : %d    cpu : %8.2lf [sec]\n",
                   (test_index+1), non_dom, (double)cpu.tms_utime/(double)HZ);
    }

    if ((test_index == nb_rays-1)) {
      if (non_dom == 0) {
        S_dom = 1;
      }
      else {
        S_dom = 0;
      }
    }
  
  };

  printf("\n#inequalities in [S s] : %6d   #non dominated inequalities : %6d\n",
                                                           nb_ineqs_S, non_dom);
  if ((nb_ineqs_S + nonnegative) > 0) {
    if (S_dom) {
        printf("\n\n[S s] %s is a dominating system\n\n", indic_nonneg);
    }
    else {
      printf("\n\n[S s] %s is not a dominating system\n\n", indic_nonneg);
    }
  }

  printf("\nOutput file:\n");
  
  scanf("%s", f_name_res);
  
  f_res = fopen(f_name_res, "w");

  fprintf(f_res, "begin\n%10d  %10d          real\n", non_dom, nb_vars);

  for (index=0; index<mac; index++) {
    i = find_index(index);
    if (bdu [i] > 0.0) {
      if ((gen_ineqs [i] [nb_vars-1]-((int)gen_ineqs [i] [nb_vars-1]) > precision) ||
          (gen_ineqs [i] [nb_vars-1]-((int)gen_ineqs [i] [nb_vars-1]) < -precision)) {
        fprintf(f_res, "%9lf    ", gen_ineqs [i] [nb_vars-1]);
      }
      else {
        fprintf(f_res, "%3d    ", ((int) gen_ineqs [i] [nb_vars-1]));        
      }
      for (j=0; j<nb_vars-1; j++)
      if ((gen_ineqs [i] [j]-((int)gen_ineqs [i] [j]) > precision) ||
          (gen_ineqs [i] [j]-((int)gen_ineqs [i] [j]) < -precision)) {
        fprintf(f_res, "%9lf ", -gen_ineqs [i] [j]);
      }
      else {
        fprintf(f_res, "%3d ", (-(int)gen_ineqs [i] [j]));
      }
      if (i<nb_rays) {
        fprintf(f_res, "     (%d)\n", (i+1));
      } else {
          fprintf(f_res, "     (S%d)\n", (i-nb_rays+1));
        }
          
    }
  }
  fprintf(f_res, "end\n");
  if (nonnegative) {
    fprintf(f_res, "nonnegative\n");
  } 
  close (f_res);
  
  times(&cpu);
  printf("total cpu time: %8.2lf [sec]\n", (double)cpu.tms_utime/(double)HZ);
 
}