/* cddcore.c:  Core Procedures for cddlib   written by Komei Fukuda, fukuda@ifor.math.ethz.ch   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);