/* cddproj.c:  Polyhedral Projections in cddlib
   written by Komei Fukuda, fukuda@cs.mcgill.ca
   Version 0.94, Aug. 4, 2005
*/

/* cddlib : C-library of the double description method for
   computing all vertices and extreme rays of the polyhedron 
   P= {x :  b - A x >= 0}.  
   Please read COPYING (GNU General Public Licence) and
   the manual cddlibman.tex for detail.
*/

#include "setoper.h"  /* set operation library header (Ver. June 1, 2000 or later) */
#include "cdd.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>

dd_MatrixPtr dd_BlockElimination(dd_MatrixPtr M, dd_colset delset, dd_ErrorType *error)
/* Eliminate the variables (columns) delset by
   the Block Elimination with dd_DoubleDescription algorithm.

   Given (where y is to be eliminated):
   c1 + A1 x + B1 y >= 0
   c2 + A2 x + B2 y =  0

   1. First construct the dual system:  z1^T B1 + z2^T B2 = 0, z1 >= 0.
   2. Compute the generators of the dual.
   3. Then take the linear combination of the original system with each generator.
   4. Remove redundant inequalies.

*/
{
  dd_MatrixPtr Mdual=NULL, Mproj=NULL, Gdual=NULL;
  dd_rowrange i,h,m,mproj,mdual,linsize;
  dd_colrange j,k,d,dproj,ddual,delsize;
  dd_colindex delindex;
  mytype temp,prod;
  dd_PolyhedraPtr dualpoly;
  dd_ErrorType err=dd_NoError;
  dd_boolean localdebug=dd_FALSE;

  *error=dd_NoError;
  m= M->rowsize;
  d= M->colsize;
  delindex=(long*)calloc(d+1,sizeof(long));
  dd_init(temp);
  dd_init(prod);

  k=0; delsize=0;
  for (j=1; j<=d; j++){
    if (set_member(j, delset)){
      k++;  delsize++;
      delindex[k]=j;  /* stores the kth deletion column index */
    }
  }
  if (localdebug) dd_WriteMatrix(stdout, M);

  linsize=set_card(M->linset);
  ddual=m+1;
  mdual=delsize + m - linsize;  /* #equalitions + dimension of z1 */

  /* setup the dual matrix */
  Mdual=dd_CreateMatrix(mdual, ddual);
  Mdual->representation=dd_Inequality;
  for (i = 1; i <= delsize; i++){
    set_addelem(Mdual->linset,i);  /* equality */
    for (j = 1; j <= m; j++) {
      dd_set(Mdual->matrix[i-1][j], M->matrix[j-1][delindex[i]-1]);
    }
  } 

  k=0;
  for (i = 1; i <= m; i++){
    if (!set_member(i, M->linset)){
      /* set nonnegativity for the dual variable associated with
         each non-linearity inequality. */
      k++;
      dd_set(Mdual->matrix[delsize+k-1][i], dd_one);  
    }
  } 
  
  /* 2. Compute the generators of the dual system. */
  dualpoly=dd_DDMatrix2Poly(Mdual, &err);
  Gdual=dd_CopyGenerators(dualpoly);

  /* 3. Take the linear combination of the original system with each generator.  */
  dproj=d-delsize;
  mproj=Gdual->rowsize;
  Mproj=dd_CreateMatrix(mproj, dproj);
  Mproj->representation=dd_Inequality;
  set_copy(Mproj->linset, Gdual->linset);

  for (i=1; i<=mproj; i++){
    k=0;
    for (j=1; j<=d; j++){
      if (!set_member(j, delset)){
        k++;  /* new index of the variable x_j  */
        dd_set(prod, dd_purezero);
        for (h = 1; h <= m; h++){
          dd_mul(temp,M->matrix[h-1][j-1],Gdual->matrix[i-1][h]); 
          dd_add(prod,prod,temp);
        }
        dd_set(Mproj->matrix[i-1][k-1],prod);
      }
    }
  }
  if (localdebug) printf("Size of the projection system: %ld x %ld\n", mproj, dproj);
  
  dd_FreePolyhedra(dualpoly);
  free(delindex);
  dd_clear(temp);
  dd_clear(prod);
  dd_FreeMatrix(Mdual);
  dd_FreeMatrix(Gdual);
  return Mproj;
}


dd_MatrixPtr dd_FourierElimination(dd_MatrixPtr M,dd_ErrorType *error)
/* Eliminate the last variable (column) from the given H-matrix using 
   the standard Fourier Elimination.
 */
{
  dd_MatrixPtr Mnew=NULL;
  dd_rowrange i,inew,ip,in,iz,m,mpos=0,mneg=0,mzero=0,mnew;
  dd_colrange j,d,dnew;
  dd_rowindex posrowindex, negrowindex,zerorowindex;
  mytype temp1,temp2;
  dd_boolean localdebug=dd_FALSE;

  *error=dd_NoError;
  m= M->rowsize;
  d= M->colsize;
  if (d<=1){
    *error=dd_ColIndexOutOfRange;
    if (localdebug) {
      printf("The number of column is too small: %ld for Fourier's Elimination.\n",d);
    }
    goto _L99;
  }

  if (M->representation==dd_Generator){
    *error=dd_NotAvailForV;
    if (localdebug) {
      printf("Fourier's Elimination cannot be applied to a V-polyhedron.\n");
    }
    goto _L99;
  }

  if (set_card(M->linset)>0){
    *error=dd_CannotHandleLinearity;
    if (localdebug) {
      printf("The Fourier Elimination function does not handle equality in this version.\n");
    }
    goto _L99;
  }

  /* Create temporary spaces to be removed at the end of this function */
  posrowindex=(long*)calloc(m+1,sizeof(long));
  negrowindex=(long*)calloc(m+1,sizeof(long));
  zerorowindex=(long*)calloc(m+1,sizeof(long));
  dd_init(temp1);
  dd_init(temp2);

  for (i = 1; i <= m; i++) {
    if (dd_Positive(M->matrix[i-1][d-1])){
      mpos++;
      posrowindex[mpos]=i;
    } else if (dd_Negative(M->matrix[i-1][d-1])) {
      mneg++;
      negrowindex[mneg]=i;
    } else {
      mzero++;
      zerorowindex[mzero]=i;
    }
  }  /*of i*/

  if (localdebug) {
    dd_WriteMatrix(stdout, M);
    printf("No of  (+  -  0) rows = (%ld, %ld, %ld)\n", mpos,mneg, mzero);
  }

  /* The present code generates so many redundant inequalities and thus
     is quite useless, except for very small examples
  */
  mnew=mzero+mpos*mneg;  /* the total number of rows after elimination */
  dnew=d-1;

  Mnew=dd_CreateMatrix(mnew, dnew);
  dd_CopyArow(Mnew->rowvec, M->rowvec, dnew);
/*  set_copy(Mnew->linset,M->linset);  */
  Mnew->numbtype=M->numbtype;
  Mnew->representation=M->representation;
  Mnew->objective=M->objective;


  /* Copy the inequalities independent of x_d to the top of the new matrix. */
  for (iz = 1; iz <= mzero; iz++){
    for (j = 1; j <= dnew; j++) {
      dd_set(Mnew->matrix[iz-1][j-1], M->matrix[zerorowindex[iz]-1][j-1]);
    }
  } 

  /* Create the new inequalities by combining x_d positive and negative ones. */
  inew=mzero;  /* the index of the last x_d zero inequality */
  for (ip = 1; ip <= mpos; ip++){
    for (in = 1; in <= mneg; in++){
      inew++;
      dd_neg(temp1, M->matrix[negrowindex[in]-1][d-1]);
      for (j = 1; j <= dnew; j++) {
        dd_LinearComb(temp2,M->matrix[posrowindex[ip]-1][j-1],temp1,\
          M->matrix[negrowindex[in]-1][j-1],\
          M->matrix[posrowindex[ip]-1][d-1]);
        dd_set(Mnew->matrix[inew-1][j-1],temp2);
      }
      dd_Normalize(dnew,Mnew->matrix[inew-1]);
    }
  } 


  free(posrowindex);
  free(negrowindex);
  free(zerorowindex);
  dd_clear(temp1);
  dd_clear(temp2);

 _L99:
  return Mnew;
}


/* end of cddproj.c */