cddio.c

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

  
  m=(*M)->rowsize;
  d=(*M)->colsize;

  if (r >= 1 && r <=m){
    (*M)->rowsize=m-1;
    dd_FreeArow(d, (*M)->matrix[r-1]);
    set_delelem((*M)->linset,r);
    /* slide the row headers */
    for (i=r; i<m; i++){
      (*M)->matrix[i-1]=(*M)->matrix[i];
      if (set_member(i+1, (*M)->linset)){
        set_delelem((*M)->linset,i+1);
        set_addelem((*M)->linset,i);
      }
    }
    success=1;
  }
  return success;
}

int dd_MatrixRowRemove2(dd_MatrixPtr *M, dd_rowrange r, dd_rowindex *newpos) /* 094 */
{
  dd_rowrange i,m;
  dd_colrange d;
  dd_boolean success=0;
  dd_rowindex roworder;
  
  m=(*M)->rowsize;
  d=(*M)->colsize;

  if (r >= 1 && r <=m){
    roworder=(long *)calloc(m+1,sizeof(long*));
    (*M)->rowsize=m-1;
    dd_FreeArow(d, (*M)->matrix[r-1]);
    set_delelem((*M)->linset,r);
    /* slide the row headers */
    for (i=1; i<r; i++) roworder[i]=i;
    roworder[r]=0; /* meaning it is removed */
    for (i=r; i<m; i++){
      (*M)->matrix[i-1]=(*M)->matrix[i];
      roworder[i+1]=i;
      if (set_member(i+1, (*M)->linset)){
        set_delelem((*M)->linset,i+1);
        set_addelem((*M)->linset,i);
      }
    }
    success=1;
  }
  return success;
}

dd_MatrixPtr dd_MatrixSubmatrix(dd_MatrixPtr M, dd_rowset delset) /* 092 */
{
  dd_MatrixPtr Msub=NULL;
  dd_rowrange i,isub=1, m,msub;
  dd_colrange d;
 
  m= M->rowsize;
  d= M->colsize;
  msub=m;
  if (m >=0 && d >=0){
    for (i=1; i<=m; i++) {
       if (set_member(i,delset)) msub-=1;
    }
    Msub=dd_CreateMatrix(msub, d);
    for (i=1; i<=m; i++){
      if (!set_member(i,delset)){
        dd_CopyArow(Msub->matrix[isub-1], M->matrix[i-1], d);
        if (set_member(i, M->linset)){
          set_addelem(Msub->linset,isub);
        }
        isub++;
      }
    }
    dd_CopyArow(Msub->rowvec, M->rowvec, d);
    Msub->numbtype=M->numbtype;
    Msub->representation=M->representation;
    Msub->objective=M->objective;
  }
  return Msub;
}

dd_MatrixPtr dd_MatrixSubmatrix2(dd_MatrixPtr M, dd_rowset delset,dd_rowindex *newpos) /* 092 */
{ /* returns a pointer to a new matrix which is a submatrix of M with rows in delset
  removed.  *newpos[i] returns the position of the original row i in the new matrix.
  It is -1 if and only if it is deleted.
  */
  
  dd_MatrixPtr Msub=NULL;
  dd_rowrange i,isub=1, m,msub;
  dd_colrange d;
  dd_rowindex roworder;

  m= M->rowsize;
  d= M->colsize;
  msub=m;
  if (m >=0 && d >=0){
    roworder=(long *)calloc(m+1,sizeof(long*));
    for (i=1; i<=m; i++) {
       if (set_member(i,delset)) msub-=1;
    }
    Msub=dd_CreateMatrix(msub, d);
    for (i=1; i<=m; i++){
      if (set_member(i,delset)){
        roworder[i]=0; /* zero means the row i is removed */
      } else {
        dd_CopyArow(Msub->matrix[isub-1], M->matrix[i-1], d);
        if (set_member(i, M->linset)){
          set_addelem(Msub->linset,isub);
        }
        roworder[i]=isub;
        isub++;
      }
    }
    *newpos=roworder;
    dd_CopyArow(Msub->rowvec, M->rowvec, d);
    Msub->numbtype=M->numbtype;
    Msub->representation=M->representation;
    Msub->objective=M->objective;
  }
  return Msub;
}


dd_MatrixPtr dd_MatrixSubmatrix2L(dd_MatrixPtr M, dd_rowset delset,dd_rowindex *newpos) /* 094 */
{ /* This is same as dd_MatrixSubmatrix2 except that the linearity rows will be shifted up
     so that they are at the top of the matrix.
  */
  dd_MatrixPtr Msub=NULL;
  dd_rowrange i,iL, iI, m,msub;
  dd_colrange d;
  dd_rowindex roworder;

  m= M->rowsize;
  d= M->colsize;
  msub=m;
  if (m >=0 && d >=0){
    roworder=(long *)calloc(m+1,sizeof(long*));
    for (i=1; i<=m; i++) {
       if (set_member(i,delset)) msub-=1;
    }
    Msub=dd_CreateMatrix(msub, d);
    iL=1; iI=set_card(M->linset)+1;  /* starting positions */
    for (i=1; i<=m; i++){
      if (set_member(i,delset)){
        roworder[i]=0; /* zero means the row i is removed */
      } else {
        if (set_member(i,M->linset)){
          dd_CopyArow(Msub->matrix[iL-1], M->matrix[i-1], d);
          set_delelem(Msub->linset,i);
          set_addelem(Msub->linset,iL);
          roworder[i]=iL;
          iL+=1;
        } else {
          dd_CopyArow(Msub->matrix[iI-1], M->matrix[i-1], d);
          roworder[i]=iI;
          iI+=1;
        }
       }
    }
    *newpos=roworder;
    dd_CopyArow(Msub->rowvec, M->rowvec, d);
    Msub->numbtype=M->numbtype;
    Msub->representation=M->representation;
    Msub->objective=M->objective;
  }
  return Msub;
}

int dd_MatrixRowsRemove(dd_MatrixPtr *M, dd_rowset delset) /* 094 */
{
  dd_MatrixPtr Msub=NULL;
  int success;
  
  Msub=dd_MatrixSubmatrix(*M, delset);
  dd_FreeMatrix(*M);
  *M=Msub;
  success=1;
  return success;
}

int dd_MatrixRowsRemove2(dd_MatrixPtr *M, dd_rowset delset,dd_rowindex *newpos) /* 094 */
{
  dd_MatrixPtr Msub=NULL;
  int success;
  
  Msub=dd_MatrixSubmatrix2(*M, delset,newpos);
  dd_FreeMatrix(*M);
  *M=Msub;
  success=1;
  return success;
}

int dd_MatrixShiftupLinearity(dd_MatrixPtr *M,dd_rowindex *newpos) /* 094 */
{
  dd_MatrixPtr Msub=NULL;
  int success;
  dd_rowset delset;
  
  set_initialize(&delset,(*M)->rowsize);  /* emptyset */
  Msub=dd_MatrixSubmatrix2L(*M, delset,newpos);
  dd_FreeMatrix(*M);
  *M=Msub;
  
  set_free(delset);
  success=1;
  return success;
}

dd_PolyhedraPtr dd_CreatePolyhedraData(dd_rowrange m, dd_colrange d)
{
  dd_rowrange i;
  dd_PolyhedraPtr poly=NULL;

  poly=(dd_PolyhedraPtr) malloc (sizeof(dd_PolyhedraType));
  poly->child       =NULL; /* this links the homogenized cone data */
  poly->m           =m;
  poly->d           =d;  
  poly->n           =-1;  /* the size of output is not known */
  poly->m_alloc     =m+2; /* the allocated row size of matrix A */
  poly->d_alloc     =d;   /* the allocated col size of matrix A */
  poly->numbtype=dd_Real;
  dd_InitializeAmatrix(poly->m_alloc,poly->d_alloc,&(poly->A));
  dd_InitializeArow(d,&(poly->c));           /* cost vector */
  poly->representation       =dd_Inequality;
  poly->homogeneous =dd_FALSE;

  poly->EqualityIndex=(int *)calloc(m+2, sizeof(int));  
    /* size increased to m+2 in 092b because it is used by the child cone, 
       This is a bug fix suggested by Thao Dang. */
    /* ith component is 1 if it is equality, -1 if it is strict inequality, 0 otherwise. */
  for (i = 0; i <= m+1; i++) poly->EqualityIndex[i]=0;

  poly->IsEmpty                 = -1;  /* initially set to -1, neither TRUE nor FALSE, meaning unknown */
  
  poly->NondegAssumed           = dd_FALSE;
  poly->InitBasisAtBottom       = dd_FALSE;
  poly->RestrictedEnumeration   = dd_FALSE;
  poly->RelaxedEnumeration      = dd_FALSE;

  poly->AincGenerated=dd_FALSE;  /* Ainc is a set array to store the input incidence. */

  return poly;
}

dd_boolean dd_InitializeConeData(dd_rowrange m, dd_colrange d, dd_ConePtr *cone)
{
  dd_boolean success=dd_TRUE;
  dd_colrange j;

  (*cone)=(dd_ConePtr)calloc(1, sizeof(dd_ConeType));

/* INPUT: A given representation of a cone: inequality */
  (*cone)->m=m;
  (*cone)->d=d;
  (*cone)->m_alloc=m+2; /* allocated row size of matrix A */
  (*cone)->d_alloc=d;   /* allocated col size of matrix A, B and Bsave */
  (*cone)->numbtype=dd_Real;
  (*cone)->parent=NULL;

/* CONTROL: variables to control computation */
  (*cone)->Iteration=0;

  (*cone)->HalfspaceOrder=dd_LexMin;

  (*cone)->ArtificialRay=NULL;
  (*cone)->FirstRay=NULL;
  (*cone)->LastRay=NULL; /* The second description: Generator */
  (*cone)->PosHead=NULL;
  (*cone)->ZeroHead=NULL;
  (*cone)->NegHead=NULL;
  (*cone)->PosLast=NULL;
  (*cone)->ZeroLast=NULL;
  (*cone)->NegLast=NULL;
  (*cone)->RecomputeRowOrder  = dd_TRUE;
  (*cone)->PreOrderedRun      = dd_FALSE;
  set_initialize(&((*cone)->GroundSet),(*cone)->m_alloc);
  set_initialize(&((*cone)->EqualitySet),(*cone)->m_alloc);
  set_initialize(&((*cone)->NonequalitySet),(*cone)->m_alloc);
  set_initialize(&((*cone)->AddedHalfspaces),(*cone)->m_alloc);
  set_initialize(&((*cone)->WeaklyAddedHalfspaces),(*cone)->m_alloc);
  set_initialize(&((*cone)->InitialHalfspaces),(*cone)->m_alloc);
  (*cone)->RayCount=0;
  (*cone)->FeasibleRayCount=0;
  (*cone)->WeaklyFeasibleRayCount=0;
  (*cone)->TotalRayCount=0;
  (*cone)->ZeroRayCount=0;
  (*cone)->EdgeCount=0;
  (*cone)->TotalEdgeCount=0;
  (*cone)->count_int=0;
  (*cone)->count_int_good=0;
  (*cone)->count_int_bad=0;
  (*cone)->rseed=1;  /* random seed for random row permutation */
 
  dd_InitializeBmatrix((*cone)->d_alloc, &((*cone)->B));
  dd_InitializeBmatrix((*cone)->d_alloc, &((*cone)->Bsave));
  dd_InitializeAmatrix((*cone)->m_alloc,(*cone)->d_alloc,&((*cone)->A));

  (*cone)->Edges
     =(dd_AdjacencyType**) calloc((*cone)->m_alloc,sizeof(dd_AdjacencyType*));
  (*cone)->InitialRayIndex=(long*)calloc(d+1,sizeof(long));
  (*cone)->OrderVector=(long*)calloc((*cone)->m_alloc+1,sizeof(long));

  (*cone)->newcol=(long*)calloc(((*cone)->d)+1,sizeof(long));
  for (j=0; j<=(*cone)->d; j++) (*cone)->newcol[j]=j;  /* identity map, initially */
  (*cone)->LinearityDim = -2; /* -2 if it is not computed */
  (*cone)->ColReduced   = dd_FALSE;
  (*cone)->d_orig = d;

/* STATES: variables to represent current state. */
/*(*cone)->Error;
  (*cone)->CompStatus;
  (*cone)->starttime;
  (*cone)->endtime;
*/
    
  return success;
}

dd_ConePtr dd_ConeDataLoad(dd_PolyhedraPtr poly)
{
  dd_ConePtr cone=NULL;
  dd_colrange d,j;
  dd_rowrange m,i;

  m=poly->m;
  d=poly->d;
  if (!(poly->homogeneous) && poly->representation==dd_Inequality){
    m=poly->m+1;
  }
  poly->m1=m;

  dd_InitializeConeData(m, d, &cone);
  cone->representation=poly->representation;

/* Points to the original polyhedra data, and reversely */
  cone->parent=poly;
  poly->child=cone;

  for (i=1; i<=poly->m; i++)
    for (j=1; j<=cone->d; j++)
      dd_set(cone->A[i-1][j-1],poly->A[i-1][j-1]);  
  
  if (poly->representation==dd_Inequality && !(poly->homogeneous)){
    dd_set(cone->A[m-1][0],dd_one);
    for (j=2; j<=d; j++) dd_set(cone->A[m-1][j-1],dd_purezero);
  }

  return cone;
}

void dd_SetLinearity(dd_MatrixPtr M, char *line)
{
  int i=0;
  dd_rowrange eqsize,var;
  char *next;
  const char ct[]=", ";  /* allows separators "," and " ". */

  next=strtok(line,ct);
  eqsize=atol(next); 
  while (i < eqsize && (next=strtok(NULL,ct))!=NULL) {
     var=atol(next);
     set_addelem(M->linset,var); i++;
  }
  if (i!=eqsize) {
    fprintf(stderr,"* Warning: there are inconsistencies in linearity setting.\n");
  }
  return;
}

dd_MatrixPtr dd_PolyFile2Matrix (FILE *f, dd_ErrorType *Error)
{
  dd_MatrixPtr M=NULL;
  dd_rowrange m_input,i;
  dd_colrange d_input,j;
  dd_RepresentationType rep=dd_Inequality;
  mytype value;
  dd_boolean found=dd_FALSE, newformat=dd_FALSE, successful=dd_FALSE, linearity=dd_FALSE;
  char command[dd_linelenmax], comsave[dd_linelenmax], numbtype[dd_wordlenmax];
  dd_NumberType NT;
#if !defined(GMPRATIONAL)
  double rvalue;
#endif

  dd_init(value);
  (*Error)=dd_NoError;
  while (!found)
  {
    if (fscanf(f,"%s",command)==EOF) {
      (*Error)=dd_ImproperInputFormat;
      goto _L99;
    }
    else {
      if (strncmp(command, "V-representation", 16)==0) {
        rep=dd_Generator; newformat=dd_TRUE;
      }
      if (strncmp(command, "H-representation", 16)==0){
        rep=dd_Inequality; newformat=dd_TRUE;
      }
      if (strncmp(command, "partial_enum", 12)==0 || 
          strncmp(command, "equality", 8)==0  ||
          strncmp(command, "linearity", 9)==0 ) {
        linearity=dd_TRUE;
        fgets(comsave,dd_linelenmax,f);
      }
      if (strncmp(command, "begin", 5)==0) found=dd_TRUE;
    }
  }
  fscanf(f, "%ld %ld %s", &m_input, &d_input, numbtype);
  fprintf(stderr,"size = %ld x %ld\nNumber Type = %s\n", m_input, d_input, numbtype);
  NT=dd_GetNumberType(numbtype);
  if (NT==dd_Unknown) {
      (*Error)=dd_ImproperInputFormat;
      goto _L99;
    } 
  M=dd_CreateMatrix(m_input, d_input);
  M->representation=rep;
  M->numbtype=NT;

  for (i = 1; i <= m_input; i++) {
    for (j = 1; j <= d_input; j++) {
      if (NT==dd_Real) {
#if defined GMPRATIONAL
        *Error=dd_NoRealNumberSupport;
        goto _L99;
#else
        fscanf(f, "%lf", &rvalue);
        dd_set_d(value, rvalue);
#endif
      } else {
        dd_fread_rational_value (f, value);
      }
      dd_set(M->matrix[i-1][j - 1],value);
      if (dd_debug) {fprintf(stderr,"a(%3ld,%5ld) = ",i,j); dd_WriteNumber(stderr,value);}
    }  /*of j*/
  }  /*of i*/
  if (fscanf(f,"%s",command)==EOF) {
   	 (*Error)=dd_ImproperInputFormat;
  	 goto _L99;
  }
  else if (strncmp(command, "end", 3)!=0) {
     if (dd_debug) fprintf(stderr,"'end' missing or illegal extra data: %s\n",command);
     (*Error)=dd_ImproperInputFormat;
     goto _L99;
  }
  
  successful=dd_TRUE;
  if (linearity) {
    dd_SetLinearity(M,comsave);
  }
  while (!feof(f)) {
    fscanf(f,"%s", command);
    dd_ProcessCommandLine(f, M, command);
    fgets(command,dd_linelenmax,f); /* skip the CR/LF */
  } 

_L99: ;
  dd_clear(value);
  /* if (f!=NULL) fclose(f); */
  return M;
}


dd_PolyhedraPtr dd_DDMatrix2Poly(dd_MatrixPtr M, dd_ErrorType *err)
{
  dd_rowrange i;
  dd_colrange j;
  dd_PolyhedraPtr poly=NULL;

  *err=dd_NoError;
  if (M->rowsize<0 || M->colsize<0){
    *err=dd_NegativeMatrixSize;
    goto _L99;
  }
  poly=dd_CreatePolyhedraData(M->rowsize, M->colsize);
  poly->representation=M->representation;
  poly->homogeneous=dd_TRUE;

  for (i = 1; i <= M->rowsize; i++) {
    if (set_member(i, M->linset)) {
      poly->EqualityIndex[i]=1;
    }
    for (j = 1; j <= M->colsize; j++) {
      dd_set(poly->A[i-1][j-1], M->matrix[i-1][j-1]);
      if (j==1 && dd_Nonzero(M->matrix[i-1][j-1])) poly->homogeneous = dd_FALSE;
    }  /*of j*/
  }  /*of i*/
  dd_DoubleDescription(poly,err);
_L99:
  return poly; 
}

dd_PolyhedraPtr dd_DDMatrix2Poly2(dd_MatrixPtr M, dd_RowOrderType horder, dd_ErrorType *err)
{
  dd_rowrange i;
  dd_colrange j;
  dd_PolyhedraPtr poly=NULL;

  *err=dd_NoError;
  if (M->rowsize<0 || M->colsize<0){
    *err=dd_NegativeMatrixSize;
    goto _L99;
  }
  poly=dd_CreatePolyhedraData(M->rowsize, M->colsize);
  poly->representation=M->representation;
  poly->homogeneous=dd_TRUE;

  for (i = 1; i <= M->rowsize; i++) {