cddio.c

CheckMate is a MATLAB-based tool for modeling, simulating and investigating properties of hybrid dyn

/* cddio.c:  Basic Input and Output Procedures for cddlib
   written by Komei Fukuda, [e-mail]
   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>

/* void dd_fread_rational_value (FILE *, mytype *); */
void dd_SetLinearity(dd_MatrixPtr, char *);

void dd_SetInputFile(FILE **f,dd_DataFileType inputfile,dd_ErrorType *Error)
{
  int opened=0,stop,quit=0;
  int i,dotpos=0,trial=0;
  char ch;
  char *tempname;
  
  
  *Error=dd_NoError;
  while (!opened && !quit) {
    fprintf(stderr,"\n>> Input file: ");
    scanf("%s",inputfile);
    ch=getchar();
    stop=dd_FALSE;
    for (i=0; i<dd_filenamelen && !stop; i++){
      ch=inputfile[i];
      switch (ch) {
        case '.': 
          dotpos=i+1;
          break;
        case ';':  case ' ':  case '\0':  case '\n':  case '\t':     
          stop=dd_TRUE;
          tempname=(char*)calloc(dd_filenamelen,sizeof(ch));
          strncpy(tempname,inputfile,i);
          strcpy(inputfile,tempname);
          free(tempname);
          break;
      }
    }
    if ( ( *f = fopen(inputfile,"r") )!= NULL) {
      fprintf(stderr,"input file %s is open\n",inputfile);
      opened=1;
      *Error=dd_NoError;
    }
    else{
      fprintf(stderr,"The file %s not found\n",inputfile);
      trial++;
      if (trial>5) {
        *Error=dd_IFileNotFound;
        quit=1;
      }
    }
  }
}

void dd_SetWriteFileName(dd_DataFileType inputfile, dd_DataFileType outfile, char cflag, dd_RepresentationType rep)
{
  char *extension;
  dd_DataFileType ifilehead="";
  unsigned int i,dotpos;
  
  switch (cflag) {
    case 'o':
      switch (rep) {
        case dd_Generator:
          extension=".ine"; break;     /* output file for ine data */
        case dd_Inequality:
          extension=".ext"; break;     /* output file for ext data */
        default: 
          extension=".xxx";break;
      }
      break;

    case 'a':         /* decide for output adjacence */
      if (rep==dd_Inequality)
        extension=".ead";       /* adjacency file for ext data */
      else
        extension=".iad";       /* adjacency file for ine data */
      break;
    case 'i':         /* decide for output incidence */
      if (rep==dd_Inequality)
        extension=".ecd";       /* ext incidence file */
      else
        extension=".icd";       /* ine incidence file */
      break;
    case 'n':         /* decide for input incidence */
      if (rep==dd_Inequality)
        extension=".icd";       /* ine incidence file */
      else
        extension=".ecd";       /* ext incidence file */
      break;
    case 'j':        /* decide for input adjacency */
      if (rep==dd_Inequality)
        extension=".iad";       /* ine adjacency file */
      else
        extension=".ead";       /* ext adjacency file */
      break;
    case 'l':
      extension=".ddl";break;   /* log file */
    case 'd':
      extension=".dex";break;   /* decomposition output */
    case 'p':
      extension="sub.ine";break;  /* preprojection sub inequality file */
    case 'v':
      extension=".solved";break;  /* verify_input file */
    case 's':
      extension=".lps";break;   /* LP solution file */
    default:
      extension=".xxx";break;
  }
  dotpos=-1;
  for (i=0; i< strlen(inputfile); i++){
    if (inputfile[i]=='.') dotpos=i;
  }
  if (dotpos>1) strncpy(ifilehead, inputfile, dotpos);
  else strcpy(ifilehead,inputfile);
  if (strlen(inputfile)<=0) strcpy(ifilehead,"tempcdd");
  strcpy(outfile,ifilehead); 
  strcat(outfile,extension); 
  if (strcmp(inputfile, outfile)==0) {
    strcpy(outfile,inputfile); 
    strcat(outfile,extension); 
  }
/*  fprintf(stderr,"outfile name = %s\n",outfile);  */
}


dd_NumberType dd_GetNumberType(char *line)
{
  dd_NumberType nt;

  if (strncmp(line, "integer", 7)==0) {
    nt = dd_Integer;
  }
  else if (strncmp(line, "rational", 8)==0) {
    nt = dd_Rational;
  }
  else if (strncmp(line, "real", 4)==0) {
    nt = dd_Real;
  }
  else { 
    nt=dd_Unknown;
  }
  return nt;
}

void dd_ProcessCommandLine(FILE *f, dd_MatrixPtr M, char *line)
{
  char newline[dd_linelenmax];
  dd_colrange j;
  mytype value;

  dd_init(value);
  if (strncmp(line, "hull", 4)==0) {
    M->representation = dd_Generator;
  }
  if (strncmp(line, "debug", 5)==0) {
    dd_debug = dd_TRUE;
#ifdef GMPRATIONAL
    ddf_debug = ddf_TRUE;
#endif
  }
  if (strncmp(line, "partial_enum", 12)==0 ||
       strncmp(line, "equality", 8)==0  ||
       strncmp(line, "linearity", 9)==0 ) {
    fgets(newline,dd_linelenmax,f);    
    dd_SetLinearity(M,newline);
  }
  if (strncmp(line, "maximize", 8)==0 ||
      strncmp(line, "minimize", 8)==0) {
    if (strncmp(line, "maximize", 8)==0) M->objective=dd_LPmax;
    else M->objective=dd_LPmin;
    for (j = 1; j <= M->colsize; j++) {
    if (M->numbtype==dd_Real) {
#if !defined(GMPRATIONAL)
        double rvalue;
        fscanf(f, "%lf", &rvalue);
        dd_set_d(value, rvalue);
#endif
      } else {
        dd_fread_rational_value (f, value);
      }
      dd_set(M->rowvec[j - 1],value);
      if (dd_debug) {fprintf(stderr,"cost(%5ld) =",j); dd_WriteNumber(stderr,value);}
    }  /*of j*/
  }
  dd_clear(value);
}

dd_boolean dd_AppendMatrix2Poly(dd_PolyhedraPtr *poly, dd_MatrixPtr M)
{
  dd_boolean success=dd_FALSE;
  dd_MatrixPtr Mpoly,Mnew=NULL;
  dd_ErrorType err;

  if ((*poly)!=NULL && (*poly)->m >=0 && (*poly)->d>=0 &&
      (*poly)->d==M->colsize && M->rowsize>0){
    Mpoly=dd_CopyInput(*poly);
    Mnew=dd_AppendMatrix(Mpoly, M);
    dd_FreePolyhedra(*poly);
    *poly=dd_DDMatrix2Poly(Mnew,&err);
    dd_FreeMatrix(Mpoly);
    dd_FreeMatrix(Mnew);
    if (err==dd_NoError) success=dd_TRUE;
  }
  return success;
}

dd_MatrixPtr dd_MatrixCopy(dd_MatrixPtr M)
{
  dd_MatrixPtr Mcopy=NULL;
  dd_rowrange m;
  dd_colrange d;

  m= M->rowsize;
  d= M->colsize;
  if (m >=0 && d >=0){
    Mcopy=dd_CreateMatrix(m, d);
    dd_CopyAmatrix(Mcopy->matrix, M->matrix, m, d);
    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
    set_copy(Mcopy->linset,M->linset);
    Mcopy->numbtype=M->numbtype;
    Mcopy->representation=M->representation;
    Mcopy->objective=M->objective;
  }
  return Mcopy;
}

dd_MatrixPtr dd_CopyMatrix(dd_MatrixPtr M)
{
  return dd_MatrixCopy(M);
}

dd_MatrixPtr dd_MatrixNormalizedCopy(dd_MatrixPtr M)
{
  dd_MatrixPtr Mcopy=NULL;
  dd_rowrange m;
  dd_colrange d;

  m= M->rowsize;
  d= M->colsize;
  if (m >=0 && d >=0){
    Mcopy=dd_CreateMatrix(m, d);
    dd_CopyNormalizedAmatrix(Mcopy->matrix, M->matrix, m, d);
    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
    set_copy(Mcopy->linset,M->linset);
    Mcopy->numbtype=M->numbtype;
    Mcopy->representation=M->representation;
    Mcopy->objective=M->objective;
  }
  return Mcopy;
}


dd_MatrixPtr dd_MatrixAppend(dd_MatrixPtr M1, dd_MatrixPtr M2)
{
  dd_MatrixPtr M=NULL;
  dd_rowrange i, m,m1,m2;
  dd_colrange j, d,d1,d2;
  
  m1=M1->rowsize;
  d1=M1->colsize;
  m2=M2->rowsize;
  d2=M2->colsize;

  m=m1+m2;
  d=d1;

  if (d1>=0 && d1==d2 && m1>=0 && m2>=0){
    M=dd_CreateMatrix(m, d);
    dd_CopyAmatrix(M->matrix, M1->matrix, m1, d);
    dd_CopyArow(M->rowvec, M1->rowvec, d);
    for (i=0; i<m1; i++){
      if (set_member(i+1,M1->linset)) set_addelem(M->linset,i+1);
    }
    for (i=0; i<m2; i++){
       for (j=0; j<d; j++)
         dd_set(M->matrix[m1+i][j],M2->matrix[i][j]);
         /* append the second matrix */
       if (set_member(i+1,M2->linset)) set_addelem(M->linset,m1+i+1);
    }
    M->numbtype=M1->numbtype;
  }
  return M;
}

dd_MatrixPtr dd_MatrixNormalizedSortedCopy(dd_MatrixPtr M,dd_rowindex *newpos)  /* 094 */
{ 
  /* Sort the rows of Amatrix lexicographically, and return a link to this sorted copy. 
  The vector newpos is allocated, where newpos[i] returns the new row index
  of the original row i (i=1,...,M->rowsize). */
  dd_MatrixPtr Mcopy=NULL,Mnorm=NULL;
  dd_rowrange m,i;
  dd_colrange d;
  dd_rowindex roworder;

  /* if (newpos!=NULL) free(newpos); */
  m= M->rowsize;
  d= M->colsize;
  roworder=(long *)calloc(m+1,sizeof(long*));
  *newpos=(long *)calloc(m+1,sizeof(long*));
  if (m >=0 && d >=0){
    Mnorm=dd_MatrixNormalizedCopy(M);
    Mcopy=dd_CreateMatrix(m, d);
    for(i=1; i<=m; i++) roworder[i]=i;
    
    dd_RandomPermutation(roworder, m, 123);
    dd_QuickSort(roworder,1,m,Mnorm->matrix,d);

    dd_PermuteCopyAmatrix(Mcopy->matrix, Mnorm->matrix, m, d, roworder);
    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
    for(i=1; i<=m; i++) {
      if (set_member(roworder[i],M->linset)) set_addelem(Mcopy->linset, i);
      (*newpos)[roworder[i]]=i;
    }
    Mcopy->numbtype=M->numbtype;
    Mcopy->representation=M->representation;
    Mcopy->objective=M->objective;
    dd_FreeMatrix(Mnorm);
  }
  free(roworder);
  return Mcopy;
}

dd_MatrixPtr dd_MatrixUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos)
{ 
  /* Remove row duplicates, and return a link to this sorted copy. 
     Linearity rows have priority over the other rows.
     It is better to call this after sorting with dd_MatrixNormalizedSortedCopy.
     The vector newpos is allocated, where *newpos[i] returns the new row index
     of the original row i (i=1,...,M->rowsize).  *newpos[i] is negative if the original
     row is dominated by -*newpos[i] and eliminated in the new copy.
  */
  dd_MatrixPtr Mcopy=NULL;
  dd_rowrange m,i,uniqrows;
  dd_rowset preferredrows;
  dd_colrange d;
  dd_rowindex roworder;

  /* if (newpos!=NULL) free(newpos); */
  m= M->rowsize;
  d= M->colsize;
  preferredrows=M->linset;
  roworder=(long *)calloc(m+1,sizeof(long*));
  if (m >=0 && d >=0){
    for(i=1; i<=m; i++) roworder[i]=i;
    dd_UniqueRows(roworder, 1, m, M->matrix, d,preferredrows, &uniqrows);

    Mcopy=dd_CreateMatrix(uniqrows, d);
    dd_PermutePartialCopyAmatrix(Mcopy->matrix, M->matrix, m, d, roworder,1,m);
    dd_CopyArow(Mcopy->rowvec, M->rowvec, d);
    for(i=1; i<=m; i++) {
      if (roworder[i]>0 && set_member(i,M->linset)) set_addelem(Mcopy->linset, roworder[i]);
    }
    Mcopy->numbtype=M->numbtype;
    Mcopy->representation=M->representation;
    Mcopy->objective=M->objective;
  }
  *newpos=roworder;
  return Mcopy;
}


dd_MatrixPtr dd_MatrixNormalizedSortedUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos) /* 094 */
{ 
  /* Sort and remove row duplicates, and return a link to this sorted copy. 
     Linearity rows have priority over the other rows.
     It is better to call this after sorting with dd_MatrixNormalizedSortedCopy.
     The vector newpos is allocated, where *newpos[i] returns the new row index
     of the original row i (i=1,...,M->rowsize).  *newpos[i] is negative if the original
     row is dominated by -*newpos[i] and eliminated in the new copy.
  */
  dd_MatrixPtr M1=NULL,M2=NULL;
  dd_rowrange m,i;
  dd_colrange d;
  dd_rowindex newpos1=NULL,newpos1r=NULL,newpos2=NULL;

  /* if (newpos!=NULL) free(newpos); */
  m= M->rowsize;
  d= M->colsize;
  *newpos=(long *)calloc(m+1,sizeof(long*));  
  newpos1r=(long *)calloc(m+1,sizeof(long*));  
  if (m>=0 && d>=0){
    M1=dd_MatrixNormalizedSortedCopy(M,&newpos1);
    for (i=1; i<=m;i++) newpos1r[newpos1[i]]=i;  /* reverse of newpos1 */
    M2=dd_MatrixUniqueCopy(M1,&newpos2);
    set_emptyset(M2->linset);
    for(i=1; i<=m; i++) {
      if (newpos2[newpos1[i]]>0){
         printf("newpos1[%ld]=%ld, newpos2[newpos1[%ld]]=%ld\n",i,newpos1[i], i,newpos2[newpos1[i]]);
         if (set_member(i,M->linset)) set_addelem(M2->linset, newpos2[newpos1[i]]);
         (*newpos)[i]=newpos2[newpos1[i]];
      } else {
         (*newpos)[i]=-newpos1r[-newpos2[newpos1[i]]];
      }
    }
  dd_FreeMatrix(M1);free(newpos1);free(newpos2);free(newpos1r);
  }
  
  return M2;
}

dd_MatrixPtr dd_MatrixSortedUniqueCopy(dd_MatrixPtr M,dd_rowindex *newpos)  /* 094 */
{ 
  /* Same as dd_MatrixNormalizedSortedUniqueCopy except that it returns a unnormalized origial data
     with original ordering.
  */
  dd_MatrixPtr M1=NULL,M2=NULL;
  dd_rowrange m,i,k,ii;
  dd_colrange d;
  dd_rowindex newpos1=NULL,newpos1r=NULL,newpos2=NULL;

  /* if (newpos!=NULL) free(newpos); */
  m= M->rowsize;
  d= M->colsize;
  *newpos=(long *)calloc(m+1,sizeof(long*));  
  newpos1r=(long *)calloc(m+1,sizeof(long*));  
  if (m>=0 && d>=0){
    M1=dd_MatrixNormalizedSortedCopy(M,&newpos1);
    for (i=1; i<=m;i++) newpos1r[newpos1[i]]=i;  /* reverse of newpos1 */
    M2=dd_MatrixUniqueCopy(M1,&newpos2);
    dd_FreeMatrix(M1);
    set_emptyset(M2->linset);
    for(i=1; i<=m; i++) {
      if (newpos2[newpos1[i]]>0){
         if (set_member(i,M->linset)) set_addelem(M2->linset, newpos2[newpos1[i]]);
         (*newpos)[i]=newpos2[newpos1[i]];
      } else {
         (*newpos)[i]=-newpos1r[-newpos2[newpos1[i]]];
      }
    }

    ii=0;
    set_emptyset(M2->linset);
    for (i = 1; i<=m ; i++) {
      k=(*newpos)[i];
      if (k>0) {
        ii+=1;
        (*newpos)[i]=ii;
        dd_CopyArow(M2->matrix[ii-1],M->matrix[i-1],d);
        if (set_member(i,M->linset)) set_addelem(M2->linset, ii);
      }
    }

    free(newpos1);free(newpos2);free(newpos1r);
  }
  
  return M2;
}

dd_MatrixPtr dd_AppendMatrix(dd_MatrixPtr M1, dd_MatrixPtr M2)
{
  return dd_MatrixAppend(M1,M2);
}

int dd_MatrixAppendTo(dd_MatrixPtr *M1, dd_MatrixPtr M2)
{
  dd_MatrixPtr M=NULL;
  dd_rowrange i, m,m1,m2;
  dd_colrange j, d,d1,d2;
  dd_boolean success=0;
  
  m1=(*M1)->rowsize;
  d1=(*M1)->colsize;
  m2=M2->rowsize;
  d2=M2->colsize;

  m=m1+m2;
  d=d1;

  if (d1>=0 && d1==d2 && m1>=0 && m2>=0){
    M=dd_CreateMatrix(m, d);
    dd_CopyAmatrix(M->matrix, (*M1)->matrix, m1, d);
    dd_CopyArow(M->rowvec, (*M1)->rowvec, d);
    for (i=0; i<m1; i++){
      if (set_member(i+1,(*M1)->linset)) set_addelem(M->linset,i+1);
    }
    for (i=0; i<m2; i++){
       for (j=0; j<d; j++)
         dd_set(M->matrix[m1+i][j],M2->matrix[i][j]);
         /* append the second matrix */
       if (set_member(i+1,M2->linset)) set_addelem(M->linset,m1+i+1);
    }
    M->numbtype=(*M1)->numbtype;
    dd_FreeMatrix(*M1);
    *M1=M;
    success=1;
  }
  return success;
}

int dd_MatrixRowRemove(dd_MatrixPtr *M, dd_rowrange r) /* 092 */
{
  dd_rowrange i,m;
  dd_colrange d;
  dd_boolean success=0;