cddcore.c

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

/* cddcore.c:  Core 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_CheckAdjacency(dd_ConePtr cone,
    dd_RayPtr *RP1, dd_RayPtr *RP2, dd_boolean *adjacent)
{
  dd_RayPtr TempRay;
  dd_boolean localdebug=dd_FALSE;
  static dd_rowset Face, Face1;
  static dd_rowrange last_m=0;
  
  if (last_m!=cone->m) {
    if (last_m>0){
      set_free(Face); set_free(Face1);
    }
    set_initialize(&Face, cone->m); 
    set_initialize(&Face1, cone->m); 
    last_m=cone->m;
  }

  if (dd_debug) localdebug=dd_TRUE;
  *adjacent = dd_TRUE;
  set_int(Face1, (*RP1)->ZeroSet, (*RP2)->ZeroSet);
  set_int(Face, Face1, cone->AddedHalfspaces);
  if (set_card(Face)< cone->d - 2) {
    *adjacent = dd_FALSE;
    if (localdebug) {
      fprintf(stderr,"non adjacent: set_card(face) %ld < %ld = cone->d.\n",
        set_card(Face),cone->d);
    }
    return;
  }
  else if (cone->parent->NondegAssumed) {
  	*adjacent = dd_TRUE;
  	return;
  }
  TempRay = cone->FirstRay;
  while (TempRay != NULL && *adjacent) {
    if (TempRay != *RP1 && TempRay != *RP2) {
    	set_int(Face1, TempRay->ZeroSet, cone->AddedHalfspaces);
      	if (set_subset(Face, Face1)) *adjacent = dd_FALSE;
    }
    TempRay = TempRay->Next;
  }
}

void dd_Eliminate(dd_ConePtr cone, dd_RayPtr*Ptr)
{
  /*eliminate the record pointed by Ptr^.Next*/
  dd_RayPtr TempPtr;
  dd_colrange j;

  TempPtr = (*Ptr)->Next;
  (*Ptr)->Next = (*Ptr)->Next->Next;
  if (TempPtr == cone->FirstRay)   /*Update the first pointer*/
    cone->FirstRay = (*Ptr)->Next;
  if (TempPtr == cone->LastRay)   /*Update the last pointer*/
    cone->LastRay = *Ptr;

  /* Added, Marc Pfetsch 010219 */
  for (j=0;j < cone->d;j++)
dd_clear(TempPtr->Ray[j]);
  dd_clear(TempPtr->ARay);

  free(TempPtr->Ray);          /* free the ray vector memory */
  set_free(TempPtr->ZeroSet);  /* free the ZeroSet memory */
  free(TempPtr);   /* free the dd_Ray structure memory */
  cone->RayCount--; 
}

void dd_SetInequalitySets(dd_ConePtr cone)
{
  dd_rowrange i;
  
  set_emptyset(cone->GroundSet);
  set_emptyset(cone->EqualitySet);
  set_emptyset(cone->NonequalitySet);  
  for (i = 1; i <= (cone->parent->m); i++){
    set_addelem(cone->GroundSet, i);
    if (cone->parent->EqualityIndex[i]==1) set_addelem(cone->EqualitySet,i);
    if (cone->parent->EqualityIndex[i]==-1) set_addelem(cone->NonequalitySet,i);
  }
}


void dd_AValue(mytype *val, dd_colrange d_size, dd_Amatrix A, mytype *p, dd_rowrange i)
{
  /*return the ith component of the vector  A x p */
  dd_colrange j;
  mytype x;

  dd_init(x); 
  dd_set(*val,dd_purezero);
 /* Changed by Marc Pfetsch 010219 */

  for (j = 0; j < d_size; j++){
    dd_mul(x,A[i - 1][j], p[j]);
    dd_add(*val, *val, x);
  }
  dd_clear(x);
}

void dd_StoreRay1(dd_ConePtr cone, mytype *p, dd_boolean *feasible)
{  /* Original ray storing routine when RelaxedEnumeration is dd_FALSE */
  dd_rowrange i,k,fii=cone->m+1;
  dd_colrange j;
  mytype temp;
  dd_RayPtr RR;
  dd_boolean localdebug=dd_debug;

  dd_init(temp);
  RR=cone->LastRay;
  *feasible = dd_TRUE;
  set_initialize(&(RR->ZeroSet),cone->m);
  for (j = 0; j < cone->d; j++){
    dd_set(RR->Ray[j],p[j]);
  }
  for (i = 1; i <= cone->m; i++) {
    k=cone->OrderVector[i];
    dd_AValue(&temp, cone->d, cone->A, p, k);
    if (localdebug) {
      fprintf(stderr,"dd_StoreRay1: dd_AValue at row %ld =",k);
      dd_WriteNumber(stderr, temp);
      fprintf(stderr,"\n");
    }
    if (dd_EqualToZero(temp)) {
      set_addelem(RR->ZeroSet, k);
      if (localdebug) {
        fprintf(stderr,"recognized zero!\n");
      }
    }
    if (dd_Negative(temp)){
      if (localdebug) {
        fprintf(stderr,"recognized negative!\n");
      }
      *feasible = dd_FALSE;
      if (fii>cone->m) fii=i;  /* the first violating inequality index */
      if (localdebug) {
        fprintf(stderr,"this ray is not feasible, neg comp = %ld\n", fii);
        dd_WriteNumber(stderr, temp);  fprintf(stderr,"\n");
      }
    }
  }
  RR->FirstInfeasIndex=fii;
  RR->feasible = *feasible;
  dd_clear(temp);
}

void dd_StoreRay2(dd_ConePtr cone, mytype *p, 
    dd_boolean *feasible, dd_boolean *weaklyfeasible)
   /* Ray storing routine when RelaxedEnumeration is dd_TRUE.
       weaklyfeasible is true iff it is feasible with
       the strict_inequality conditions deleted. */
{
  dd_RayPtr RR;
  dd_rowrange i,k,fii=cone->m+1;
  dd_colrange j;
  mytype temp;
  dd_boolean localdebug=dd_debug;

  dd_init(temp);
  RR=cone->LastRay;
  if (dd_debug) localdebug=dd_TRUE;
  *feasible = dd_TRUE;
  *weaklyfeasible = dd_TRUE;
  set_initialize(&(RR->ZeroSet),cone->m);
  for (j = 0; j < cone->d; j++){
    dd_set(RR->Ray[j],p[j]);
  }
  for (i = 1; i <= cone->m; i++) {
    k=cone->OrderVector[i];
    dd_AValue(&temp, cone->d, cone->A, p, k);
    if (dd_EqualToZero(temp)){
      set_addelem(RR->ZeroSet, k);
      if (cone->parent->EqualityIndex[k]==-1) 
        *feasible=dd_FALSE;  /* strict inequality required */
    }
/*    if (temp < -zero){ */
    if (dd_Negative(temp)){
      *feasible = dd_FALSE;
      if (fii>cone->m && cone->parent->EqualityIndex[k]>=0) {
        fii=i;  /* the first violating inequality index */
        *weaklyfeasible=dd_FALSE;
      }
    }
  }
  RR->FirstInfeasIndex=fii;
  RR->feasible = *feasible;
  dd_clear(temp);
}


void dd_AddRay(dd_ConePtr cone, mytype *p)
{  
  dd_boolean feasible, weaklyfeasible;
  dd_colrange j;

  if (cone->FirstRay == NULL) {
    cone->FirstRay = (dd_RayPtr) malloc(sizeof(dd_RayType));
    cone->FirstRay->Ray = (mytype *) calloc(cone->d, sizeof(mytype));
    for (j=0; j<cone->d; j++) dd_init(cone->FirstRay->Ray[j]);
    dd_init(cone->FirstRay->ARay);
    if (dd_debug)
      fprintf(stderr,"Create the first ray pointer\n");
    cone->LastRay = cone->FirstRay;
    cone->ArtificialRay->Next = cone->FirstRay;
  } else {
    cone->LastRay->Next = (dd_RayPtr) malloc(sizeof(dd_RayType));
    cone->LastRay->Next->Ray = (mytype *) calloc(cone->d, sizeof(mytype));
    for (j=0; j<cone->d; j++) dd_init(cone->LastRay->Next->Ray[j]);
    dd_init(cone->LastRay->Next->ARay);
    if (dd_debug) fprintf(stderr,"Create a new ray pointer\n");
    cone->LastRay = cone->LastRay->Next;
  }
  cone->LastRay->Next = NULL;
  cone->RayCount++;
  cone->TotalRayCount++;
  if (dd_debug) {
    if (cone->TotalRayCount % 100 == 0) {
      fprintf(stderr,"*Rays (Total, Currently Active, Feasible) =%8ld%8ld%8ld\n",
	 cone->TotalRayCount, cone->RayCount, cone->FeasibleRayCount);
    }
  }
  if (cone->parent->RelaxedEnumeration){
    dd_StoreRay2(cone, p, &feasible, &weaklyfeasible);
    if (weaklyfeasible) (cone->WeaklyFeasibleRayCount)++;
  } else {
    dd_StoreRay1(cone, p, &feasible);
    if (feasible) (cone->WeaklyFeasibleRayCount)++;
    /* weaklyfeasible is equiv. to feasible in this case. */
  }
  if (!feasible) return;
  else {
    (cone->FeasibleRayCount)++;
  }
}

void dd_AddArtificialRay(dd_ConePtr cone)
{  
  dd_Arow zerovector;
  dd_colrange j,d1;
  dd_boolean feasible;

  if (cone->d<=0) d1=1; else d1=cone->d;
  dd_InitializeArow(d1, &zerovector);
  if (cone->ArtificialRay != NULL) {
    fprintf(stderr,"Warning !!!  FirstRay in not nil.  Illegal Call\n");
    free(zerovector); /* 086 */
    return;
  }
  cone->ArtificialRay = (dd_RayPtr) malloc(sizeof(dd_RayType));
  cone->ArtificialRay->Ray = (mytype *) calloc(d1, sizeof(mytype));
  for (j=0; j<d1; j++) dd_init(cone->ArtificialRay->Ray[j]);
  dd_init(cone->ArtificialRay->ARay);

  if (dd_debug) fprintf(stderr,"Create the artificial ray pointer\n");

  cone->LastRay=cone->ArtificialRay;
  dd_StoreRay1(cone, zerovector, &feasible);  
    /* This stores a vector to the record pointed by cone->LastRay */
  cone->ArtificialRay->Next = NULL;
  for (j = 0; j < d1; j++){
    dd_clear(zerovector[j]);
  }
  free(zerovector); /* 086 */
}

void dd_ConditionalAddEdge(dd_ConePtr cone, 
    dd_RayPtr Ray1, dd_RayPtr Ray2, dd_RayPtr ValidFirstRay)
{
  long it,it_row,fii1,fii2,fmin,fmax;
  dd_boolean adjacent,lastchance;
  dd_RayPtr TempRay,Rmin,Rmax;
  dd_AdjacencyType *NewEdge;
  dd_boolean localdebug=dd_FALSE;
  dd_rowset ZSmin, ZSmax;
  static dd_rowset Face, Face1;
  static dd_rowrange last_m=0;
  
  if (last_m!=cone->m) {
    if (last_m>0){
      set_free(Face); set_free(Face1);
    }
    set_initialize(&Face, cone->m);
    set_initialize(&Face1, cone->m);
    last_m=cone->m;
  }
  
  fii1=Ray1->FirstInfeasIndex;
  fii2=Ray2->FirstInfeasIndex;
  if (fii1<fii2){
    fmin=fii1; fmax=fii2;
    Rmin=Ray1;
    Rmax=Ray2;
  }
  else{
    fmin=fii2; fmax=fii1;
    Rmin=Ray2;
    Rmax=Ray1;
  }
  ZSmin = Rmin->ZeroSet;
  ZSmax = Rmax->ZeroSet;
  if (localdebug) {
    fprintf(stderr,"dd_ConditionalAddEdge: FMIN = %ld (row%ld)   FMAX=%ld\n",
      fmin, cone->OrderVector[fmin], fmax);
  }
  if (fmin==fmax){
    if (localdebug) fprintf(stderr,"dd_ConditionalAddEdge: equal FII value-> No edge added\n");
  }
  else if (set_member(cone->OrderVector[fmin],ZSmax)){
    if (localdebug) fprintf(stderr,"dd_ConditionalAddEdge: No strong separation -> No edge added\n");
  }
  else {  /* the pair will be separated at the iteration fmin */
    lastchance=dd_TRUE;
    /* flag to check it will be the last chance to store the edge candidate */
    set_int(Face1, ZSmax, ZSmin);
    (cone->count_int)++;
    if (localdebug){
      fprintf(stderr,"Face: ");
      for (it=1; it<=cone->m; it++) {
        it_row=cone->OrderVector[it];
        if (set_member(it_row, Face1)) fprintf(stderr,"%ld ",it_row);
      }
      fprintf(stderr,"\n");
    }
    for (it=cone->Iteration+1; it < fmin && lastchance; it++){
      it_row=cone->OrderVector[it];
      if (cone->parent->EqualityIndex[it_row]>=0 && set_member(it_row, Face1)){
        lastchance=dd_FALSE;
        (cone->count_int_bad)++;
        if (localdebug){
          fprintf(stderr,"There will be another chance iteration %ld (row %ld) to store the pair\n", it, it_row);
        }
      }
    }
    if (lastchance){
      adjacent = dd_TRUE;
      (cone->count_int_good)++;
      /* adjacent checking */
      set_int(Face, Face1, cone->AddedHalfspaces);
      if (localdebug){
        fprintf(stderr,"Check adjacency\n");
        fprintf(stderr,"AddedHalfspaces: "); set_fwrite(stderr,cone->AddedHalfspaces);
        fprintf(stderr,"Face: ");
        for (it=1; it<=cone->m; it++) {
          it_row=cone->OrderVector[it];
          if (set_member(it_row, Face)) fprintf(stderr,"%ld ",it_row);
        }
        fprintf(stderr,"\n");
      }
      if (set_card(Face)< cone->d - 2) {
        adjacent = dd_FALSE;
      }
      else if (cone->parent->NondegAssumed) {
    	adjacent = dd_TRUE;
      }
      else{
        TempRay = ValidFirstRay;  /* the first ray for adjacency checking */
        while (TempRay != NULL && adjacent) {
          if (TempRay != Ray1 && TempRay != Ray2) {
            set_int(Face1, TempRay->ZeroSet, cone->AddedHalfspaces);
            if (set_subset(Face, Face1)) {
              if (localdebug) set_fwrite(stderr,Face1);
              adjacent = dd_FALSE;
            }
          }
          TempRay = TempRay->Next;
        }
      }
      if (adjacent){
        if (localdebug) fprintf(stderr,"The pair is adjacent and the pair must be stored for iteration %ld (row%ld)\n",
          fmin, cone->OrderVector[fmin]);
        NewEdge=(dd_AdjacencyPtr) malloc(sizeof *NewEdge);
        NewEdge->Ray1=Rmax;  /* save the one remains in iteration fmin in the first */
        NewEdge->Ray2=Rmin;  /* save the one deleted in iteration fmin in the second */
        NewEdge->Next=NULL;
        (cone->EdgeCount)++; 
        (cone->TotalEdgeCount)++;
        if (cone->Edges[fmin]==NULL){
          cone->Edges[fmin]=NewEdge;
          if (localdebug) fprintf(stderr,"Create a new edge list of %ld\n", fmin);
        }else{
          NewEdge->Next=cone->Edges[fmin];
          cone->Edges[fmin]=NewEdge;
        }
      }
    }
  }
}

void dd_CreateInitialEdges(dd_ConePtr cone)
{
  dd_RayPtr Ptr1, Ptr2;
  dd_rowrange fii1,fii2;
  long count=0;
  dd_boolean adj,localdebug=dd_FALSE;

  cone->Iteration=cone->d;  /* CHECK */
  if (cone->FirstRay ==NULL || cone->LastRay==NULL){
    /* fprintf(stderr,"Warning: dd_ CreateInitialEdges called with NULL pointer(s)\n"); */
    goto _L99;
  }
  Ptr1=cone->FirstRay;
  while(Ptr1!=cone->LastRay && Ptr1!=NULL){
    fii1=Ptr1->FirstInfeasIndex;
    Ptr2=Ptr1->Next;
    while(Ptr2!=NULL){
      fii2=Ptr2->FirstInfeasIndex;
      count++;
      if (localdebug) fprintf(stderr,"dd_ CreateInitialEdges: edge %ld \n",count);
      dd_CheckAdjacency(cone, &Ptr1, &Ptr2, &adj);
      if (fii1!=fii2 && adj) 
        dd_ConditionalAddEdge(cone, Ptr1, Ptr2, cone->FirstRay);
      Ptr2=Ptr2->Next;
    }
    Ptr1=Ptr1->Next;
  }
_L99:;  
}


void dd_UpdateEdges(dd_ConePtr cone, dd_RayPtr RRbegin, dd_RayPtr RRend)
/* This procedure must be called after the ray list is sorted
   by dd_EvaluateARay2 so that FirstInfeasIndex's are monotonically
   increasing.
*/
{
  dd_RayPtr Ptr1, Ptr2begin, Ptr2;
  dd_rowrange fii1;
  dd_boolean ptr2found,quit,localdebug=dd_FALSE;
  long count=0,pos1, pos2;
  float workleft,prevworkleft=110.0,totalpairs;

  totalpairs=(cone->ZeroRayCount-1.0)*(cone->ZeroRayCount-2.0)+1.0;
  Ptr2begin = NULL; 
  if (RRbegin ==NULL || RRend==NULL){
    if (1) fprintf(stderr,"Warning: dd_UpdateEdges called with NULL pointer(s)\n");
    goto _L99;
  }
  Ptr1=RRbegin;
  pos1=1;
  do{
    ptr2found=dd_FALSE;
    quit=dd_FALSE;
    fii1=Ptr1->FirstInfeasIndex;
    pos2=2;
    for (Ptr2=Ptr1->Next; !ptr2found && !quit; Ptr2=Ptr2->Next,pos2++){
      if  (Ptr2->FirstInfeasIndex > fii1){
        Ptr2begin=Ptr2;
        ptr2found=dd_TRUE;
      }
      else if (Ptr2==RRend) quit=dd_TRUE;
    }
    if (ptr2found){
      quit=dd_FALSE;
      for (Ptr2=Ptr2begin; !quit ; Ptr2=Ptr2->Next){
        count++;
        if (localdebug) fprintf(stderr,"dd_UpdateEdges: edge %ld \n",count);
        dd_ConditionalAddEdge(cone, Ptr1,Ptr2,RRbegin);
        if (Ptr2==RRend || Ptr2->Next==NULL) quit=dd_TRUE;
      }
    }
    Ptr1=Ptr1->Next;
    pos1++;
    workleft = 100.0 * (cone->ZeroRayCount-pos1) * (cone->ZeroRayCount - pos1-1.0) / totalpairs;
    if (cone->ZeroRayCount>=500 && dd_debug && pos1%10 ==0 && prevworkleft-workleft>=10 ) {
      fprintf(stderr,"*Work of iteration %5ld(/%ld): %4ld/%4ld => %4.1f%% left\n",
	     cone->Iteration, cone->m, pos1, cone->ZeroRayCount, workleft);
      prevworkleft=workleft;
    }    
  }while(Ptr1!=RRend && Ptr1!=NULL);
_L99:;  
}

void dd_FreeDDMemory0(dd_ConePtr cone)
{
  dd_RayPtr Ptr, PrevPtr;
  long count;
  dd_colrange j;
  dd_boolean localdebug=dd_FALSE;
  
  /* THIS SHOULD BE REWRITTEN carefully */
  PrevPtr=cone->ArtificialRay;
  if (PrevPtr!=NULL){
    count=0;
    for (Ptr=cone->ArtificialRay->Next; Ptr!=NULL; Ptr=Ptr->Next){
      /* Added Marc Pfetsch 2/19/01 */
      for (j=0;j < cone->d;j++)
dd_clear(PrevPtr->Ray[j]);
      dd_clear(PrevPtr->ARay);

      free(PrevPtr->Ray);
      free(PrevPtr->ZeroSet);
      free(PrevPtr);
      count++;
      PrevPtr=Ptr;
    };
    cone->FirstRay=NULL;
    /* Added Marc Pfetsch 010219 */
    for (j=0;j < cone->d;j++)
dd_clear(cone->LastRay->Ray[j]);
    dd_clear(cone->LastRay->ARay);

    free(cone->LastRay->Ray);
    cone->LastRay->Ray = NULL;
    set_free(cone->LastRay->ZeroSet);
    cone->LastRay->ZeroSet = NULL;
    free(cone->LastRay);
    cone->LastRay = NULL;
    cone->ArtificialRay=NULL;
    if (localdebug) fprintf(stderr,"%ld ray storage spaces freed\n",count);