/***************************************************************************
                          leastcore.c  -  description
                             -------------------
    begin                : Wed Mar 13 2002
    copyright            : (C) 2002 by Sven Klauke
    email                : sklauke@wiwi.uni-bielefeld.de
                           Institute of Mathematical Economics (IMW)
                           University of Bielefeld, Germany
                           http://www.wiwi.uni-bielefeld.de/~imw/
***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
#include "gmp.h"
#include "setoper.h"
#include "cdd.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "balance.h"
#include <string.h>
#include <time.h>

dd_boolean SetReadFile(FILE **f, dd_DataFileType fname)
{
  dd_boolean success=dd_FALSE;

  if ( (*f = fopen(fname, "r")) != NULL){
    //printf("output file %s is open\n",fname);
    success=dd_TRUE;
  }
  else{
    printf("The input file %s cannot be opened\n",fname);
  }
  return success;
}

dd_boolean SetWriteFile(FILE **f, dd_DataFileType fname)
{
  dd_boolean success=dd_FALSE;

  if ( (*f = fopen(fname, "w")) != NULL){
    //printf("output file %s is open\n",fname);
    success=dd_TRUE;
  }
  else{
    printf("The output file %s cannot be opened\n",fname);
  }
  return success;
}

int calc_leastcore(int nplayer,int* p[], int flag_redundant)
{
  int i,j;
  int ncoalition=(int)pow(2,nplayer);
  dd_ErrorType err=dd_NoError;
  dd_MatrixPtr A,B,B2;
  dd_NumberType numb;
  dd_rowrange m;
  dd_rowset redrows,linrows;
  dd_colrange n,d;
  dd_LPSolverType solver=dd_DualSimplex;
  dd_LPPtr lp;
  dd_LPSolutionPtr lps;
  mpq_t z;
  dd_PolyhedraPtr poly;
  dd_MatrixPtr vert;
  dd_DataFileType outputfile;
  FILE *writing;
  dd_Arow point;

  dd_set_global_constants();

  numb=dd_Rational;
  m=ncoalition-1;
  n=nplayer+2;
  A=dd_CreateMatrix(m,n);
  for(i=0;i<ncoalition-1;i++)
    {
      for(j=0;j<nplayer+2;j++)
    	dd_set_si(A->matrix[i][j],p[i][j]);
    }
  for(i=0;i<nplayer+1;i++)
    dd_set_si(A->rowvec[i],0);
  dd_set_si(A->rowvec[nplayer+1],1);
  set_addelem(A->linset,ncoalition-1);
  A->representation=dd_Inequality;
  A->objective=dd_LPmin;
  //printf("die least-core-berechnungs-matrix:\n");
  //dd_WriteMatrix(stdout,A);
  lp=dd_Matrix2LP(A,&err);
  dd_FreeMatrix(A);
  //printf("solving lp ... ");
  dd_LPSolve(lp,solver,&err);
  //printf("done\n");
  lps=dd_CopyLPSolution(lp);
  dd_FreeLPData(lp);
  m=ncoalition-1;
  n=nplayer+1;
  B=dd_CreateMatrix(m,n);
  mpq_init(z);
  for(i=0;i<ncoalition-2;i++)
    {
      mpq_set_si(z,p[i][0],1);
      mpq_canonicalize(z);
      mpq_add(z,z,lps->optvalue);
      mpq_canonicalize(z);
      dd_set(B->matrix[i][0],z);
      for(j=0;j<nplayer;j++)
    	dd_set_si(B->matrix[i][j+1],p[i][j+1]);
    }
  dd_set_si(B->matrix[ncoalition-2][0],p[ncoalition-2][0]);
  for(j=0;j<nplayer;j++)
    dd_set_si(B->matrix[ncoalition-2][j+1],-1);
  set_addelem(B->linset,ncoalition-1);
  B->representation=dd_Inequality;
  //SetWriteFile(&writing,"optsol");
  //dd_WriteMatrix(writing,B);
  //fclose(writing);
  //printf("optsol written\n");
  if(flag_redundant==1)
    {
      printf("identifying redundant rows\n");
      if (B->representation==dd_Generator) d=B->colsize+1; else d=B->colsize;
      redrows=dd_RedundantRows(B,&err);
      printf("identified these redundant rows:\n");
      set_fwrite(stdout,redrows);
      B2=dd_MatrixSubmatrix(B,redrows);
      linrows=dd_ImplicitLinearityRows(B2,&err);
      printf("implicit linearity rows (after removal of redundant rows):\n");
      set_fwrite(stdout,linrows);
      set_uni(B2->linset,B2->linset,linrows);
    }
  if(flag_redundant==1)
    {
      //printf("region of optimal solutions:\n");
      //dd_WriteMatrix(stdout,B2);
      poly=dd_DDMatrix2Poly(B2,&err);
      dd_FreeMatrix(B2);
    }
  else
    {
      //dd_WriteMatrix(stdout,B);
      poly=dd_DDMatrix2Poly(B,&err);
    }
  dd_FreeMatrix(B);
  vert=dd_CopyGenerators(poly);
  SetWriteFile(&writing, "lc");
  dd_WriteMatrix(writing,vert);
  fclose(writing);
  dd_FreeLPSolution(lps);
  dd_FreePolyhedra(poly);
  //dd_FreeMatrix(vert);
  return vert->rowsize;
}

int calc_lp(int nplayer,int ncoalition, int* p[], int flag_redundant)
{
  int i,j;
  dd_ErrorType err=dd_NoError;
  dd_MatrixPtr Mat_New_Ineq,B,C,Mat_Prev_Vertices,B2;
  dd_NumberType numb;
  dd_rowset redrows,linrows;
  dd_rowrange m;
  dd_colrange n,d;
  dd_LPSolverType solver=dd_DualSimplex;
  dd_LPPtr lp;
  dd_LPSolutionPtr lps;
  mpq_t z;
  dd_PolyhedraPtr poly,Poly_Prev_Vert;
  dd_MatrixPtr vert;
  dd_DataFileType outputfile;
  FILE *writing;
  FILE *reading;

  dd_set_global_constants();

  numb=dd_Rational;
  m=ncoalition;
  n=nplayer+2;
  Mat_New_Ineq=dd_CreateMatrix(m,n);
  for(i=0;i<ncoalition;i++)
    {
      for(j=0;j<nplayer+2;j++)
    	dd_set_si(Mat_New_Ineq->matrix[i][j],p[i][j]);
    }
  for(i=0;i<nplayer+1;i++)
    dd_set_si(Mat_New_Ineq->rowvec[i],0);
  dd_set_si(Mat_New_Ineq->rowvec[nplayer+1],1);
  set_addelem(Mat_New_Ineq->linset,ncoalition);
  Mat_New_Ineq->representation=dd_Inequality;
  Mat_New_Ineq->objective=dd_LPmin;
  // die ecken aus "lp" einlesen und als ungleichungen einfügen
  SetReadFile(&reading, "lc");
  //printf("reading previous vertices\n");
  Mat_Prev_Vertices=dd_PolyFile2Matrix(reading,&err);
  //printf("previous vertices:\n");
  //dd_WriteMatrix(stdout,Mat_Prev_Vertices);
  // zu prev_vertices muß noch die n+1. Koordinate zugefügt werden
  C=dd_CreateMatrix(Mat_Prev_Vertices->rowsize,Mat_Prev_Vertices->colsize+1);
  mpq_init(z);
  for(i=0;i<Mat_Prev_Vertices->rowsize;i++)
    {
      for(j=0;j<Mat_Prev_Vertices->colsize;j++)
    	mpq_set(C->matrix[i][j],Mat_Prev_Vertices->matrix[i][j]);
      dd_set_si(C->matrix[i][Mat_Prev_Vertices->colsize],0);
    }
  C->representation=dd_Generator;
  Poly_Prev_Vert=dd_DDMatrix2Poly(C,&err);
  C=dd_CopyInequalities(Poly_Prev_Vert);
  dd_MatrixRowRemove(&C,C->rowsize);
  dd_MatrixAppendTo(&Mat_New_Ineq,C);

  for(i=0;i<nplayer+1;i++)
    dd_set_si(Mat_New_Ineq->rowvec[i],0);
  dd_set_si(Mat_New_Ineq->rowvec[nplayer+1],1);
  set_addelem(Mat_New_Ineq->linset,ncoalition);
  Mat_New_Ineq->representation=dd_Inequality;
  Mat_New_Ineq->objective=dd_LPmin;
  //printf("solving lp\n");
  lp=dd_Matrix2LP(Mat_New_Ineq,&err);
  dd_FreeMatrix(Mat_New_Ineq);
  dd_LPSolve(lp,solver,&err);
  //printf("done\n");
  lps=dd_CopyLPSolution(lp);
  dd_FreeLPData(lp);
  m=ncoalition-1;
  n=nplayer+1;
  B=dd_CreateMatrix(m,n);
  for(i=0;i<ncoalition-1;i++)
    {
      mpq_set_si(z,p[i][0],1);
      mpq_canonicalize(z);
      mpq_add(z,z,lps->optvalue);
      mpq_canonicalize(z);
      dd_set(B->matrix[i][0],z);
      for(j=0;j<nplayer;j++)
    	dd_set_si(B->matrix[i][j+1],p[i][j+1]);
    }
  B->representation=dd_Inequality;
  poly=dd_DDMatrix2Poly(Mat_Prev_Vertices,&err);
  vert=dd_CopyInequalities(poly);
  dd_MatrixAppendTo(&B,vert);
  B->representation=dd_Inequality;
  if(flag_redundant==1)
    {
      printf("identifying redundant rows\n");
      if (B->representation==dd_Generator) d=B->colsize+1; else d=B->colsize;
      redrows=dd_RedundantRows(B,&err);
      printf("identified these redundant rows:\n");
      set_fwrite(stdout,redrows);
      B2=dd_MatrixSubmatrix(B,redrows);
      linrows=dd_ImplicitLinearityRows(B2,&err);
      printf("implicit linearity rows (after removal of redundant rows):\n");
      set_fwrite(stdout,linrows);
      set_uni(B2->linset,B2->linset,linrows);
    }
  if(flag_redundant==1)
    {
      printf("region of optimal solutions:\n");
      dd_WriteMatrix(stdout,B2);
      poly=dd_DDMatrix2Poly(B2,&err);
      dd_FreeMatrix(B2);
    }
  else
    {
      //dd_WriteMatrix(stdout,B);
      poly=dd_DDMatrix2Poly(B,&err);
    }
  dd_FreeMatrix(B);
  vert=dd_CopyGenerators(poly);
  //dd_WriteMatrix(stdout,vert);
  SetWriteFile(&writing, "lc");
  dd_WriteMatrix(writing,vert);
  fclose(writing);
  dd_FreeLPSolution(lps);
  dd_FreePolyhedra(poly);
  //dd_FreeMatrix(vert);
  return vert->rowsize;
}