/*
Bohdan L. Kaluzny: "fourier.c"
compile: gcc -o fourier fourier.c
run: ./fourier <inputfile> <outputfile> <redundancy executable>

for example: "./fourier project1.ine project1res.ext ../Research/lrslib-040/redund" 

input file: H-Rep of convex polyhedron (see www.cs.mcgill.ca/~avis/lrs.htm) 
            with "project N n_1 n_2 n_3 ... n_n" after "end"
where N is the dimension to project onto.  n_i's are the variables to keep!
output file: Program outputs H-Rep after projection with no redundant constraints.

the heart: fourier elimination with redundancy removal after each iteration!
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

FILE *fp_input;
FILE *fp_output;
int input_A[100000][100];
int matrix_T[50000][1000];
int output_A[50000][100];
int Col[100];
char name[100];
int iterations;
int num_redund;
int m_output;

/**************************************************************/
void compute_T(int * mT, int *nT, int mA, int proj)
{
  int i,j,k,l,count,a1,a2;
  count = 1;
  //initialize T matrix
  for(i=1;i<=*mT;i++)
    {
      for(j=1;j<=*nT;j++)
	{
	  matrix_T[i][j]=0;
	}
    }
  for(i=1;i<=mA;i++)
    {
      if(input_A[i][proj] == 0)
	{ //we create e row to T
	  for(j=1;j<=mA;j++)
	    {
	      if(j==i)
		matrix_T[count][j]=1;
	      else
		matrix_T[count][j]=0;
	    }
	  count++;
	}
      else if(input_A[i][proj] > 0)
	{ //we create a row in T from I- x I+
	  a1 = input_A[i][proj];
	  for(k=i;k<=mA;k++)
	    { //look for a row with a < 0
	      if(input_A[k][proj] < 0)
		{ //now create a row in T
		  a2 = (-1)*input_A[k][proj];
		  for(j=1;j<=mA;j++)
		    {
		      if(j==i)
			matrix_T[count][j] = a2;
		      else if(j==k)
			matrix_T[count][j] = a1;
		      else
			matrix_T[count][j] = 0;		       
		    }
		  count++;
		}
	    }
	}
      else if(input_A[i][proj] < 0)
	{ //we create a row in T from I- x I+
	  a1 = (-1)*input_A[i][proj];
	  for(k=i;k<=mA;k++)
	    { //look for a row with a > 0
	      if(input_A[k][proj] > 0)
		{ //now create a row in T
		  a2 = input_A[k][proj];
		  for(j=1;j<=mA;j++)
		    {
		      if(j==i)
			matrix_T[count][j] = a2;
		      else if(j==k)
			matrix_T[count][j] = a1;
		      else
			matrix_T[count][j] = 0;		       
		    }
		  count++;
		}
	    }
	}
    }
  *mT = count-1;
  *nT = mA;
}
/**************************************************************/
void matrix_multiply(int m1,int m2,int n2, int *oldm, int *oldn)
{
  int i,j,k;
  
  for(i=0;i<=*oldm;i++)
    for(j=0;j<=*oldn;j++)
      output_A[i][j]=0;

  for(i=1;i<=m1;i++)
    {
      for(j=0;j<n2;j++)
	{
	  for(k=1;k<=m2;k++)
	    {
	      output_A[i][j] += matrix_T[i][k]*input_A[k][j];
	    }
	}
    }
  *oldm=m1;
  *oldn=n2;
}
/**************************************************************/
void print_output(int m,int n)
{
  int i,j;
  
  fp_output = fopen("lrs_fou.ine","w");    //Output file
  fprintf(fp_output,"H-Representation\nbegin\n%d %d integer\n",m,n);
  for(i=1;i<=m;i++)
    {
      for(j=0;j<n;j++)
	{
	  fprintf(fp_output,"%d ", output_A[i][j]);
	}
      fprintf(fp_output,"\n");
    }
  fprintf(fp_output,"end");
  m_output = m;
}
/**************************************************************/
void print_last_output(int m,int n)
{
  int i,j;
  
  fp_output = fopen("lrs_fou.ine","w");    //Output file
  fprintf(fp_output,"H-Representation\nbegin\n%d %d integer\n",m,n);
  for(i=1;i<=m;i++)
    {
      for(j=0;j<n;j++)
	{
	  fprintf(fp_output,"%d ", input_A[i][j]);
	}
      fprintf(fp_output,"\n");
    }
  fprintf(fp_output,"end");
}
/**************************************************************/
void read_input(int *m, int *n)
{
  int i,j,temp;

  fscanf(fp_input, "%s", name);
  while(strcmp(name,"begin")!=0)
    {
     fscanf(fp_input, "%s", name);  //any garbage before "begin"
     //num_redund++;
    }
  //num_redund = num_redund-3;
  //if(num_redund < 0)
  //   num_redund = 0;
  fscanf(fp_input,"%d %d %s", m, n, name);        
  for(i=1;i<=*m;i++)
    { 
      for(j=0;j<*n;j++)
	{
          fscanf(fp_input,"%d",&temp);
	  input_A[i][j] = temp;
	}
    }
  
  fscanf(fp_input,"%s", name); //Read "end"
  if(m_output > 0)
    {  
      num_redund += m_output-*m;
    }
}
/**************************************************************/
int main(int argc, char *argv[])
{
  int m,n,mT,nT,i,j,k,proj,oldm,oldn;
  char system_call[100];
  
  num_redund = 0;
  m_output = 0;
  mT = 1;
  nT = 1;
  oldm = 1;
  oldn = 1;
  fp_input = fopen(argv[1],"r");     //Input file
  read_input(&m,&n);
  fscanf(fp_input,"%s %d", name, &iterations);
  for(i=1;i<n;i++)
    {
       Col[i] = 0;
    }
  for(i=1;i<=iterations;i++)
    {
       fscanf(fp_input, "%d", &proj);
       Col[proj] = 1;
    }
  fclose(fp_input);

  for(i=1;i<=(n-iterations-1);i++)
    {
      //find next variable to remove
      proj = 0;
      j = 1;
      while(proj == 0 && j < n)
        {
	  if(Col[j] == 0)
 	    {
              Col[j] = 1;
              proj = j;
            }
          j++;
	}
      if(proj != 0)
        {
      	  printf("\nComputing projection matrix...\n");
      	  compute_T(&mT,&nT,m,proj);
          printf("Matrix multiplication...\n");
          matrix_multiply(mT,m,n, &oldm, &oldn);
          printf("Printing output to file...\n");
	  printf("External call for redundancy removal...\n");
          print_output(mT,n);
          fclose(fp_output);
	  sprintf(system_call,"rm -f lrs_fou.ext");
          system(system_call);
	  sprintf(system_call,"%s lrs_fou.ine lrs_fou.ext", argv[3]);
	  system(system_call);
      
          if (i < (n-iterations-1)) 
           {
            //re-initialization
	    sprintf(system_call,"rm -f lrs_fou.ine");
	    system(system_call);
	    printf("Reading file...\n");
	    fp_input = fopen("lrs_fou.ext","r");     //Input file
	    read_input(&m,&n);
	    fclose(fp_input);
           }
	  else
	    {
	      printf("Reading file...\n");
	      fp_input = fopen("lrs_fou.ext","r");     //Input file
	      read_input(&m,&n);
	      fclose(fp_input);
	      sprintf(system_call,"rm -f lrs_fou.ext");
	      system(system_call);
	      print_last_output(m,n);
	    }
        }
   }
  sprintf(system_call,"mv lrs_fou.ine %s", argv[2]);
  system(system_call);
  printf("\nProject Computed:  \n\t%d redundant constraints removed in process\n\tplease see %s for output\n\n",num_redund,argv[2]);
}
/**************************************************************/