cddio.c

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

void dd_WriteLPTimes(FILE *f, dd_LPPtr lp)
{
  dd_WriteTimes(f,lp->starttime,lp->endtime);
}

void dd_WriteLPStats(FILE *f)
{
  time_t currenttime;
  
  time(&currenttime);
  
  fprintf(f, "\n*--- Statistics of pivots ---\n");
#if defined GMPRATIONAL
  fprintf(f, "* f0 = %ld (float basis finding pivots)\n",ddf_statBApivots);
  fprintf(f, "* fc = %ld (float CC pivots)\n",ddf_statCCpivots);
  fprintf(f, "* f1 = %ld (float dual simplex phase I pivots)\n",ddf_statDS1pivots);
  fprintf(f, "* f2 = %ld (float dual simplex phase II pivots)\n",ddf_statDS2pivots);
  fprintf(f, "* f3 = %ld (float anticycling CC pivots)\n",ddf_statACpivots);
  fprintf(f, "* e0 = %ld (exact basis finding pivots)\n",dd_statBApivots);
  fprintf(f, "* ec = %ld (exact CC pivots)\n",dd_statCCpivots);
  fprintf(f, "* e1 = %ld (exact dual simplex phase I pivots)\n",dd_statDS1pivots);
  fprintf(f, "* e2 = %ld (exact dual simplex phase II pivots)\n",dd_statDS2pivots);
  fprintf(f, "* e3 = %ld (exact anticycling CC pivots)\n",dd_statACpivots);
  fprintf(f, "* e4 = %ld (exact basis verification pivots)\n",dd_statBSpivots);
#else
  fprintf(f, "f0 = %ld (float basis finding pivots)\n",dd_statBApivots);
  fprintf(f, "fc = %ld (float CC pivots)\n",dd_statCCpivots);
  fprintf(f, "f1 = %ld (float dual simplex phase I pivots)\n",dd_statDS1pivots);
  fprintf(f, "f2 = %ld (float dual simplex phase II pivots)\n",dd_statDS2pivots);
  fprintf(f, "f3 = %ld (float anticycling CC pivots)\n",dd_statACpivots);
#endif
 dd_WriteLPMode(f);
  dd_WriteTimes(f,dd_statStartTime, currenttime);
}

void dd_WriteLPMode(FILE *f)
{
  fprintf(f, "\n* LP solver: ");
  switch (dd_choiceLPSolverDefault) {
    case dd_DualSimplex:
      fprintf(f, "DualSimplex\n");
      break;
    case dd_CrissCross:
      fprintf(f, "Criss-Cross\n");
      break;
    default: break;
  }
  
  fprintf(f, "* Redundancy cheking solver: ");
  switch (dd_choiceRedcheckAlgorithm) {
    case dd_DualSimplex:
      fprintf(f, "DualSimplex\n");
      break;
    case dd_CrissCross:
      fprintf(f, "Criss-Cross\n");
      break;
    default: break;
  }
  
  fprintf(f, "* Lexicographic pivot: ");
  if (dd_choiceLexicoPivotQ)  fprintf(f, " on\n"); 
  else fprintf(f, " off\n"); 

}


void dd_WriteRunningMode(FILE *f, dd_PolyhedraPtr poly)
{
  if (poly->child!=NULL){
    fprintf(f,"* roworder: ");
    switch (poly->child->HalfspaceOrder) {

    case dd_MinIndex:
      fprintf(f, "minindex\n");
      break;

    case dd_MaxIndex:
      fprintf(f, "maxindex\n");
      break;

    case dd_MinCutoff:
      fprintf(f, "mincutoff\n");
      break;

    case dd_MaxCutoff:
      fprintf(f, "maxcutoff\n");
    break;

    case dd_MixCutoff:
      fprintf(f, "mixcutoff\n");
      break;

    case dd_LexMin:
      fprintf(f, "lexmin\n");
      break;

    case dd_LexMax:
      fprintf(f, "lexmax\n");
      break;

    case dd_RandomRow:
      fprintf(f, "random  %d\n",poly->child->rseed);
      break;

    default: break;
    }
  }
}


void dd_WriteCompletionStatus(FILE *f, dd_ConePtr cone)
{
  if (cone->Iteration<cone->m && cone->CompStatus==dd_AllFound) {
    fprintf(f,"*Computation completed at Iteration %4ld.\n", cone->Iteration);
  } 
  if (cone->CompStatus == dd_RegionEmpty) {
    fprintf(f,"*Computation completed at Iteration %4ld because the region found empty.\n",cone->Iteration);
  }   
}

void dd_WritePolyFile(FILE *f, dd_PolyhedraPtr poly)
{
  dd_WriteAmatrix(f,poly->A,poly->m,poly->d);
}


void dd_WriteErrorMessages(FILE *f, dd_ErrorType Error)
{
  switch (Error) {

  case dd_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 dd_IFileNotFound:
    fprintf(f, "*Input Error: Specified input file does not exist.\n");
    break;

  case dd_OFileNotOpen:
    fprintf(f, "*Output Error: Specified output file cannot be opened.\n");
    break;

  case dd_NegativeMatrixSize:
    fprintf(f, "*Input Error: Input matrix has a negative size:\n");
    fprintf(f, "*Please check rowsize or colsize.\n");
    break;

  case dd_ImproperInputFormat:
    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 dd_EmptyVrepresentation:
    fprintf(f, "*Input Error: V-representation is empty:\n");
    fprintf(f, "*cddlib does not accept this trivial case for which output can be any inconsistent system.\n");
    break;

  case dd_EmptyHrepresentation:
    fprintf(f, "*Input Error: H-representation is empty.\n");
    break;

  case dd_EmptyRepresentation:
    fprintf(f, "*Input Error: Representation is empty.\n");
    break;

  case dd_NoLPObjective:
    fprintf(f, "*LP Error: No LP objective (max or min) is set.\n");
    break;

  case dd_NoRealNumberSupport:
    fprintf(f, "*LP Error: The binary (with GMP Rational) does not support Real number input.\n");
    fprintf(f, "         : Use a binary compiled without -DGMPRATIONAL option.\n");
    break;

 case dd_NotAvailForH:
    fprintf(f, "*Error: A function is called with H-rep which does not support an H-representation.\n");
    break;

 case dd_NotAvailForV:
    fprintf(f, "*Error: A function is called with V-rep which does not support an V-representation.\n");
    break;

 case dd_CannotHandleLinearity:
    fprintf(f, "*Error: The function called cannot handle linearity.\n");
    break;

 case dd_RowIndexOutOfRange:
    fprintf(f, "*Error: Specified row index is out of range\n");
    break;

 case dd_ColIndexOutOfRange:
    fprintf(f, "*Error: Specified column index is out of range\n");
    break;

 case dd_LPCycling:
    fprintf(f, "*Error: Possibly an LP cycling occurs.  Use the Criss-Cross method.\n");
    break;
    
 case dd_NumericallyInconsistent:
    fprintf(f, "*Error: Numerical inconsistency is found.  Use the GMP exact arithmetic.\n");
    break;
    
  case dd_NoError:
    fprintf(f,"*No Error found.\n");
    break;
  }
}


dd_SetFamilyPtr dd_CopyIncidence(dd_PolyhedraPtr poly)
{
  dd_SetFamilyPtr F=NULL;
  dd_bigrange k;
  dd_rowrange i;

  if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
  if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
  F=dd_CreateSetFamily(poly->n, poly->m1);
  for (i=1; i<=poly->m1; i++)
    for (k=1; k<=poly->n; k++)
      if (set_member(k,poly->Ainc[i-1])) set_addelem(F->set[k-1],i);
_L99:;
  return F;
}

dd_SetFamilyPtr dd_CopyInputIncidence(dd_PolyhedraPtr poly)
{
  dd_rowrange i;
  dd_SetFamilyPtr F=NULL;

  if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
  if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
  F=dd_CreateSetFamily(poly->m1, poly->n);
  for(i=0; i< poly->m1; i++){
    set_copy(F->set[i], poly->Ainc[i]);
  }
_L99:;
  return F;
}

dd_SetFamilyPtr dd_CopyAdjacency(dd_PolyhedraPtr poly)
{
  dd_RayPtr RayPtr1,RayPtr2;
  dd_SetFamilyPtr F=NULL;
  long pos1, pos2;
  dd_bigrange lstart,k,n;
  set_type linset,allset;
  dd_boolean adj;

  if (poly->n==0 && poly->homogeneous && poly->representation==dd_Inequality){
    n=1; /* the origin (the unique vertex) should be output. */
  } else n=poly->n;
  set_initialize(&linset, n);
  set_initialize(&allset, n);
  if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
  F=dd_CreateSetFamily(n, n);
  if (n<=0) goto _L99;
  poly->child->LastRay->Next=NULL;
  for (RayPtr1=poly->child->FirstRay, pos1=1;RayPtr1 != NULL; 
				RayPtr1 = RayPtr1->Next, pos1++){
    for (RayPtr2=poly->child->FirstRay, pos2=1; RayPtr2 != NULL; 
					RayPtr2 = RayPtr2->Next, pos2++){
      if (RayPtr1!=RayPtr2){
        dd_CheckAdjacency(poly->child, &RayPtr1, &RayPtr2, &adj);
        if (adj){
          set_addelem(F->set[pos1-1], pos2);
        }
      }
    }
  }
  lstart=poly->n - poly->ldim + 1;
  set_compl(allset,allset);  /* allset is set to the ground set. */
  for (k=lstart; k<=poly->n; k++){
    set_addelem(linset,k);     /* linearity set */
    set_copy(F->set[k-1],allset);  /* linearity generator is adjacent to all */
  }
  for (k=1; k<lstart; k++){
    set_uni(F->set[k-1],F->set[k-1],linset);
     /* every generator is adjacent to all linearity generators */
  }
_L99:;
  set_free(allset); set_free(linset);
  return F;
}

dd_SetFamilyPtr dd_CopyInputAdjacency(dd_PolyhedraPtr poly)
{
  dd_rowrange i,j;
  dd_SetFamilyPtr F=NULL;

  if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;
  if (poly->AincGenerated==dd_FALSE) dd_ComputeAinc(poly);
  F=dd_CreateSetFamily(poly->m1, poly->m1);
  for (i=1; i<=poly->m1; i++){
    for (j=1; j<=poly->m1; j++){
      if (i!=j && dd_InputAdjacentQ(poly, i, j)) {
        set_addelem(F->set[i-1],j);
      }
    }
  }
_L99:;
  return F;
}

dd_MatrixPtr dd_CopyOutput(dd_PolyhedraPtr poly)
{
  dd_RayPtr RayPtr;
  dd_MatrixPtr M=NULL;
  dd_rowrange i=0,total;
  dd_colrange j,j1;
  mytype b;
  dd_RepresentationType outputrep=dd_Inequality;
  dd_boolean outputorigin=dd_FALSE;

  dd_init(b);
  total=poly->child->LinearityDim + poly->child->FeasibleRayCount;
  
  if (poly->child->d<=0 || poly->child->newcol[1]==0) total=total-1;
  if (poly->representation==dd_Inequality) outputrep=dd_Generator;
  if (total==0 && poly->homogeneous && poly->representation==dd_Inequality){
    total=1;
    outputorigin=dd_TRUE;
    /* the origin (the unique vertex) should be output. */
  }
  if (poly->child==NULL || poly->child->CompStatus!=dd_AllFound) goto _L99;

  M=dd_CreateMatrix(total, poly->d);
  RayPtr = poly->child->FirstRay;
  while (RayPtr != NULL) {
    if (RayPtr->feasible) {
     dd_CopyRay(M->matrix[i], poly->d, RayPtr, outputrep, poly->child->newcol);
      i++;  /* 086 */
    }
    RayPtr = RayPtr->Next;
  }
  for (j=2; j<=poly->d; j++){
    if (poly->child->newcol[j]==0){ 
       /* original column j is dependent on others and removed for the cone */
      dd_set(b,poly->child->Bsave[0][j-1]);
      if (outputrep==dd_Generator && dd_Positive(b)){
        dd_set(M->matrix[i][0],dd_one);  /* dd_Normalize */
        for (j1=1; j1<poly->d; j1++) 
          dd_div(M->matrix[i][j1],(poly->child->Bsave[j1][j-1]),b);
      } else {
        for (j1=0; j1<poly->d; j1++)
          dd_set(M->matrix[i][j1],poly->child->Bsave[j1][j-1]);
      }
      set_addelem(M->linset, i+1);
      i++;
    }     
  }
  if (outputorigin){ 
    /* output the origin for homogeneous H-polyhedron with no rays. */
    dd_set(M->matrix[0][0],dd_one);
    for (j=1; j<poly->d; j++){
      dd_set(M->matrix[0][j],dd_purezero);
    }
  }
  dd_MatrixIntegerFilter(M);
  if (poly->representation==dd_Inequality)
    M->representation=dd_Generator;
  else
    M->representation=dd_Inequality;
_L99:;
  dd_clear(b);
  return M;
}

dd_MatrixPtr dd_CopyInput(dd_PolyhedraPtr poly)
{
  dd_MatrixPtr M=NULL;
  dd_rowrange i;

  M=dd_CreateMatrix(poly->m, poly->d);
  dd_CopyAmatrix(M->matrix, poly->A, poly->m, poly->d);
  for (i=1; i<=poly->m; i++) 
    if (poly->EqualityIndex[i]==1) set_addelem(M->linset,i);
  dd_MatrixIntegerFilter(M);
  if (poly->representation==dd_Generator)
    M->representation=dd_Generator;
  else
    M->representation=dd_Inequality;
  return M;
}

dd_MatrixPtr dd_CopyGenerators(dd_PolyhedraPtr poly)
{
  dd_MatrixPtr M=NULL;

  if (poly->representation==dd_Generator){
    M=dd_CopyInput(poly);
  } else {
    M=dd_CopyOutput(poly);
  }
  return M;
}

dd_MatrixPtr dd_CopyInequalities(dd_PolyhedraPtr poly)
{
  dd_MatrixPtr M=NULL;

  if (poly->representation==dd_Inequality){
    M=dd_CopyInput(poly);
  } else {
    M=dd_CopyOutput(poly);
  }
  return M;
}

/****************************************************************************************/
/*  rational number (a/b) read is taken from Vinci by Benno Bueeler and Andreas Enge    */
/****************************************************************************************/
void dd_sread_rational_value (char *s, mytype value)
   /* reads a rational value from the specified string "s" and assigns it to "value"    */
   
{
   char     *numerator_s=NULL, *denominator_s=NULL, *position;
   int      sign = 1;
   double   numerator, denominator;
#if defined GMPRATIONAL
   mpz_t znum, zden;
#else
   double rvalue;
#endif
  
   /* determine the sign of the number */
   numerator_s = s;
   if (s [0] == '-')
   {  sign = -1;
      numerator_s++;
   }
   else if (s [0] == '+')
      numerator_s++;
      
   /* look for a sign '/' and eventually split the number in numerator and denominator */
   position = strchr (numerator_s, '/');
   if (position != NULL)
   {  *position = '\0'; /* terminates the numerator */
      denominator_s = position + 1;
   };

   /* determine the floating point values of numerator and denominator */
   numerator=atol (numerator_s);
  
   if (position != NULL)
   { 
     denominator=atol (denominator_s);  
   }
   else denominator = 1;

/* 
   fprintf(stderr,"\nrational_read: numerator %f\n",numerator);
   fprintf(stderr,"rational_read: denominator %f\n",denominator);
   fprintf(stderr,"rational_read: sign %d\n",sign); 
*/

#if defined GMPRATIONAL
   mpz_init_set_str(znum,numerator_s,10);
   if (sign<0) mpz_neg(znum,znum);
   mpz_init(zden); mpz_set_ui(zden,1);
   if (denominator_s!=NULL) mpz_init_set_str(zden,denominator_s,10);
   mpq_set_num(value,znum); mpq_set_den(value,zden);
   mpq_canonicalize(value);
   mpz_clear(znum); mpz_clear(zden);
   /* num=(long)sign * numerator; */
   /* den=(unsigned long) denominator; */
   /* mpq_set_si(value, num, den); */
#elif defined GMPFLOAT
   rvalue=sign * numerator/ (signed long) denominator;
   mpf_set_d(value, rvalue);
#else
   rvalue=sign * numerator/ (signed long) denominator;
   ddd_set_d(value, rvalue);
#endif
   if (dd_debug) {
     fprintf(stderr,"rational_read: "); 
     dd_WriteNumber(stderr,value); fprintf(stderr,"\n");
   }
}
   

void dd_fread_rational_value (FILE *f, mytype value)
   /* reads a rational value from the specified file "f" and assigns it to "value"      */
   
{
   char     number_s [255];
   mytype rational_value;
   
   dd_init(rational_value);
   fscanf(f, "%s ", number_s);
   dd_sread_rational_value (number_s, rational_value);
   dd_set(value,rational_value);
   dd_clear(rational_value);
}
   
/****************************************************************************************/


/* end of cddio.c */