cddcore.c

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

  }
  cone->Edges[cone->Iteration]=NULL;
  
  dd_DeleteNegativeRays(cone);
    
  set_addelem(cone->AddedHalfspaces, hnew);

  if (cone->Iteration<cone->m){
    if (cone->ZeroHead!=NULL && cone->ZeroHead!=cone->LastRay){
      if (cone->ZeroRayCount>200 && dd_debug) fprintf(stderr,"*New edges being scanned...\n");
      dd_UpdateEdges(cone, cone->ZeroHead, cone->LastRay);
    }
  }

  if (cone->RayCount==cone->WeaklyFeasibleRayCount) cone->CompStatus=dd_AllFound;
_L99:;
}


void dd_SelectNextHalfspace0(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*A natural way to choose the next hyperplane.  Simply the largest index*/
  long i;
  dd_boolean determined;

  i = cone->m;
  determined = dd_FALSE;
  do {
    if (set_member(i, excluded))
      i--;
    else
      determined = dd_TRUE;
  } while (!determined && i>=1);
  if (determined) 
    *hnext = i;
  else
    *hnext = 0;
}

void dd_SelectNextHalfspace1(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*Natural way to choose the next hyperplane.  Simply the least index*/
  long i;
  dd_boolean determined;

  i = 1;
  determined = dd_FALSE;
  do {
    if (set_member(i, excluded))
      i++;
    else
      determined = dd_TRUE;
  } while (!determined && i<=cone->m);
  if (determined) 
    *hnext = i;
  else 
    *hnext=0;
}

void dd_SelectNextHalfspace2(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*Choose the next hyperplane with maximum infeasibility*/
  long i, fea, inf, infmin, fi=0;   /*feasibility and infeasibility numbers*/

  infmin = cone->RayCount + 1;
  for (i = 1; i <= cone->m; i++) {
    if (!set_member(i, excluded)) {
      dd_FeasibilityIndices(&fea, &inf, i, cone);
      if (inf < infmin) {
	infmin = inf;
	fi = fea;
	*hnext = i;
      }
    }
  }
  if (dd_debug) {
    fprintf(stderr,"*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infmin, fi);
  }
}

void dd_SelectNextHalfspace3(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*Choose the next hyperplane with maximum infeasibility*/
  long i, fea, inf, infmax, fi=0;   /*feasibility and infeasibility numbers*/
  dd_boolean localdebug=dd_debug;

  infmax = -1;
  for (i = 1; i <= cone->m; i++) {
    if (!set_member(i, excluded)) {
      dd_FeasibilityIndices(&fea, &inf, i, cone);
      if (inf > infmax) {
	infmax = inf;
	fi = fea;
	*hnext = i;
      }
    }
  }
  if (localdebug) {
    fprintf(stderr,"*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infmax, fi);
  }
}

void dd_SelectNextHalfspace4(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*Choose the next hyperplane with the most unbalanced cut*/
  long i, fea, inf, max, tmax, fi=0, infi=0;
      /*feasibility and infeasibility numbers*/

  max = -1;
  for (i = 1; i <= cone->m; i++) {
    if (!set_member(i, excluded)) {
      dd_FeasibilityIndices(&fea, &inf, i, cone);
      if (fea <= inf)
        tmax = inf;
      else
        tmax = fea;
      if (tmax > max) {
        max = tmax;
        fi = fea;
        infi = inf;
        *hnext = i;
      }
    }
  }
  if (!dd_debug)
    return;
  if (max == fi) {
    fprintf(stderr,"*infeasible rays (min) =%5ld, #feas rays =%5ld\n", infi, fi);
  } else {
    fprintf(stderr,"*infeasible rays (max) =%5ld, #feas rays =%5ld\n", infi, fi);
  }
}

void dd_SelectNextHalfspace5(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*Choose the next hyperplane which is lexico-min*/
  long i, minindex;
  mytype *v1, *v2;

  minindex = 0;
  v1 = NULL;
  for (i = 1; i <= cone->m; i++) {
    if (!set_member(i, excluded)) {
	  v2 = cone->A[i - 1];
      if (minindex == 0) {
	    minindex = i;
	    v1=v2;
      } else if (dd_LexSmaller(v2,v1,cone->d)) {
        minindex = i;
	    v1=v2;
      }
    }
  }
  *hnext = minindex;
}


void dd_SelectNextHalfspace6(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hnext)
{
  /*Choose the next hyperplane which is lexico-max*/
  long i, maxindex;
  mytype *v1, *v2;

  maxindex = 0;
  v1 = NULL;
  for (i = 1; i <= cone->m; i++) {
    if (!set_member(i, excluded)) {
      v2= cone->A[i - 1];
      if (maxindex == 0) {
        maxindex = i;
        v1=v2;
      } else if (dd_LexLarger(v2, v1, cone->d)) {
        maxindex = i;
        v1=v2;
     }
    }
  }
  *hnext = maxindex;
}

void dd_UniqueRows(dd_rowindex OV, long p, long r, dd_Amatrix A, long dmax, dd_rowset preferred, long *uniqrows)
{
 /* Select a subset of rows of A (with range [p, q] up to dimension dmax) by removing duplicates.
    When a row subset preferred is nonempty, those row indices in the set have priority.  If
    two priority rows define the same row vector, one is chosen.
    For a selected unique row i, OV[i] returns a new position of the unique row i. 
    For other nonuniqu row i, OV[i] returns a negative of the original row j dominating i.
    Thus the original contents of OV[p..r] will be rewritten.  Other components remain the same.
    *uniqrows returns the number of unique rows.
*/
  long i,iuniq,j;
  mytype *x;
  dd_boolean localdebug=dd_FALSE;
  
  if (p<=0 || p > r) goto _L99;
  iuniq=p; j=1;  /* the first unique row candidate */
  x=A[p-1];
  OV[p]=j;  /* tentative row index of the candidate */
  for (i=p+1;i<=r; i++){
    if (!dd_LexEqual(x,A[i-1],dmax)) {
      /* a new row vector found. */
      iuniq=i;
      j=j+1;
      OV[i]=j;    /* Tentatively register the row i.  */
      x=A[i-1];
    } else {
      /* rows are equal */
      if (set_member(i, preferred) && !set_member(iuniq, preferred)){
        OV[iuniq]=-i;  /* the row iuniq is dominated by the row i */
        iuniq=i;  /* the row i is preferred.  Change the candidate. */
        OV[i]=j;  /* the row i is tentatively registered. */
        x=A[i-1];
      } else {
        OV[i]=-iuniq;  /* the row iuniq is dominated by the row i */
      }
    }
  }
  *uniqrows=j;
  if (localdebug){
    printf("The number of unique rows are %ld\n with order vector",*uniqrows);
    for (i=p;i<=r; i++){
      printf(" %ld:%ld ",i,OV[i]);
    }
    printf("\n");
  }
  _L99:;
}

long dd_Partition(dd_rowindex OV, long p, long r, dd_Amatrix A, long dmax)
{
  mytype *x;
  long i,j,ovi;
  
  x=A[OV[p]-1];

  i=p-1;
  j=r+1;
  while (dd_TRUE){
    do{
      j--;
    } while (dd_LexLarger(A[OV[j]-1],x,dmax));
    do{
      i++;
    } while (dd_LexSmaller(A[OV[i]-1],x,dmax));
    if (i<j){
      ovi=OV[i];
      OV[i]=OV[j];
      OV[j]=ovi;
    }
    else{
      return j;
    }
  }
}

void dd_QuickSort(dd_rowindex OV, long p, long r, dd_Amatrix A, long dmax)
{
  long q;
  
  if (p < r){
    q = dd_Partition(OV, p, r, A, dmax);
    dd_QuickSort(OV, p, q, A, dmax);
    dd_QuickSort(OV, q+1, r, A, dmax);
  }
}


#ifndef RAND_MAX 
#define RAND_MAX 32767 
#endif

void dd_RandomPermutation(dd_rowindex OV, long t, unsigned int seed)
{
  long k,j,ovj;
  double u,xk,r,rand_max=(double) RAND_MAX;
  dd_boolean localdebug=dd_FALSE;

  srand(seed);
  for (j=t; j>1 ; j--) {
    r=rand();
    u=r/rand_max;
    xk=(double)(j*u +1);
    k=(long)xk;
    if (localdebug) fprintf(stderr,"u=%g, k=%ld, r=%g, randmax= %g\n",u,k,r,rand_max);
    ovj=OV[j];
    OV[j]=OV[k];
    OV[k]=ovj;
    if (localdebug) fprintf(stderr,"row %ld is exchanged with %ld\n",j,k); 
  }
}

void dd_ComputeRowOrderVector(dd_ConePtr cone)
{
  long i,itemp;

  cone->OrderVector[0]=0;
  switch (cone->HalfspaceOrder){
  case dd_MaxIndex:
    for(i=1; i<=cone->m; i++) cone->OrderVector[i]=cone->m-i+1;
    break;

  case dd_MinIndex: 
    for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;    
    break;

  case dd_LexMin: case dd_MinCutoff: case dd_MixCutoff: case dd_MaxCutoff:
    for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;
    dd_RandomPermutation(cone->OrderVector, cone->m, cone->rseed);
    dd_QuickSort(cone->OrderVector, 1, cone->m, cone->A, cone->d);
    break;

  case dd_LexMax:
    for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;
    dd_RandomPermutation(cone->OrderVector, cone->m, cone->rseed);
    dd_QuickSort(cone->OrderVector, 1, cone->m, cone->A, cone->d);
    for(i=1; i<=cone->m/2;i++){   /* just reverse the order */
      itemp=cone->OrderVector[i];
      cone->OrderVector[i]=cone->OrderVector[cone->m-i+1];
      cone->OrderVector[cone->m-i+1]=itemp;
    }
    break;

  case dd_RandomRow:
    for(i=1; i<=cone->m; i++) cone->OrderVector[i]=i;
    dd_RandomPermutation(cone->OrderVector, cone->m, cone->rseed);
    break;

  }
}

void dd_UpdateRowOrderVector(dd_ConePtr cone, dd_rowset PriorityRows)
/* Update the RowOrder vector to shift selected rows
in highest order.
*/
{
  dd_rowrange i,j,k,j1=0,oj=0;
  long rr;
  dd_boolean found, localdebug=dd_FALSE;
  
  if (dd_debug) localdebug=dd_TRUE;
  found=dd_TRUE;
  rr=set_card(PriorityRows);
  if (localdebug) set_fwrite(stderr,PriorityRows);
  for (i=1; i<=rr; i++){
    found=dd_FALSE;
    for (j=i; j<=cone->m && !found; j++){
      oj=cone->OrderVector[j];
      if (set_member(oj, PriorityRows)){
        found=dd_TRUE;
        if (localdebug) fprintf(stderr,"%ldth in sorted list (row %ld) is in PriorityRows\n", j, oj);
        j1=j;
      }
    }
    if (found){
      if (j1>i) {
        /* shift everything lower: ov[i]->cone->ov[i+1]..ov[j1-1]->cone->ov[j1] */
        for (k=j1; k>=i; k--) cone->OrderVector[k]=cone->OrderVector[k-1];
        cone->OrderVector[i]=oj;
        if (localdebug){
          fprintf(stderr,"OrderVector updated to:\n");
          for (j = 1; j <= cone->m; j++) fprintf(stderr," %2ld", cone->OrderVector[j]);
          fprintf(stderr,"\n");
        }
      }
    } else {
      fprintf(stderr,"UpdateRowOrder: Error.\n");
      goto _L99;
    }
  }
_L99:;
}

void dd_SelectPreorderedNext(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hh)
{
  dd_rowrange i,k;
  
  *hh=0;
  for (i=1; i<=cone->m && *hh==0; i++){
    k=cone->OrderVector[i];
    if (!set_member(k, excluded)) *hh=k ;
  }
}

void dd_SelectNextHalfspace(dd_ConePtr cone, dd_rowset excluded, dd_rowrange *hh)
{
  if (cone->PreOrderedRun){
    if (dd_debug) {
      fprintf(stderr,"debug dd_SelectNextHalfspace: Use PreorderNext\n");
    }
    dd_SelectPreorderedNext(cone, excluded, hh);
  }
  else {
    if (dd_debug) {
      fprintf(stderr,"debug dd_SelectNextHalfspace: Use DynamicOrderedNext\n");
    }

    switch (cone->HalfspaceOrder) {

    case dd_MaxIndex:
      dd_SelectNextHalfspace0(cone, excluded, hh);
      break;

    case dd_MinIndex:
      dd_SelectNextHalfspace1(cone, excluded, hh);
      break;

    case dd_MinCutoff:
      dd_SelectNextHalfspace2(cone, excluded, hh);
      break;

    case dd_MaxCutoff:
      dd_SelectNextHalfspace3(cone, excluded, hh);
      break;

    case dd_MixCutoff:
      dd_SelectNextHalfspace4(cone, excluded, hh);
      break;

    default:
      dd_SelectNextHalfspace0(cone, excluded, hh);
      break;
    }
  }
}

dd_boolean dd_Nonnegative(mytype val)
{
/*  if (val>=-dd_zero) return dd_TRUE;  */
  if (dd_cmp(val,dd_minuszero)>=0) return dd_TRUE;
  else return dd_FALSE;
}

dd_boolean dd_Nonpositive(mytype val)
{
/*  if (val<=dd_zero) return dd_TRUE;  */
  if (dd_cmp(val,dd_zero)<=0) return dd_TRUE;
  else return dd_FALSE;
}

dd_boolean dd_Positive(mytype val)
{
  return !dd_Nonpositive(val);
}

dd_boolean dd_Negative(mytype val)
{
  return !dd_Nonnegative(val);
}

dd_boolean dd_EqualToZero(mytype val)
{
  return (dd_Nonnegative(val) && dd_Nonpositive(val));
}

dd_boolean dd_Nonzero(mytype val)
{
  return (dd_Positive(val) || dd_Negative(val));
}

dd_boolean dd_Equal(mytype val1,mytype val2)
{
  return (!dd_Larger(val1,val2) && !dd_Smaller(val1,val2));
}

dd_boolean dd_Larger(mytype val1,mytype val2)
{
  mytype temp;
  dd_boolean answer;

  dd_init(temp);
  dd_sub(temp,val1, val2);
  answer=dd_Positive(temp);
  dd_clear(temp);
  return answer;
}

dd_boolean dd_Smaller(mytype val1,mytype val2)
{
  return dd_Larger(val2,val1);
}

void dd_abs(mytype absval, mytype val)
{
  if (dd_Negative(val)) dd_neg(absval,val);
  else dd_set(absval,val); 
}

void dd_LinearComb(mytype lc, mytype v1, mytype c1, mytype v2, mytype c2)
/*  lc := v1 * c1 + v2 * c2   */
{
  mytype temp;

  dd_init(temp);
  dd_mul(lc,v1,c1);
  dd_mul(temp,v2,c2); 
  dd_add(lc,lc,temp);
  dd_clear(temp);
}

void dd_InnerProduct(mytype prod, dd_colrange d, dd_Arow v1, dd_Arow v2)
{
  mytype temp;
  dd_colrange j;
  dd_boolean localdebug=dd_FALSE;

  dd_init(temp);
  dd_set_si(prod, 0);
  for (j = 0; j < d; j++){
    dd_mul(temp,v1[j],v2[j]); 
    dd_add(prod,prod,temp);
  }
  if (localdebug) {
    fprintf(stderr,"InnerProduct:\n"); 
    dd_WriteArow(stderr, v1, d);
    dd_WriteArow(stderr, v2, d);
    fprintf(stderr,"prod ="); 
    dd_WriteNumber(stderr, prod);
    fprintf(stderr,"\n");
  }
  
  dd_clear(temp);
}

/* end of cddcore.c */