cddio.c

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

    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_DoubleDescription2(poly, horder, err);
_L99:
  return poly; 
}

void dd_MatrixIntegerFilter(dd_MatrixPtr M)
{   /* setting an almost integer to the integer. */
  dd_rowrange i;
  dd_colrange j;
  mytype x;

  dd_init(x);
  for (i=0; i< M->rowsize; i++)
    for (j=0; j< M->colsize; j++){
       dd_SnapToInteger(x, M->matrix[i][j]);
       dd_set(M->matrix[i][j],x);
    }
  dd_clear(x);
}

void dd_CopyRay(mytype *a, dd_colrange d_origsize, dd_RayPtr RR, 
  dd_RepresentationType rep, dd_colindex reducedcol)
{
  long j,j1;
  mytype b;

  dd_init(b);
  for (j = 1; j <= d_origsize; j++){
    j1=reducedcol[j];
    if (j1>0){
      dd_set(a[j-1],RR->Ray[j1-1]); 
        /* the original column j is mapped to j1, and thus
           copy the corresponding component */
    } else {
      dd_set(a[j-1],dd_purezero);  
        /* original column is redundant and removed for computation */
    }
  }

  dd_set(b,a[0]);
  if (rep==dd_Generator && dd_Nonzero(b)){
    dd_set(a[0],dd_one);
    for (j = 2; j <= d_origsize; j++)
       dd_div(a[j-1],a[j-1],b);    /* normalization for generators */
  }
  dd_clear(b);
}

void dd_WriteRay(FILE *f, dd_colrange d_origsize, dd_RayPtr RR, dd_RepresentationType rep, dd_colindex reducedcol)
{
  dd_colrange j;
  static dd_colrange d_last=0;
  static dd_Arow a;

  if (d_last< d_origsize){
    if (d_last>0) free(a);
    dd_InitializeArow(d_origsize+1, &a);
    d_last=d_origsize+1;
  }

 dd_CopyRay(a, d_origsize, RR, rep, reducedcol);
  for (j = 0; j < d_origsize; j++) dd_WriteNumber(f, a[j]);
  fprintf(f, "\n");
}

void dd_WriteArow(FILE *f, dd_Arow a, dd_colrange d)
{
  dd_colrange j;

  for (j = 0; j < d; j++) dd_WriteNumber(f, a[j]);
  fprintf(f, "\n");
}

void dd_WriteAmatrix(FILE *f, dd_Amatrix A, long rowmax, long colmax)
{
  long i,j;

  if (A==NULL){
    fprintf(f, "WriteAmatrix: The requested matrix is empty\n");
    goto _L99;
  }
  fprintf(f, "begin\n");
#if defined GMPRATIONAL
  fprintf(f, " %ld %ld rational\n",rowmax, colmax);
#else
  fprintf(f, " %ld %ld real\n",rowmax, colmax);
#endif
  for (i=1; i <= rowmax; i++) {
    for (j=1; j <= colmax; j++) {
      dd_WriteNumber(f, A[i-1][j-1]);
    }
    fprintf(f,"\n");
  }
  fprintf(f, "end\n");
_L99:;
}

void dd_WriteBmatrix(FILE *f, dd_colrange d_size, dd_Bmatrix B)
{
  dd_colrange j1, j2;

  if (B==NULL){
    fprintf(f, "WriteBmatrix: The requested matrix is empty\n");
    goto _L99;
  }
  for (j1 = 0; j1 < d_size; j1++) {
    for (j2 = 0; j2 < d_size; j2++) {
      dd_WriteNumber(f, B[j1][j2]);
    }  /*of j2*/
    putc('\n', f);
  }  /*of j1*/
  putc('\n', f);
_L99:;
}

void dd_WriteSetFamily(FILE *f, dd_SetFamilyPtr F)
{
  dd_bigrange i;

  if (F==NULL){
    fprintf(f, "WriteSetFamily: The requested family is empty\n");
    goto _L99;
  }
  fprintf(f,"begin\n");
  fprintf(f,"  %ld    %ld\n", F->famsize, F->setsize);
  for (i=0; i<F->famsize; i++) {
    fprintf(f, " %ld %ld : ", i+1, set_card(F->set[i]));
    set_fwrite(f, F->set[i]);
  }
  fprintf(f,"end\n");
_L99:;
}

void dd_WriteSetFamilyCompressed(FILE *f, dd_SetFamilyPtr F)
{
  dd_bigrange i,card;

  if (F==NULL){
    fprintf(f, "WriteSetFamily: The requested family is empty\n");
    goto _L99;
  }
  fprintf(f,"begin\n");
  fprintf(f,"  %ld    %ld\n", F->famsize, F->setsize);
  for (i=0; i<F->famsize; i++) {
    card=set_card(F->set[i]);
    if (F->setsize - card >= card){
      fprintf(f, " %ld %ld : ", i+1, card);
      set_fwrite(f, F->set[i]);
    } else {
      fprintf(f, " %ld %ld : ", i+1, -card);
      set_fwrite_compl(f, F->set[i]);
    }
  }
  fprintf(f,"end\n");
_L99:;
}

void dd_WriteMatrix(FILE *f, dd_MatrixPtr M)
{
  dd_rowrange i, linsize;

  if (M==NULL){
    fprintf(f, "WriteMatrix: The requested matrix is empty\n");
    goto _L99;
  }
  switch (M->representation) {
    case dd_Inequality:
      fprintf(f, "H-representation\n"); break;
    case dd_Generator:
      fprintf(f, "V-representation\n"); break;
    case dd_Unspecified:
      break;
  }
  linsize=set_card(M->linset);
  if (linsize>0) {
    fprintf(f, "linearity %ld ", linsize);
    for (i=1; i<=M->rowsize; i++) 
      if (set_member(i, M->linset)) fprintf(f, " %ld", i);
    fprintf(f, "\n");
  }
  dd_WriteAmatrix(f, M->matrix, M->rowsize, M->colsize);
  if (M->objective!=dd_LPnone){
    if (M->objective==dd_LPmax)
      fprintf(f, "maximize\n");
    else
      fprintf(f, "minimize\n");      
    dd_WriteArow(f, M->rowvec, M->colsize);
  }
_L99:;
}

void dd_WriteLP(FILE *f, dd_LPPtr lp)
{
  if (lp==NULL){
    fprintf(f, "WriteLP: The requested lp is empty\n");
    goto _L99;
  }
  fprintf(f, "H-representation\n");

  dd_WriteAmatrix(f, lp->A, (lp->m)-1, lp->d);
  if (lp->objective!=dd_LPnone){
    if (lp->objective==dd_LPmax)
      fprintf(f, "maximize\n");
    else
      fprintf(f, "minimize\n");      
    dd_WriteArow(f, lp->A[lp->objrow-1], lp->d);
  }
_L99:;
}


void dd_SnapToInteger(mytype y, mytype x)
{
 /* this is broken.  It does nothing. */
   dd_set(y,x);
}


void dd_WriteReal(FILE *f, mytype x)
{
  long ix1,ix2,ix;
  double ax;

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

void dd_WriteNumber(FILE *f, mytype x)
{
#if defined GMPRATIONAL
  mpz_t zn,zd;

  mpz_init(zn); mpz_init(zd);
  mpq_canonicalize(x);
  mpq_get_num(zn,x);
  mpq_get_den(zd,x);
  fprintf(f," ");
  if (mpz_sgn(zn)==0){
    fprintf(f,"0");
  } else if (mpz_cmp_ui(zd,1U)==0){
    mpz_out_str(f,10,zn);
  } else {
    mpz_out_str(f,10,zn);fprintf(f,"/");mpz_out_str(f,10,zd);
  }
  mpz_clear(zn); mpz_clear(zd);
#else
  dd_WriteReal(f, x);
#endif
}


void dd_WriteIncidence(FILE *f, dd_PolyhedraPtr poly)
{
  dd_SetFamilyPtr I;

  switch (poly->representation) {
    case dd_Inequality:
       fprintf(f, "ecd_file: Incidence of generators and inequalities\n");
      break;
    case dd_Generator:
       fprintf(f, "icd_file: Incidence of inequalities and generators\n");
      break;

    default:
      break;
  }
  I=dd_CopyIncidence(poly);
  dd_WriteSetFamilyCompressed(f,I);
  dd_FreeSetFamily(I);
}

void dd_WriteAdjacency(FILE *f, dd_PolyhedraPtr poly)
{
  dd_SetFamilyPtr A;

  switch (poly->representation) {
    case dd_Inequality:
       fprintf(f, "ead_file: Adjacency of generators\n");
      break;
    case dd_Generator:
       fprintf(f, "iad_file: Adjacency of inequalities\n");
      break;

    default:
      break;
  }
  A=dd_CopyAdjacency(poly);
  dd_WriteSetFamilyCompressed(f,A);
  dd_FreeSetFamily(A);
}


void dd_ComputeAinc(dd_PolyhedraPtr poly)
{
/* This generates the input incidence array poly->Ainc, and
   two sets: poly->Ared, poly->Adom. 
*/
  dd_bigrange k;
  dd_rowrange i,m1;
  dd_colrange j;
  dd_boolean redundant;
  dd_MatrixPtr M=NULL;
  mytype sum,temp;

  dd_init(sum); dd_init(temp);
  if (poly->AincGenerated==dd_TRUE) goto _L99;

  M=dd_CopyOutput(poly);
  poly->n=M->rowsize;
  m1=poly->m1;  
   /* this number is same as poly->m, except when
      poly is given by nonhomogeneous inequalty:
      !(poly->homogeneous) && poly->representation==Inequality,
      it is poly->m+1.   See dd_ConeDataLoad.
   */
  poly->Ainc=(set_type*)calloc(m1, sizeof(set_type));
  for(i=1; i<=m1; i++) set_initialize(&(poly->Ainc[i-1]),poly->n);
  set_initialize(&(poly->Ared), m1); 
  set_initialize(&(poly->Adom), m1); 

  for (k=1; k<=poly->n; k++){
    for (i=1; i<=poly->m; i++){
      dd_set(sum,dd_purezero);
      for (j=1; j<=poly->d; j++){
        dd_mul(temp,poly->A[i-1][j-1],M->matrix[k-1][j-1]);
        dd_add(sum,sum,temp);
      }
      if (dd_EqualToZero(sum)) {
        set_addelem(poly->Ainc[i-1], k);
      }
    }
    if (!(poly->homogeneous) && poly->representation==dd_Inequality){
      if (dd_EqualToZero(M->matrix[k-1][0])) {
        set_addelem(poly->Ainc[m1-1], k);  /* added infinity inequality (1,0,0,...,0) */
      }
    }
  }

  for (i=1; i<=m1; i++){
    if (set_card(poly->Ainc[i-1])==M->rowsize){
      set_addelem(poly->Adom, i);
    }  
  }
  for (i=m1; i>=1; i--){
    if (set_card(poly->Ainc[i-1])==0){
      redundant=dd_TRUE;
      set_addelem(poly->Ared, i);
    }else {
      redundant=dd_FALSE;
      for (k=1; k<=m1; k++) {
        if (k!=i && !set_member(k, poly->Ared)  && !set_member(k, poly->Adom) && 
            set_subset(poly->Ainc[i-1], poly->Ainc[k-1])){
          if (!redundant){
            redundant=dd_TRUE;
          }
          set_addelem(poly->Ared, i);
        }
      }
    }
  }
  dd_FreeMatrix(M);
  poly->AincGenerated=dd_TRUE;
_L99:;
  dd_clear(sum);  dd_clear(temp);
}

dd_boolean dd_InputAdjacentQ(dd_PolyhedraPtr poly, 
  dd_rowrange i1, dd_rowrange i2)
/* Before calling this function, RedundantSet must be 
   a set of row indices whose removal results in a minimal
   nonredundant system to represent the input polyhedron,
   DominantSet must be the set of row indices which are
   active at every extreme points/rays.
*/
{
  dd_boolean adj=dd_TRUE;
  dd_rowrange i;
  static set_type common;
  static long lastn=0;

  if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
  if (lastn!=poly->n){
    if (lastn >0) set_free(common);
    set_initialize(&common, poly->n);
    lastn=poly->n;
  }
  if (set_member(i1, poly->Ared) || set_member(i2, poly->Ared)){
    adj=dd_FALSE;
    goto _L99;
  }
  if (set_member(i1, poly->Adom) || set_member(i2, poly->Adom)){
      /* dominant inequality is considered adjacencent to all others. */
    adj=dd_TRUE;
    goto _L99;
  }
  set_int(common, poly->Ainc[i1-1], poly->Ainc[i2-1]);
  i=0;
  while (i<poly->m1 && adj==dd_TRUE){ 
    i++; 
    if (i!=i1 && i!=i2 && !set_member(i, poly->Ared) &&
        !set_member(i, poly->Adom) && set_subset(common,poly->Ainc[i-1])){
      adj=dd_FALSE;
    }
  }
_L99:;
  return adj;
} 


void dd_WriteInputIncidence(FILE *f, dd_PolyhedraPtr poly)
{
  dd_SetFamilyPtr I;

  if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
  switch (poly->representation) {
  case dd_Inequality:
    fprintf(f,"icd_file: Incidence of inequalities and generators\n");
    break;

  case dd_Generator:
   fprintf(f,"ecd_file: Incidence of generators and inequalities\n");
    break;

  default:
    break;
  }

  I=dd_CopyInputIncidence(poly);
  dd_WriteSetFamilyCompressed(f,I);
  dd_FreeSetFamily(I);
}

void dd_WriteInputAdjacency(FILE *f, dd_PolyhedraPtr poly)
{
  dd_SetFamilyPtr A;

  if (poly->AincGenerated==dd_FALSE){
    dd_ComputeAinc(poly);
  }
  switch (poly->representation) {
  case dd_Inequality:
    fprintf(f, "iad_file: Adjacency of inequalities\n");
    break;

  case dd_Generator:
    fprintf(f, "ead_file: Adjacency of generators\n");
    break;

  default:
    break;
  }
  A=dd_CopyInputAdjacency(poly);
  dd_WriteSetFamilyCompressed(f,A);
  dd_FreeSetFamily(A);
}


void dd_WriteProgramDescription(FILE *f)
{
  fprintf(f, "* cddlib: a double description library:%s\n", dd_DDVERSION);
  fprintf(f, "* compiled for %s arithmetic.\n", dd_ARITHMETIC);
  fprintf(f,"* %s\n",dd_COPYRIGHT);
}

void dd_WriteTimes(FILE *f, time_t starttime, time_t endtime)
{
  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 started at %s",asctime(localtime(&starttime)));
  fprintf(f,"*             ended   at %s",asctime(localtime(&endtime)));
  fprintf(f,"* Total processor time = %ld seconds\n",ptime);
  fprintf(f,"*                      = %ld h %ld m %ld s\n",ptime_hour, ptime_minu, ptime_sec);
}

void dd_WriteDDTimes(FILE *f, dd_PolyhedraPtr poly)
{
  dd_WriteTimes(f,poly->child->starttime,poly->child->endtime);
}