#define COPYRIGHT   "Copyright (C) 1994, Komei Fukuda, fukuda@dma.epfl.ch"
#define DDVERSION   "Version C0.51c (March 15, 1994)"

/*  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., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/* cdd.c : C-Implementation of the double description method for
   computing all vertices and extremal rays of the polyhedron 
   P= {x :  b - A x >= 0}.  
   Please read COPYING (GNU General Public Licence) and
   the manual cddman.tex for details.
*/

/* The first version C0.21 was created , November 10, 1993 
   with Dave Gillespie's p2c translator 
   from the Pascal program pdd.p written by Komei Fukuda. 
*/

#include "setoper.h" /* set operation library header (Jan. 23,1994 version or later) */
#include "cdddef.h"
#include "cdd.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
/* #include <profile.h>   THINK C PROFILER */
/* #include <console.h>   THINK C PROFILER */

long minput, ninput;   /*size of input data [b -A] */
long mm, nn;   /*size of the homogenous system to be solved by dd*/
long projdim;  /*dimension of orthogonal preprojection */
colset projvars;   /*set of variables spanning the space of preprojection, 
     i.e. the remaining variables are to be removed*/
rowset MarkedSet, GroundSet, Face, Face1;
rowrange Iteration, hh;
rowindex OrderVector;  /* the permutation vector to store a preordered row indeces */
rowset AddedHyperplanes, InitialHyperplanes;
long RayCount, FeasibleRayCount, TotalRayCount, VertexCount, EdgeCount, ZeroRayCount;
long count_int=0,count_int_good=0,count_int_bad=0;
boolean DynamicWriteOn, DynamicRayWriteOn, LogWriteOn, debug;
Amatrix AA;
Bmatrix InitialRays;
colrange RHScol;   /* LP RHS column */
rowrange OBJrow;   /* LP OBJ row */
LPStatusType LPStatus;
Arow LPcost;  /* LP cost vector to be maximized  */
RayRecord *ArtificialRay, *FirstRay, *LastRay;
RayRecord *PosHead, *ZeroHead, *NegHead, *PosLast, *ZeroLast, *NegLast;
AdjacencyRecord *Edges[MMAX];  /* adjacency relation storage for iteration k */
boolean RecomputeRowOrder, found, inputsuccessful;
HyperplaneOrderType HyperplaneOrder;
AdjacencyTestType AdjacencyTest;
NumberType Number;
InequalityType Inequality;
boolean NondegAssumed;   /* Nondegeneacy preknowledge flag */
boolean InitBasisAtBottom;  /* if it is on, the initial Basis will be selected at bottom */
boolean PartialEnumeration; /* Partial enumeration Switch (TRUE if it is restricted on the intersection of MarkedSet hyperplanes) */
boolean PreOrderedRun; 
  /* TRUE if the rows are ordered before execution & all necessary adjacencies are stored */
CompStatusType CompStatus;  /* Computation Status */
ConversionType Conversion;
IncidenceOutputType IncidenceOutput;
AdjacencyOutputType AdjacencyOutput;
ErrorType Error;
FileInputModeType FileInputMode;
DataFileType inputfile,ifilehead,ifiletail,
  outputfile,projfile,icdfile,adjfile,logfile;
FILE *reading, *writing, *writing_proj, *writing_icd, *writing_adj,*writing_log;
time_t starttime, endtime;

void DefaultOptionSetup(void)
{
  debug = FALSE;
  DynamicWriteOn = TRUE;
  DynamicRayWriteOn = TRUE;
  LogWriteOn = FALSE;
  HyperplaneOrder = LexMin;
  AdjacencyTest = Combinatorial;
  NondegAssumed = FALSE;
  RecomputeRowOrder=TRUE;
  PreOrderedRun=TRUE;
  Conversion = IneToExt;
    
  IncidenceOutput = IncOff;
  AdjacencyOutput = AdjOff;
  InitBasisAtBottom = FALSE;
}

void SetInputFile(FILE **f, boolean *success)
{
  boolean opened=FALSE,stop, quit=FALSE;
  int i,dotpos=0, semipos=0, ifilelen=0;
  char ch;
  char *tempname;
  
  *success=FALSE;
  while (!opened && !quit) {
    if (FileInputMode!=Auto){
      printf("\n>> Input file (*.ine) : ");
      scanf("%s",inputfile);
      ch=getchar();
    }
    stop=FALSE;
    for (i=0; i<filenamelen && !stop; i++){
      ch=inputfile[i];
      switch (ch) {
        case '.': 
          dotpos=i+1;
          break;
        case ';':  case ' ':  case '\0':  case '\n':  case '\t':     
          if (ch==';'){
            semipos=i+1;
            FileInputMode=SemiAuto;   
            /* semicolon at the end of the filename
            -> output file names will be creeated with defaults. */
          }
          stop=TRUE;
          tempname=(char*)calloc(filenamelen,sizeof ch);
          strncpy(tempname, inputfile, i);
          strcpy(inputfile,tempname);
          break;
      }
    }
    if (dotpos>0){
      strncpy(ifilehead, inputfile, dotpos-1);
    }
    if (debug){
      printf("inputfile name: %s\n", inputfile);  
      printf("inputfile name head: %s\n", ifilehead);  
      printf("semicolon pos: %d\n", semipos);
    }  
    if ( ( *f = fopen(inputfile, "r") )!= NULL) {
      if (DynamicWriteOn) printf("input file %s is open\n", inputfile);
      opened=TRUE;
      *success=TRUE;
    }
    else{
      printf("The file %s not found\n",inputfile);
      if (FileInputMode==Auto) {
        quit=TRUE;
      }
    }
  }
}

void SetWriteFile(FILE **f, DataFileType fname, char cflag, char *fscript)
{
  boolean quit=FALSE;
  char *extension;
  
  switch (cflag) {
    case 'o':
      extension=".ext";break;
    case 'a':
      extension=".adj";break;
    case 'i':
      extension=".icd";break;
    case 'l':
      extension=".ddl";break;
    case 'p':
      extension="sub.ine";break;
    default:
      extension=".xxx";break;
  }
  if (FileInputMode==Manual){
    while (!quit) {
      printf("\n>> %s file name (*%s)   : ",fscript, extension);
      scanf("%s",fname);
      if (fname[0]==';'|| fname[0]<'0'){
        quit=TRUE;  /* default file name */
      } 
      else if (strcmp(inputfile, fname)!=0){
        *f = fopen(fname, "w");
        if (DynamicWriteOn) printf("write file %s is open\n",fname);
        goto _L99;
      }
      else {
        printf("%s file %s must have a name different from inputfile.\n",fscript,fname);
      }
    }
  }
  /* Auto or SemiAuto FileInput */
  strcpy(fname,ifilehead); 
  strcat(fname,extension); 
  if (strcmp(inputfile, fname)==0) {
    strcpy(fname,inputfile); 
    strcat(fname,extension); 
  }
  *f = fopen(fname, "w");
  if (DynamicWriteOn) printf("%s file %s is open\n",fscript,fname);
_L99:;
}


void SetNumberType(char *line)
{
  if (strncmp(line, "integer", 7)==0) {
    Number = Integer;
    return;
  }
  else if (strncmp(line, "rational", 8)==0) {
    Number = Rational;
    Error=ImproperInputFormat;  /* Rational Input not supported */
    return;
  }
  else if (strncmp(line, "real", 4)==0) {
    Number = Integer;
    return;
  }
  else { 
    Number=Unknown;
    Error=ImproperInputFormat;
  }
}


long Cardinality(long *S)
{
  rowrange i;
  long car;

  car = 0;
  for (i = 1; i <= mm; i++) {
    if (set_member(i, S))
      car++;
  }
  return car;
}


void WriteReal(FILE *f, double x)
{
  long ix1,ix2,ix;

  ix1= fabs(x) * 10000. + 0.5;
  ix2= (fabs(x) + 0.5);
  ix2= ix2*10000;
  if ( ix1 == ix2) {
    if (x>0) {
      ix = x + 0.5;
    } else {
      ix = -x + 0.5;
      ix = -ix;
    }
    fprintf(f, " %2ld", ix);
  } else
    fprintf(f, " % .9E", x);
}


void WriteSetElements(FILE *f, long *S)
{
  rowrange i;

  for (i = 1; i <= mm; i++) {
    if (set_member(i, S))
      fprintf(f, " %4ld", i);
  }
}


void WriteIncidence(FILE *f, RayRecord *RR)
{
  rowset cset;
  long zcar;

  set_initialize(&cset,mm);
  zcar = Cardinality(RR->ZeroSet);
  switch (IncidenceOutput) {

  case IncCardinality:
    fprintf(f, "%8ld", zcar);
    break;

  case IncSet:
    if (mm - zcar >= zcar) {
      fprintf(f, "%8ld", zcar);
      WriteSetElements(f, RR->ZeroSet);
    } else {
      set_diff(cset, GroundSet, RR->ZeroSet);
      fprintf(f, "%8ld", zcar - mm);
      WriteSetElements(f, cset);
    }
    break;

  case IncOff:
    break;
  }
  putc('\n', f);
  set_free(cset);
}



void CheckAdjacency1(RayRecord **RP1, RayRecord **RP2,
			    boolean *adjacent)
{
  long rank;

  *adjacent = TRUE;
  set_int(Face1, (*RP1)->ZeroSet, (*RP2)->ZeroSet);
  set_int(Face, Face1, AddedHyperplanes);
  if (debug)
    printf("Check adjacency\n");
  if (Cardinality(Face)< nn - 2) {
    *adjacent = FALSE;
    return;
  }
  else if (NondegAssumed) {
  	*adjacent = TRUE;
  	return;
  }
  ComputeRank(AA,Face,&rank);
  if (rank < nn - 2){
    *adjacent = FALSE;
  }
}


void CheckAdjacency2(RayRecord **RP1, RayRecord **RP2,
			    boolean *adjacent)
{
  RayRecord *TempRay;

  *adjacent = TRUE;
  set_int(Face1, (*RP1)->ZeroSet, (*RP2)->ZeroSet);
  set_int(Face, Face1, AddedHyperplanes);
  if (debug){
    printf("Check adjacency\n");
  }
  if (Cardinality(Face)< nn - 2) {
    *adjacent = FALSE;
    return;
  }
  else if (NondegAssumed) {
  	*adjacent = TRUE;
  	return;
  }
  TempRay = FirstRay;
  while (TempRay != NULL && *adjacent) {
    if (TempRay != *RP1 && TempRay != *RP2) {
    	set_int(Face1, TempRay->ZeroSet, AddedHyperplanes);
      	if (set_subset(Face, Face1))
			*adjacent = FALSE;
    }
    TempRay = TempRay->Next;
  }
}



void Eliminate(RayRecord **Ptr)
{
  /*eliminate the record pointed by Ptr^.Next*/
  RayRecord *TempPtr;

  if (debug) {
    printf("            Delete:");
    WriteRayRecord(stdout, (*Ptr)->Next);
  }
  TempPtr = (*Ptr)->Next;
  (*Ptr)->Next = (*Ptr)->Next->Next;
  if (TempPtr == FirstRay)   /*Update the first pointer*/
    FirstRay = (*Ptr)->Next;
  if (TempPtr == LastRay)   /*Update the last pointer*/
    LastRay = *Ptr;
  free(TempPtr->Ray);          /* free the ray vector memory */
  set_free(TempPtr->ZeroSet);  /* free the ZeroSet memory */
  free(TempPtr);   /* free the RayRecord structure memory */
  RayCount--; 
}


void SelectNextHyperplane0(long *excluded, rowrange *hnext)
{
  /*A natural way to choose the next hyperplane.  Simply the largest index*/
  long i;
  boolean determined;

  i = mm;
  determined = FALSE;
  do {
    if (set_member(i, excluded))
      i--;
    else
      determined = TRUE;
  } while (!determined && i>=1);
  if (determined) 
    *hnext = i;
  else
    *hnext = 0;
}


void SelectNextHyperplane1(long *excluded, rowrange *hnext)
{
  /*Natural way to choose the next hyperplane.  Simply the least index*/
  long i;
  boolean determined;

  i = 1;
  determined = FALSE;
  do {
    if (set_member(i, excluded))
      i++;
    else
      determined = TRUE;
  } while (!determined && i<=mm);
  if (determined) 
    *hnext = i;
  else 
    *hnext=0;
}


void SelectNextHyperplane2(long *excluded, rowrange *hnext)
{
  /*Choose the next hyperplane with maximum infeasibility*/
  long i, fea, inf, infmin, fi=0;   /*feasibility and infeasibility numbers*/

  infmin = RayCount + 1;
  for (i = 1; i <= mm; i++) {
    if (!set_member(i, excluded)) {
      FeasibilityIndices(&fea, &inf, i);
      if (inf < infmin) {
	infmin = inf;
	fi = fea;
	*hnext = i;
      }
    }
  }
  if (DynamicWriteOn) {
    printf("*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infmin, fi);
    fprintf(writing, "*infeasible rays (min) =%5ld, #feas rays =%5ld\n",
	    infmin, fi);
  }
}


void SelectNextHyperplane3(long *excluded, rowrange *hnext)
{
  /*Choose the next hyperplane with maximum infeasibility*/
  long i, fea, inf, infmax, fi=0;   /*feasibility and infeasibility numbers*/

  infmax = -1;
  for (i = 1; i <= mm; i++) {
    if (!set_member(i, excluded)) {
      FeasibilityIndices(&fea, &inf, i);
      if (inf > infmax) {
	infmax = inf;
	fi = fea;
	*hnext = i;
      }
    }
  }
  if (DynamicWriteOn) {
    printf("*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infmax, fi);
    fprintf(writing, "*infeasible rays (max) =%5ld, #feas rays =%5ld\n",
	    infmax, fi);
  }
}


void SelectNextHyperplane4(long *excluded, rowrange *hnext)
{
  /*Choose the next hyperplane with the most unbalanced cut*/
  long i, fea, inf, max, tmax, fi=0, infi=0;
      /*feasibility and infeasibility numbers*/

  max = -1;
  for (i = 1; i <= mm; i++) {
    if (!set_member(i, excluded)) {
      FeasibilityIndices(&fea, &inf, i);
      if (fea <= inf)
	tmax = inf;
      else
	tmax = fea;
      if (tmax > max) {
	max = tmax;
	fi = fea;
	infi = inf;
	*hnext = i;
      }
    }
  }
  if (!DynamicWriteOn)
    return;
  if (max == fi) {
    printf("*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infi, fi);
    fprintf(writing, "*infeasible rays (min) =%5ld, #feas rays =%5ld\n",
	    infi, fi);
  } else {
    printf("*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infi, fi);
    fprintf(writing, "*infeasible rays (max) =%5ld, #feas rays =%5ld\n",
	    infi, fi);
  }
}



void SelectNextHyperplane5(long *excluded, rowrange *hnext)
{
  /*Choose the next hyperplane which is lexico-min*/
  long i, minindex;
  double *v1, *v2;

  minindex = 0;
  v1 = NULL;
  for (i = 1; i <= mm; i++) {
    if (!set_member(i, excluded)) {
	  v2 = AA[i - 1];
      if (minindex == 0) {
	    minindex = i;
	    v1=v2;
      } else if (LexSmaller(v2,v1)) {
        minindex = i;
	    v1=v2;
      }
    }
  }
  *hnext = minindex;
}


void SelectNextHyperplane6(long *excluded, rowrange *hnext)
{
  /*Choose the next hyperplane which is lexico-max*/
  long i, maxindex;
  double *v1, *v2;

  maxindex = 0;
  v1 = NULL;
  for (i = 1; i <= mm; i++) {
    if (!set_member(i, excluded)) {
      v2= AA[i - 1];
      if (maxindex == 0) {
        maxindex = i;
        v1=v2;
      } else if (LexLarger(v2, v1)) {
        maxindex = i;
        v1=v2;
     }
    }
  }
  *hnext = maxindex;
}

long Partition(rowindex OV, long p, long r)
{
  double *x;
  long i,j,ovi;
  
  x=AA[OV[p]-1];
  i=p-1;
  j=r+1;
  while (TRUE){
    do{
      j--;
    } while (LexLarger(AA[OV[j]-1],x));
    do{
      i++;
    } while (LexSmaller(AA[OV[i]-1],x));
    if (i<j){
      ovi=OV[i];
      OV[i]=OV[j];
      OV[j]=ovi;
    }
    else{
      return j;
    }
  }
}

void QuickSort(rowindex OV, long p, long r)
{
  long q;
  
  if (p < r){
    q = Partition(OV, p, r);
    QuickSort(OV, p, q);
    QuickSort(OV, q+1, r);
  }
}

void ComputeRowOrderVector(rowindex OV, HyperplaneOrderType ho)
{
  long i,itemp;
  boolean exchanged=TRUE;
  
  OV[0]=0;
  switch (ho){
  case MaxIndex:
    for(i=1; i<=mm; i++) OV[i]=mm-i+1;
    break;

  case MinIndex: case MinCutoff:  case MixCutoff:  case MaxCutoff:
    for(i=1; i<=mm; i++) OV[i]=i;
    break;

  case LexMin:
    for(i=1; i<=mm; i++) OV[i]=i;
    while (exchanged){
      exchanged=FALSE;
      for (i=1;i<=mm-1; i++){
        if (LexLarger(AA[OV[i]-1], AA[OV[i+1]-1])){
          if (debug) printf("LexMin Order: %ld row and %ld row exchanged\n",OV[i],OV[i+1]);
          itemp=OV[i];
          OV[i]=OV[i+1];
          OV[i+1]=itemp;
          exchanged=TRUE;
        }
      }
    }
    break;

  case LexMax:
    for(i=1; i<=mm; i++) OV[i]=i;
    while (exchanged){
      exchanged=FALSE;
      for (i=1;i<=mm-1; i++){
        if (LexSmaller(AA[OV[i]-1], AA[OV[i+1]-1])){
          if (debug) printf("LexMax Order: %ld row and %ld row exchanged\n",OV[i],OV[i+1]);
          itemp=OV[i];
          OV[i]=OV[i+1]; 
          OV[i+1]=itemp;
          exchanged=TRUE;
        }
      }
    }
    break;
  }
}


void ComputeRowOrderVector2(rowindex OV, HyperplaneOrderType ho)
{
  long i,itemp;
  boolean exchanged=TRUE;
  
  OV[0]=0;
  switch (ho){
  case MaxIndex:
    for(i=1; i<=mm; i++) OV[i]=mm-i+1;
    break;

  case MinIndex: case MinCutoff:  case MixCutoff:  case MaxCutoff:
    for(i=1; i<=mm; i++) OV[i]=i;
    break;

  case LexMin:
    for(i=1; i<=mm; i++) OV[i]=i;
    QuickSort(OV, 1, mm);
    break;

  case LexMax:
    for(i=1; i<=mm; i++) OV[i]=i;
    QuickSort(OV, 1, mm);
    for(i=1; i<=mm/2;i++){   /* just reverse the order */
      itemp=OV[i];
      OV[i]=OV[mm-i+1];
      OV[mm-i+1]=itemp;
    }
    break;
  }
}

void UpdateRowOrderVector(long *PriorityRows)
/* Update the RowOrder vector to shift selected rows
in highest order.
*/
{
  rowrange i=0,i1,j=0,j1,ki=0,kj=0;
  boolean found, found2, localdebug=FALSE;
  
  if (localdebug) set_write(PriorityRows);
  found=TRUE;
  do {
    found=FALSE;
    for (i1=1; i1<=nn && !found; i1++){
      ki=OrderVector[i1];
      if (localdebug) printf("priority %ld,  row %ld\n", i1, ki);
      if (!set_member(ki, PriorityRows)){
        found=TRUE;
        i=i1;
      }
    }
    if (found){
      found2=FALSE;
      for (j1=i+1; j1<=mm && !found2 ; j1++){
        kj=OrderVector[j1];
        if (set_member(kj, PriorityRows)){
          found2=TRUE;
          j=j1;
        }
      }
      if (found2){
        if (localdebug){
          printf("exchanged with priority %ld  row %ld\n", j, kj);
        }
        OrderVector[i]=kj;
        OrderVector[j]=ki;
      }
      else{
        printf("Error: UpdateRowOrderVector\n");
        goto _L99;
      }
    }
  }while (found);
_L99:;
}

void SelectPreorderedNext(long *excluded, rowindex OV, rowrange *hnext)
{
  rowrange i,k;
  
  *hnext=0;
  for (i=1; i<=mm && *hnext==0; i++){
    k=OV[i];
    if (!set_member(k, excluded)) *hnext=k ;
  }
}

void SelectNextHyperplane(HyperplaneOrderType ho, 
         long *excluded, rowrange *hh, boolean *RefreshOrderVector)
{
  if (*RefreshOrderVector||OrderVector==NULL){
    if (OrderVector!=NULL) free(OrderVector);
    OrderVector=(long *)calloc(mm+1, sizeof *OrderVector);
    ComputeRowOrderVector(OrderVector, ho);
    *RefreshOrderVector=FALSE;
  }  
  switch (ho) {

  case MinCutoff:
    SelectNextHyperplane2(excluded,hh);
    break;

  case MaxCutoff:
    SelectNextHyperplane3(excluded, hh);
    break;

  case MixCutoff:
    SelectNextHyperplane4(excluded, hh);
    break;

  default:
    SelectPreorderedNext(excluded, OrderVector, hh);
    break;
  }
}



void WriteProgramDescription(FILE *f)
{
  fprintf(f, "* cdd: Double Description Method C-Code:%s\n", DDVERSION);
  fprintf(f,"* %s\n",COPYRIGHT);
}

void WriteRunningMode(FILE *f)
{
  if (Conversion==IneToExt || Conversion==ExtToIne){ 
    switch (HyperplaneOrder) {

    case MinIndex:
      fprintf(f, "*HyperplaneOrder: MinIndex\n");
      break;

    case MaxIndex:
      fprintf(f, "*HyperplaneOrder: MaxIndex\n");
      break;

    case MinCutoff:
      fprintf(f, "*HyperplaneOrder: MinCutoff\n");
      break;

    case MaxCutoff:
      fprintf(f, "*HyperplaneOrder: MaxCutoff\n");
    break;

    case MixCutoff:
      fprintf(f, "*HyperplaneOrder: MixCutoff\n");
      break;

    case LexMin:
      fprintf(f, "*HyperplaneOrder: LexMin\n");
      break;

    case LexMax:
      fprintf(f, "*HyperplaneOrder: LexMax\n");
      break;
    }
    if (NondegAssumed) {
      fprintf(f, "*Degeneracy preknowledge for computation: NondegenerateAssumed\n");
    }
    else {
      fprintf(f, "*Degeneracy preknowledge for computation: None (possible degeneracy)\n");
    }
  }
  switch (Conversion) {
    case ExtToIne:
      fprintf(f, "*Hull computation is chosen.\n");
      break;
    
    case IneToExt:
      fprintf(f, "*Vertex/Ray enumeration is chosen.\n");
      break;
    
    case LPmax:
      fprintf(f, "*LP (maximization) is chosen.\n");
      break;

    case LPmin:
      fprintf(f, "*LP (minimization) is chosen.\n");
      break;

    case Projection:
      fprintf(f, "*Preprojection is chosen.\n");
      break;
  
    default: break;
  }
  if (PartialEnumeration) {
    fprintf(f, "*Partial enumeration is chosen.  The permanently active rows are:");
    WriteSetElements(f,MarkedSet);
    fprintf(f,"\n");
  }
}

void WriteCompletionStatus(FILE *f)
{
  if (Iteration<mm && CompStatus==AllFound) {
    fprintf(f,"*Computation completed at Iteration %4ld.\n", Iteration);
  } 
  if (CompStatus == RegionEmpty) {
    fprintf(f,"*Computation completed at Iteration %4ld because the region found empty.\n", Iteration);
  }   
}


void WriteTimes(FILE *f)
{ 
  long ptime,ptime_sec,ptime_minu, ptime_hour;
  
  /* ptime=difftime(endtime,starttime); */   /* This function is ANSI standard, but not available sometime */
  ptime=endtime-starttime;      /* This is to replace the line above, but it may not give correct time in seconds */ 
  ptime_hour=ptime/3600;
  ptime_minu=(ptime-ptime_hour*3600)/60;
  ptime_sec=ptime%60;
  fprintf(f, "*Computation starts     at %s", asctime(localtime(&starttime)));
  fprintf(f, "*            terminates at %s", asctime(localtime(&endtime)));
  fprintf(f, "*Total processor time = %ld seconds\n", ptime);
  fprintf(f, "*                     = %ld hour %ld min %ld sec\n", ptime_hour,ptime_minu,ptime_sec);
}


void WriteAdjacency(FILE *f)
{
  RayRecord *RayPtr1, *RayPtr2;
  long pos1, pos2, degree;
  boolean adj;
  node *headnode, *tailnode, *newnode, *prevnode;

  headnode=NULL; tailnode=NULL;
  switch (Conversion) {
  case IneToExt:
    if (AdjacencyOutput==OutputAdjacency)
      fprintf(f,"*Adjacency List of output (=vertices/rays)\n");
    break;
  case ExtToIne:
    if (AdjacencyOutput==OutputAdjacency)
      fprintf(f,"*Adjacency List of output (=inequalities=facets)\n");
      break;

  default:
    break;
  }
  fprintf(f, "*cdd input file : %s   (%4ld  x %4ld)\n",
	  inputfile, minput, ninput);
  fprintf(f, "*cdd output file: %s\n", outputfile);
  fprintf(f,"begin\n");
  fprintf(f,"  %ld\n",RayCount);
  if (RayCount==0){
    goto _L99;
  }
  LastRay->Next=NULL;
  for (RayPtr1=FirstRay, pos1=1;RayPtr1 != NULL; RayPtr1 = RayPtr1->Next, pos1++){
    for (RayPtr2=FirstRay, pos2=1,degree=0; RayPtr2 != NULL; RayPtr2 = RayPtr2->Next, pos2++){
      if (RayPtr1!=RayPtr2){
        CheckAdjacency2(&RayPtr1, &RayPtr2, &adj);
        if (adj) {
          degree++;
          if (degree==1){
            newnode=(node *)malloc(sizeof *newnode);
            newnode->key=pos2;
            newnode->next=NULL;
            headnode=newnode;
            tailnode=newnode;
          }
          else{
            newnode=(node *)malloc(sizeof *newnode);
            newnode->key=pos2;
            newnode->next=NULL;
            tailnode->next=newnode;
            tailnode=newnode;
          }
        }
      }
    }
    fprintf(f," %ld  %ld :", pos1, degree);
    for (newnode=headnode; newnode!=NULL; newnode=newnode->next, free(prevnode)){
      prevnode=newnode;
      fprintf(f," %ld", newnode->key);
    }
    fprintf(f,"\n");
  }
_L99:;
  fprintf(f,"end\n");
}

void WriteDDResult(void)
{
  RayRecord *TempPtr;

  if (!debug) writing=freopen(outputfile,"w",writing);
  time(&endtime);
  WriteProgramDescription(writing);
  fprintf(writing, "*Input File:%.*s   (%4ld  x %4ld)\n",
	  filenamelen, inputfile, minput, ninput);
  WriteRunningMode(writing);
  WriteCompletionStatus(writing);
  WriteTimes(writing);
  if (Conversion == ExtToIne)
    fprintf(writing,
      "*Since hull computation is chosen, the output is a minimal inequality system\n");
  fprintf(writing, "*FINAL RESULT:\n");
  if (DynamicWriteOn)
    printf("*Computation complete.\n");
  if (Conversion == IneToExt) {
    if (DynamicWriteOn)
      printf("*Number of Vertices =%8ld,   Rays =%8ld\n",
	     VertexCount, RayCount - VertexCount);
    fprintf(writing, "*Number of Vertices =%8ld,   Rays =%8ld\n",
	    VertexCount, RayCount - VertexCount);
  } else {
    if (DynamicWriteOn)
      printf("*Number of Facets =%8ld\n", RayCount);
    fprintf(writing, "*Number of Facets =%8ld\n", RayCount);
  }
  fprintf(writing, "begin\n");
  switch (Inequality) {
  case ZeroRHS:
    fprintf(writing, " %8ld  %5ld    real\n", RayCount, nn + 1);
    break;
  case NonzeroRHS:
    fprintf(writing, " %8ld  %5ld    real\n", RayCount, nn);
    break;
  }
  if (IncidenceOutput == IncSet) {
    writing_icd=freopen(icdfile,"w",writing_icd);
    switch (Conversion) {
    case IneToExt:
      fprintf(writing_icd,
	    "*Incidences of output(=vertices/rays) and input (=hyperplanes)\n");
      fprintf(writing_icd,
        "*   for each output, #incidence and the set of hyperplanes containing it\n");
      fprintf(writing_icd,
	    "*   or its complement with its cardinality with minus sign\n");
      break;
    case ExtToIne:
      fprintf(writing_icd,
	    "*Incidences of output(=facets) and input (=points)\n");
      fprintf(writing_icd,
        "*   for each output, #incidence and the set of points lying on it\n");
      fprintf(writing_icd,
	    "*   or its complement with its cardinality with minus sign\n");
      break;

    default:
      break;
    }
    fprintf(writing_icd, "*cdd input file : %s   (%4ld  x %4ld)\n",
	  inputfile, minput, ninput);
    fprintf(writing_icd, "*cdd output file: %s\n", outputfile);
    fprintf(writing_icd, "begin\n");
    fprintf(writing_icd, "%8ld%5ld%5ld\n", RayCount, minput, mm);
  }
  TempPtr = FirstRay;
  while (TempPtr != NULL) {
    WriteRayRecord(writing, TempPtr);
    if (IncidenceOutput == IncSet)
      WriteIncidence(writing_icd, TempPtr);
    TempPtr = TempPtr->Next;
  }
  fprintf(writing, "end\n");
  if (DynamicWriteOn) {
    WriteCompletionStatus(stdout);
    WriteTimes(stdout);
  }
  if (LogWriteOn) {
    fprintf(writing_log, "end\n");
    WriteRunningMode(writing_log);
    WriteCompletionStatus(writing_log);
    fprintf(writing_log, "*set_intersection total#, effective#, loss# = %ld   %ld   %ld\n", 
      count_int, count_int_good, count_int_bad);
    WriteTimes(writing_log);
  }
  if (IncidenceOutput == IncSet)
    fprintf(writing_icd, "end\n");
  if (AdjacencyOutput != AdjOff){
    if (DynamicWriteOn) printf("Writing the adjacency file...\n");
    WriteAdjacency(writing_adj);
  }
}

void WriteProjRayRecord(FILE *f, RayRecord *RR, long *dbrow)
{
  long i,j,k;
  static double *vec;
  static rowset dbset;
  static long mprev=0, nprev=0;

  if (dbset==NULL || mprev<mm){
     set_initialize(&dbset,mm);  
     /* initialize only for the first time or when a larger space is needed */
     mprev=mm;
  }
  if (vec==NULL || nprev<nn) {
    vec=(double *)calloc(nn, sizeof *vec);
    nprev=nn;
  }
  for (j = 1; j <= mm-nn; j++){
    i=dbrow[j];
    set_addelem(dbset,i);
    if (debug) printf("index %ld is added to dbset\n",i);
    vec[i-1]=0;
    for (k=1; k<=nn; k++) {
      vec[i-1]+= (RR->Ray[k-1])*AA[j-1][k-1];
    }
    if (debug) printf("vec[ %ld]= %.5E \n",i-1, vec[i-1]);
  }
  i=1;
  for (j = 1; j <= mm; j++){
    if (!set_member(j,dbset)){
      vec[j-1]=RR->Ray[i-1];
      i++;
    }
  }
  fprintf(f, " %2d", 0);
  for (j = 0; j < mm; j++)
    WriteReal(f, vec[j]);
  putc('\n', f);
}



void WriteProjResult(long *dbrow)
{
  RayRecord *TempPtr;

  if (!debug) writing=freopen(outputfile,"w",writing);
  time(&endtime);
  WriteProgramDescription(writing);
  fprintf(writing, "*Input File:%.*s   (%4ld  x %4ld)\n",
	  filenamelen, inputfile, minput, ninput);
  WriteRunningMode(writing);
  WriteCompletionStatus(writing);
  WriteTimes(writing);
  fprintf(writing, "*FINAL RESULT:\n");
  if (DynamicWriteOn)
    printf("*Computation complete.\n");
  if (DynamicWriteOn)
     printf("*Number of Vertices =%8ld,   Rays =%8ld\n",
	   VertexCount, RayCount - VertexCount);
  fprintf(writing, "*Number of Vertices =%8ld,   Rays =%8ld\n",
	 VertexCount, RayCount - VertexCount);
  fprintf(writing, "begin\n");
  fprintf(writing, " %8ld  %5ld    real\n", RayCount, mm + 1);
  TempPtr = FirstRay;
  while (TempPtr != NULL) {
    WriteProjRayRecord(writing, TempPtr, dbrow);
    TempPtr = TempPtr->Next;
  }
  fprintf(writing, "end\n");
  if (DynamicWriteOn) {
    WriteCompletionStatus(stdout);
    WriteTimes(stdout);
  }
  if (LogWriteOn) {
    fprintf(writing_log, "end\n");
    WriteRunningMode(writing_log);
    WriteCompletionStatus(writing_log);
    WriteTimes(writing_log);
  }
  if (AdjacencyOutput != AdjOff){
    if (DynamicWriteOn) printf("Writing the adjacency file...\n");
    WriteAdjacency(writing_adj);
  }
}

void InitialWriting(void)
{
  if (LogWriteOn) {
    fprintf(writing_log, "*Input File:%.*s   (%4ld  x %4ld)\n",
	  filenamelen, inputfile, minput, ninput);
	fprintf(writing_log,"*Initial set of hyperplanes: ");
    WriteSetElements(writing_log, AddedHyperplanes);
    fprintf(writing_log,"\n");
    fprintf(writing_log, "begin\n");
    fprintf(writing_log, "%5ld %3d\n", mm - nn, 5);
  }
}

void DDInit(void)
{
  long i;

  set_initialize(&InitialHyperplanes,mm);
  set_initialize(&AddedHyperplanes,mm);
  set_initialize(&GroundSet, mm);
  set_initialize(&Face, mm);   /* used in CheckAdjacency  */
  set_initialize(&Face1, mm);  /* used in CheckAdjacency  */
  OrderVector=(long *)calloc(mm+1, sizeof *OrderVector);
  ComputeRowOrderVector(OrderVector, HyperplaneOrder);
  RecomputeRowOrder=FALSE;
  InitializeBmatrix(InitialRays);
  for (i = 1; i <= mm; i++){
    set_addelem(GroundSet, i);
  }
  EdgeCount=0; /*total edge count */
  ZeroRayCount=0;
  if (HyperplaneOrder==MaxCutoff||HyperplaneOrder==MinCutoff) PreOrderedRun=FALSE;
}

void DDMain(void)
{
  Iteration = nn + 1;
  while (Iteration <= mm) {
    SelectNextHyperplane(HyperplaneOrder, AddedHyperplanes, &hh, &RecomputeRowOrder);
    if (DynamicWriteOn) {
      fprintf(writing,
	      "*----------  Iteration =%3ld :   add  row # %3ld ----------\n",
	      Iteration, hh);
      printf("*----------  Iteration =%3ld :   add  row # %3ld ----------\n",
	     Iteration, hh);
    }
    if (PreOrderedRun) 
      AddNewHyperplane2(hh);
    else 
      AddNewHyperplane1(hh);
    if (CompStatus==AllFound||CompStatus==RegionEmpty) {
      set_addelem(AddedHyperplanes, hh);
      goto _L99;
    }
    if (LogWriteOn)
      fprintf(writing_log, "%3ld %5ld %6ld %6ld %6ld\n",
	      Iteration, hh, TotalRayCount, RayCount, FeasibleRayCount);
    set_addelem(AddedHyperplanes, hh);
    Iteration++;
  }
  _L99:;
}


void WriteErrorMessages(FILE *f)
{
  switch (Error) {

  case LowColumnRank:
    if (Conversion==IneToExt) {
      fprintf(f,"*Input Error: Input matrix (b, -A) is not column full rank => no vertices and rays.\n");
      break;
    } else {
      fprintf(f,"*Input Error: Input matrix A is not column full rank.=> no vertices and rays.\n");
      break;
    }
 
  case DimensionTooLarge:
    fprintf(f, "*Input Error: Input matrix is too large:\n");
    fprintf(f, "*Please increase MMAX and/or NMAX in the source code and recompile.\n");
    break;

  case DependentMarkedSet:
    fprintf(f, "*Input Error: Marked rows are linearly dependent.\n");
    fprintf(f, "*Please select independent rows for partial enumeration.\n");
    break;

  case FileNotFound:
    fprintf(f, "*Input Error: Specified input file does not exist.\n");
    break;

  case ImproperInputFormat:
    if (Number == Rational) {
     }
    else {
      fprintf(f,"*Input Error: Input format is not correct.\n");
      fprintf(f,"*Format:\n");
      fprintf(f," begin\n");
      fprintf(f,"   m   n  NumberType(real, rational or integer)\n");
      fprintf(f,"   b  -A\n");
      fprintf(f," end\n");
    }
    break;

  case None:
    fprintf(f,"*No Error found.\n");
    break;
  }
}

void Initialization(int ARGC, char *ARGV[])
/* Initialization of global variables */
{
  Error=None;
  CompStatus=InProgress;
  if (ARGC>1){
    FileInputMode=Auto;
    strcpy(inputfile,ARGV[1]);
  }
  else{
    FileInputMode=Manual;
  }
}

void InitialDataSetup(void)
{
  long j, r;
  Arow Vector;
  rowset ZSet;

  time(&starttime);
  RayCount = 0;
  TotalRayCount = 0;
  FeasibleRayCount = 0;
  VertexCount = 0;
  RecomputeRowOrder=FALSE;
  ArtificialRay = NULL;
  FirstRay = NULL;
  LastRay = NULL;
  set_initialize(&ZSet,mm);
  AddArtificialRay();
  Iteration = nn;   /*Initially,we have already  nn  hyperplanes */
  set_copy(AddedHyperplanes, InitialHyperplanes);
  UpdateRowOrderVector(InitialHyperplanes);
  for (r = 0; r < nn; r++) {
    for (j = 0; j < nn; j++)
      Vector[j] = InitialRays[j][r];
    Normalize(Vector);
    ZeroIndexSet(Vector, ZSet);
    if (set_subset(MarkedSet, ZSet)){
      if (debug) {
        printf("add an initial ray with zero set:");
        set_write(ZSet);
      }
      AddRay(Vector);
    }
  }
  CreateInitialEdges();
  set_free(ZSet);
}

void DDEnumerate(void)
{
  if (IncidenceOutput == IncSet)
    SetWriteFile(&writing_icd, icdfile, 'i', "incidence");
  if (AdjacencyOutput != AdjOff)
    SetWriteFile(&writing_adj, adjfile, 'a', "adjacency");
  if (LogWriteOn)
    SetWriteFile(&writing_log, logfile, 'l', "log");
  DDInit();
  FindInitialRays(InitialHyperplanes, InitialRays, &found);
  if (found) {
    InitialDataSetup();
    InitialWriting();
    DDMain();
    WriteDDResult();
  } else {
    WriteErrorMessages(stdout);
    WriteErrorMessages(writing);
  }
}

void PreProjection(void)
{
  rowset subrows1,subrows2,DBrows;
  colset subcols1,subcols2;  /* subcols1:projvars,  subcols2:rest */
  rowrange i;
  colrange j,k;
  colindex pivrow;
  Bmatrix DBinv;  /* dual basis matrix inverse */
  long DBrank;
  
  if (IncidenceOutput == IncSet)
    SetWriteFile(&writing_icd, icdfile, 'i', "incidence");
  if (AdjacencyOutput != AdjOff)
    SetWriteFile(&writing_adj, adjfile, 'a', "adjacency");
  if (LogWriteOn)
    SetWriteFile(&writing_log, logfile, 'l', "log");
  set_initialize(&subrows1,mm);
  set_initialize(&subrows2,mm);
  set_initialize(&DBrows,mm);
  set_initialize(&subcols1,nn);  /* subcol1 : projvar & RHS columns */
  set_initialize(&subcols2,nn);  /* subcol2 : remaining columns */
  SetWriteFile(&writing_proj, projfile, 'p', "preprojection variable subsystem");
  for (j=1;j<=nn;j++){
    if (set_member(j,projvars) || (j==1 && Inequality==NonzeroRHS))
      set_addelem(subcols1,j);
    else
      set_addelem(subcols2,j);
  }
  for (i=1; i<=mm; i++) set_addelem(subrows1,i);
  if (DynamicWriteOn){
    WriteSubMatrixOfAA(stdout,subrows1,subcols1,Inequality);
  }
  WriteSubMatrixOfAA(writing_proj,subrows1,subcols1,Inequality);
  Inequality=ZeroRHS;
  ReduceAA(subrows1,subcols2);
    /* Extract the submatrix of AA index by subcols2. 
       subcols2 is changed to a consecutive sequence starting from 1 */
  if (debug) {
    WriteAmatrix(stdout,AA,mm,nn,NonzeroRHS);
    WriteAmatrix(writing,AA,mm,nn,NonzeroRHS);
  }
  InitializeBmatrix(DBinv);
  FindBasis(AA,MinIndex,DBrows,pivrow,DBinv,&DBrank);
    /* DBrows stores the rows associated with a dual basis */
  if (debug){
    printf("rank of the new (deletion col) matrix is %ld\n", DBrank);
    printf("dual basis rows ="); set_write(DBrows);
  }
  for (j=1;j<=nn;j++) fprintf(writing,"pivot row at col %ld = %ld\n",j, pivrow[j]);
  set_diff(subrows2,subrows1,DBrows); 
    /* subrows2 stores the rows not in DBrows */
  for (j=1; j<=nn;j++){
    if (pivrow[j]==0) {
      set_delelem(subcols2,j);
      fprintf(writing,"Warning: col %ld is a linear combination of the other colums. The column linear dependency must be deleted for ray computation\n",j);
      for (k=j; k<=nn-1; k++){ /* shifting all pivrow information */
        pivrow[j]=pivrow[j+1];
      }
      pivrow[nn]=0;
      nn--;
    }
  }
  if (debug)  {
    printf("rows for ray enumeration:");set_write(subrows2);
    printf("cols for ray enumeration:");set_write(subcols2);
  }
  ReduceAA(subrows2,subcols2); 
    /* subrows2 is changed to a consecutive sequence starting from 1 */
  DualizeAA(DBinv);
  if (Error==DimensionTooLarge) goto _L99;
  if (debug) {
    WriteAmatrix(stdout,AA,mm,nn,ZeroRHS);
    WriteAmatrix(writing,AA,mm,nn,ZeroRHS);
  }
  if (DynamicWriteOn) {
    WriteRunningMode(stdout);
    WriteRunningMode(writing);
  }
  DDInit();
  FindInitialRays(InitialHyperplanes, InitialRays, &found);
  if (found) {
    InitialDataSetup();
    InitialWriting();
    DDMain();
    WriteProjResult(pivrow);
  } else {
    _L99:;
    WriteErrorMessages(stdout);
    WriteErrorMessages(writing);
  }
  set_free(subrows1);
  set_free(subrows2);
  set_free(DBrows);
  set_free(subcols1);
  set_free(subcols2);
}

void LPMain(void)
{
  colindex NBIndex;  /* NBIndex[s] stores the nonbasic variable in column s */ 
  Arow LPsol, LPdsol;  /*  LP solution and the dual solution (basic var only) */
  rowrange re;  /* evidence row when LP is inconsistent */
  colrange se;  /* evidence col when LP is dual-inconsistent */
  double ov;  /* LP optimum value */
  long LPiter;

  if (Inequality==ZeroRHS){
    printf("Sorry, LP optimization is not implemented for RHS==0.\n");
    goto _L99;
  }
  if (Conversion==LPmax){
    CrissCrossMaximize(AA, InitialRays, OBJrow, RHScol, 
      &LPStatus, &ov, LPsol, LPdsol,NBIndex, &re, &se, &LPiter);
  }
  else if (Conversion==LPmin){
    CrissCrossMinimize(AA, InitialRays, OBJrow, RHScol, 
      &LPStatus, &ov, LPsol, LPdsol,NBIndex, &re, &se, &LPiter);
  }
  WriteLPResult(writing, LPStatus, ov, LPsol, LPdsol, NBIndex, re, se, LPiter);
  if (DynamicWriteOn)
    WriteLPResult(stdout,LPStatus, ov, LPsol, LPdsol, NBIndex, re, se, LPiter);
_L99:;
}

void InteriorFindMain(void)
{
  colindex NBIndex;  /* NBIndex[s] stores the nonbasic variable in column s */ 
  Arow LPsol, LPdsol;  /*  LP solution and the dual solution (basic var only) */
  rowrange re;  /* evidence row when LP is inconsistent */
  colrange se;  /* evidence col when LP is dual-inconsistent */
  double ov;  /* LP optimum value */
  long LPiter;

  if (Inequality==ZeroRHS){
    printf("Sorry, find_interior is not implemented for RHS==0.\n");
    goto _L99;
  }
  EnlargeAAforInteriorFinding();
  OBJrow=mm; RHScol=1;
  CrissCrossMaximize(AA, InitialRays, OBJrow, RHScol, 
    &LPStatus, &ov, LPsol, LPdsol,NBIndex, &re, &se, &LPiter);
  WriteLPResult(writing, LPStatus, ov, LPsol, LPdsol, NBIndex, re, se, LPiter);
  if (DynamicWriteOn)
    WriteLPResult(stdout,LPStatus, ov, LPsol, LPdsol, NBIndex, re, se, LPiter);
_L99:;
}

void WriteHeading(void)
{
  WriteProgramDescription(stdout);
  printf("------------------------------------------\n");
  printf(" Vertex & Extremal Ray Enumeration for\n");
  printf(" the polyhedron P = { x :  b - A x >= 0 }\n");
  printf("------------------------------------------\n");
}


void main(int argc, char *argv[])
{
  writing_log = NULL;
  writing_icd = NULL;
  writing = NULL;
  reading = NULL;
  WriteHeading();
  DefaultOptionSetup();
  Initialization(argc, argv);
  AmatrixInput(&inputsuccessful);

 /* InitProfile(200,200);                THINK C PROFILER */
 /* cecho2file("cdd profile", 0,stdout); THINK C PROFILER */
 /* _trace=0;                            THINK C PROFILER */

  if (inputsuccessful) {
    SetWriteFile(&writing,outputfile,'o',"output");
    if (DynamicWriteOn) {
      WriteRunningMode(stdout);
      WriteRunningMode(writing);
    }
    switch (Conversion) {
    case ExtToIne: case IneToExt: /* vertex/facets enumeration is chosen */
      DDEnumerate();
      break;
    
    case LPmax:  case LPmin:      /* LP is chosen */
      LPMain();
      break;

    case Projection:              /* preprojection is chosen */
      PreProjection();
      break;

   case InteriorFind:      /* Interior point search is chosen */
      InteriorFindMain();
      break;
  
    default: break;
    }
  } else {
    WriteErrorMessages(stdout);
    if (writing!=NULL) WriteErrorMessages(writing);
  }
  if (writing != NULL){
    fclose(writing);
    if (DynamicWriteOn) printf("closing the file %s\n",outputfile);
  }
  if (writing_icd != NULL){
    fclose(writing_icd);
    if (DynamicWriteOn) printf("closing the file %s\n",icdfile);
  }
  if (writing_adj != NULL){
    fclose(writing_adj);
    if (DynamicWriteOn) printf("closing the file %s\n",adjfile);
  }
  if (writing_log != NULL){
    fclose(writing_log);
    if (DynamicWriteOn) printf("closing the file %s\n",logfile);
  }
  /* DumpProfile();   THINK C PROFILER */
  /* exit(0);         THINK C PROFILER */
}

/* end of cdd.c */

