cddcore.c

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

  }
/* must add (by Sato) */
  free(cone->Edges);
  
  set_free(cone->GroundSet); 
  set_free(cone->EqualitySet); 
  set_free(cone->NonequalitySet); 
  set_free(cone->AddedHalfspaces); 
  set_free(cone->WeaklyAddedHalfspaces); 
  set_free(cone->InitialHalfspaces);
  free(cone->InitialRayIndex);
  free(cone->OrderVector);
  free(cone->newcol);

/* Fixed by Shawn Rusaw.  Originally it was cone->d instead of cone->d_alloc */
  dd_FreeBmatrix(cone->d_alloc,cone->B);
  dd_FreeBmatrix(cone->d_alloc,cone->Bsave);

/* Fixed by Marc Pfetsch 010219*/
  dd_FreeAmatrix(cone->m_alloc,cone->d_alloc,cone->A);
  cone->A = NULL;

  free(cone);
}

void dd_FreeDDMemory(dd_PolyhedraPtr poly)
{
  dd_FreeDDMemory0(poly->child);
  poly->child=NULL;
}

void dd_FreePolyhedra(dd_PolyhedraPtr poly)
{
  dd_bigrange i;

  if ((poly)->child != NULL) dd_FreeDDMemory(poly);
  dd_FreeAmatrix((poly)->m_alloc,poly->d_alloc, poly->A);
  dd_FreeArow((poly)->d_alloc,(poly)->c);
  free((poly)->EqualityIndex);
  if (poly->AincGenerated){
    for (i=1; i<=poly->m1; i++){
      set_free(poly->Ainc[i-1]);
    }
    free(poly->Ainc);
    set_free(poly->Ared);
    set_free(poly->Adom);
    poly->Ainc=NULL;
  }

  free(poly);
}

void dd_Normalize(dd_colrange d_size, mytype *V)
{
  long j,jmin=0;
  mytype temp,min;
  dd_boolean nonzerofound=dd_FALSE;

  if (d_size>0){
    dd_init(min);  dd_init(temp);
    dd_abs(min,V[0]);  jmin=0; /* set the minmizer to 0 */
    if (dd_Positive(min)) nonzerofound=dd_TRUE;
    for (j = 1; j < d_size; j++) {
      dd_abs(temp,V[j]);
      if (dd_Positive(temp)){
        if (!nonzerofound || dd_Smaller(temp,min)){
          nonzerofound=dd_TRUE;
          dd_set(min, temp);  jmin=j;
        }
      }
    }
    if (dd_Positive(min)){
      for (j = 0; j < d_size; j++) dd_div(V[j], V[j], min);
    }
    dd_clear(min); dd_clear(temp);
  }
}


void dd_ZeroIndexSet(dd_rowrange m_size, dd_colrange d_size, dd_Amatrix A, mytype *x, dd_rowset ZS)
{
  dd_rowrange i;
  mytype temp;

  /* Changed by Marc Pfetsch 010219 */
  dd_init(temp);
  set_emptyset(ZS);
  for (i = 1; i <= m_size; i++) {
    dd_AValue(&temp, d_size, A, x, i);
    if (dd_EqualToZero(temp)) set_addelem(ZS, i);
  }

  /* Changed by Marc Pfetsch 010219 */
  dd_clear(temp);
}

void dd_CopyBmatrix(dd_colrange d_size, dd_Bmatrix T, dd_Bmatrix TCOPY)
{
  dd_rowrange i;
  dd_colrange j;

  for (i=0; i < d_size; i++) {
    for (j=0; j < d_size; j++) {
      dd_set(TCOPY[i][j],T[i][j]);
    }
  }
}


void dd_CopyArow(mytype *acopy, mytype *a, dd_colrange d)
{
  dd_colrange j;

  for (j = 0; j < d; j++) {
    dd_set(acopy[j],a[j]);
  }
}

void dd_CopyNormalizedArow(mytype *acopy, mytype *a, dd_colrange d)
{
  dd_CopyArow(acopy, a, d);
  dd_Normalize(d,acopy);
}

void dd_CopyAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d)
{
  dd_rowrange i;

  for (i = 0; i< m; i++) {
    dd_CopyArow(Acopy[i],A[i],d);
  }
}

void dd_CopyNormalizedAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d)
{
  dd_rowrange i;

  for (i = 0; i< m; i++) {
    dd_CopyNormalizedArow(Acopy[i],A[i],d);
  }
}

void dd_PermuteCopyAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d, dd_rowindex roworder)
{
  dd_rowrange i;

  for (i = 1; i<= m; i++) {
    dd_CopyArow(Acopy[i-1],A[roworder[i]-1],d);
  }
}

void dd_PermutePartialCopyAmatrix(mytype **Acopy, mytype **A, dd_rowrange m, dd_colrange d, dd_rowindex roworder,dd_rowrange p, dd_rowrange q)
{
 /* copy the rows of A whose roworder is positive.  roworder[i] is the row index of the copied row. */
  dd_rowrange i,k;

  k=0;
  for (i = 1; i<= m; i++) {
    if (roworder[i]>0) dd_CopyArow(Acopy[roworder[i]-1],A[i-1],d);
  }
}

void dd_InitializeArow(dd_colrange d,dd_Arow *a)
{
  dd_colrange j;

  if (d>0) *a=(mytype*) calloc(d,sizeof(mytype));
  for (j = 0; j < d; j++) {
      dd_init((*a)[j]);
  }
}

void dd_InitializeAmatrix(dd_rowrange m,dd_colrange d,dd_Amatrix *A)
{
  dd_rowrange i;

  if (m>0) (*A)=(mytype**) calloc(m,sizeof(mytype*));
  for (i = 0; i < m; i++) {
    dd_InitializeArow(d,&((*A)[i]));
  }
}

void dd_FreeAmatrix(dd_rowrange m,dd_colrange d,dd_Amatrix A)
{
  dd_rowrange i;
  dd_colrange j;

  for (i = 0; i < m; i++) {
    for (j = 0; j < d; j++) {
      dd_clear(A[i][j]);
    }
  }
  if (A!=NULL) {
    for (i = 0; i < m; i++) {
      free(A[i]);
    }
    free(A);
  }
}

void dd_FreeArow(dd_colrange d, dd_Arow a)
{
  dd_colrange j;

  for (j = 0; j < d; j++) {
    dd_clear(a[j]);
  }
  free(a);
}


void dd_InitializeBmatrix(dd_colrange d,dd_Bmatrix *B)
{
  dd_colrange i,j;

  (*B)=(mytype**) calloc(d,sizeof(mytype*));
  for (j = 0; j < d; j++) {
    (*B)[j]=(mytype*) calloc(d,sizeof(mytype));
  }
  for (i = 0; i < d; i++) {
    for (j = 0; j < d; j++) {
      dd_init((*B)[i][j]);
    }
  }
}

void dd_FreeBmatrix(dd_colrange d,dd_Bmatrix B)
{
  dd_colrange i,j;

  for (i = 0; i < d; i++) {
    for (j = 0; j < d; j++) {
      dd_clear(B[i][j]);
    }
  }
  if (B!=NULL) {
    for (j = 0; j < d; j++) {
      free(B[j]);
    }
    free(B);
  }
}

dd_SetFamilyPtr dd_CreateSetFamily(dd_bigrange fsize, dd_bigrange ssize)
{
  dd_SetFamilyPtr F;
  dd_bigrange i,f0,f1,s0,s1;

  if (fsize<=0) {
    f0=0; f1=1;  
    /* if fsize<=0, the fsize is set to zero and the created size is one */
  } else {
    f0=fsize; f1=fsize;
  }
  if (ssize<=0) {
    s0=0; s1=1;  
    /* if ssize<=0, the ssize is set to zero and the created size is one */
  } else {
    s0=ssize; s1=ssize;
  }

  F=(dd_SetFamilyPtr) malloc (sizeof(dd_SetFamilyType));
  F->set=(set_type*) calloc(f1,sizeof(set_type));
  for (i=0; i<f1; i++) {
    set_initialize(&(F->set[i]), s1);
  }
  F->famsize=f0;
  F->setsize=s0;
  return F;
}


void dd_FreeSetFamily(dd_SetFamilyPtr F)
{
  dd_bigrange i,f1;

  if (F!=NULL){
    if (F->famsize<=0) f1=1; else f1=F->famsize; 
      /* the smallest created size is one */
    for (i=0; i<f1; i++) {
      set_free(F->set[i]);
    }
    free(F->set);
    free(F);
  }
}

dd_MatrixPtr dd_CreateMatrix(dd_rowrange m_size,dd_colrange d_size)
{
  dd_MatrixPtr M;
  dd_rowrange m0,m1;
  dd_colrange d0,d1;

  if (m_size<=0){ 
    m0=0; m1=1;  
    /* if m_size <=0, the number of rows is set to zero, the actual size is 1 */
  } else {
    m0=m_size; m1=m_size;
  }
  if (d_size<=0){ 
    d0=0; d1=1;  
    /* if d_size <=0, the number of cols is set to zero, the actual size is 1 */
  } else {
    d0=d_size; d1=d_size;
  }
  M=(dd_MatrixPtr) malloc (sizeof(dd_MatrixType));
  dd_InitializeAmatrix(m1,d1,&(M->matrix));
  dd_InitializeArow(d1,&(M->rowvec));
  M->rowsize=m0;
  set_initialize(&(M->linset), m1);
  M->colsize=d0;
  M->objective=dd_LPnone;
  M->numbtype=dd_Unknown;
  M->representation=dd_Unspecified;
  return M;
}

void dd_FreeMatrix(dd_MatrixPtr M)
{
  dd_rowrange m1;
  dd_colrange d1;

  if (M!=NULL) {
    if (M->rowsize<=0) m1=1; else m1=M->rowsize;
    if (M->colsize<=0) d1=1; else d1=M->colsize;
    if (M!=NULL) {
      dd_FreeAmatrix(m1,d1,M->matrix);
      dd_FreeArow(d1,M->rowvec);
      set_free(M->linset);
      free(M);
    }
  }
}

void dd_SetToIdentity(dd_colrange d_size, dd_Bmatrix T)
{
  dd_colrange j1, j2;

  for (j1 = 1; j1 <= d_size; j1++) {
    for (j2 = 1; j2 <= d_size; j2++) {
      if (j1 == j2)
        dd_set(T[j1 - 1][j2 - 1],dd_one);
      else
        dd_set(T[j1 - 1][j2 - 1],dd_purezero);
    }
  }
}

void dd_ColumnReduce(dd_ConePtr cone)
{
  dd_colrange j,j1=0;
  dd_rowrange i;
  dd_boolean localdebug=dd_FALSE;

  for (j=1;j<=cone->d;j++) {
    if (cone->InitialRayIndex[j]>0){
      j1=j1+1;
      if (j1<j) {
        for (i=1; i<=cone->m; i++) dd_set(cone->A[i-1][j1-1],cone->A[i-1][j-1]);
        cone->newcol[j]=j1;
        if (localdebug){
          fprintf(stderr,"shifting the column %ld to column %ld\n", j, j1);
        }
          /* shifting forward */
      }
    } else {
      cone->newcol[j]=0;
      if (localdebug) {
        fprintf(stderr,"a generator (or an equation) of the linearity space: ");
        for (i=1; i<=cone->d; i++) dd_WriteNumber(stderr, cone->B[i-1][j-1]);
        fprintf(stderr,"\n");
      }
    }
  }
  cone->d=j1;  /* update the dimension. cone->d_orig remembers the old. */
  dd_CopyBmatrix(cone->d_orig, cone->B, cone->Bsave);  
    /* save the dual basis inverse as Bsave.  This matrix contains the linearity space generators. */
  cone->ColReduced=dd_TRUE;
}

long dd_MatrixRank(dd_MatrixPtr M, dd_rowset ignoredrows, dd_colset ignoredcols, dd_rowset *rowbasis, dd_colset *colbasis)
{
  dd_boolean stop, chosen, localdebug=dd_debug;
  dd_rowset NopivotRow,PriorityRow;
  dd_colset ColSelected;
  dd_Bmatrix B;
  dd_rowindex roworder;
  dd_rowrange r;
  dd_colrange s;
  long rank;

  rank = 0;
  stop = dd_FALSE;
  set_initialize(&ColSelected, M->colsize);
  set_initialize(&NopivotRow, M->rowsize);
  set_initialize(rowbasis, M->rowsize);
  set_initialize(colbasis, M->colsize);
  set_initialize(&PriorityRow, M->rowsize);
  set_copy(NopivotRow,ignoredrows);
  set_copy(ColSelected,ignoredcols);
  dd_InitializeBmatrix(M->colsize, &B);
  dd_SetToIdentity(M->colsize, B);
  roworder=(long *)calloc(M->rowsize+1,sizeof(long*));
  for (r=0; r<M->rowsize; r++){roworder[r+1]=r+1;
  }

  do {   /* Find a set of rows for a basis */
      dd_SelectPivot2(M->rowsize, M->colsize,M->matrix,B,dd_MinIndex,roworder,
       PriorityRow,M->rowsize, NopivotRow, ColSelected, &r, &s, &chosen);
      if (dd_debug && chosen) 
        fprintf(stderr,"Procedure dd_MatrixRank: pivot on (r,s) =(%ld, %ld).\n", r, s);
      if (chosen) {
        set_addelem(NopivotRow, r);
        set_addelem(*rowbasis, r);
        set_addelem(ColSelected, s);
        set_addelem(*colbasis, s);
        rank++;
        dd_GaussianColumnPivot(M->rowsize, M->colsize, M->matrix, B, r, s);
        if (localdebug) dd_WriteBmatrix(stderr,M->colsize,B);
      } else {
        stop=dd_TRUE;
      }
      if (rank==M->colsize) stop = dd_TRUE;
  } while (!stop);
  dd_FreeBmatrix(M->colsize,B);
  free(roworder);
  set_free(ColSelected);
  set_free(NopivotRow);
  set_free(PriorityRow);
  return rank;
}


void dd_FindBasis(dd_ConePtr cone, long *rank)
{
  dd_boolean stop, chosen, localdebug=dd_debug;
  dd_rowset NopivotRow;
  dd_colset ColSelected;
  dd_rowrange r;
  dd_colrange j,s;

  *rank = 0;
  stop = dd_FALSE;
  for (j=0;j<=cone->d;j++) cone->InitialRayIndex[j]=0;
  set_emptyset(cone->InitialHalfspaces);
  set_initialize(&ColSelected, cone->d);
  set_initialize(&NopivotRow, cone->m);
  set_copy(NopivotRow,cone->NonequalitySet);
  dd_SetToIdentity(cone->d, cone->B);
  do {   /* Find a set of rows for a basis */
      dd_SelectPivot2(cone->m, cone->d,cone->A,cone->B,cone->HalfspaceOrder,cone->OrderVector,
       cone->EqualitySet,cone->m, NopivotRow, ColSelected, &r, &s, &chosen);
      if (dd_debug && chosen) 
        fprintf(stderr,"Procedure dd_FindBasis: pivot on (r,s) =(%ld, %ld).\n", r, s);
      if (chosen) {
        set_addelem(cone->InitialHalfspaces, r);
        set_addelem(NopivotRow, r);
        set_addelem(ColSelected, s);
        cone->InitialRayIndex[s]=r;    /* cone->InitialRayIndex[s] stores the corr. row index */
        (*rank)++;
        dd_GaussianColumnPivot(cone->m, cone->d, cone->A, cone->B, r, s);
        if (localdebug) dd_WriteBmatrix(stderr,cone->d,cone->B);
      } else {
        stop=dd_TRUE;
      }
      if (*rank==cone->d) stop = dd_TRUE;
  } while (!stop);
  set_free(ColSelected);
  set_free(NopivotRow);
}


void dd_FindInitialRays(dd_ConePtr cone, dd_boolean *found)
{
  dd_rowset CandidateRows;
  dd_rowrange i;
  long rank;
  dd_RowOrderType roworder_save=dd_LexMin;

  *found = dd_FALSE;
  set_initialize(&CandidateRows, cone->m);
  if (cone->parent->InitBasisAtBottom==dd_TRUE) {
    roworder_save=cone->HalfspaceOrder;
    cone->HalfspaceOrder=dd_MaxIndex;
    cone->PreOrderedRun=dd_FALSE;
  }
  else cone->PreOrderedRun=dd_TRUE;
  if (dd_debug) dd_WriteBmatrix(stderr, cone->d, cone->B);
  for (i = 1; i <= cone->m; i++)
    if (!set_member(i,cone->NonequalitySet)) set_addelem(CandidateRows, i);
    /*all rows not in NonequalitySet are candidates for initial cone*/
  dd_FindBasis(cone, &rank);
  if (dd_debug) dd_WriteBmatrix(stderr, cone->d, cone->B);
  if (dd_debug) fprintf(stderr,"dd_FindInitialRays: rank of Amatrix = %ld\n", rank);
  cone->LinearityDim=cone->d - rank;
  if (dd_debug) fprintf(stderr,"Linearity Dimension = %ld\n", cone->LinearityDim);
  if (cone->LinearityDim > 0) {
     dd_ColumnReduce(cone);
     dd_FindBasis(cone, &rank);
  }
  if (!set_subset(cone->EqualitySet,cone->InitialHalfspaces)) {
    if (dd_debug) {
      fprintf(stderr,"Equality set is dependent. Equality Set and an initial basis:\n");
      set_fwrite(stderr,cone->EqualitySet);
      set_fwrite(stderr,cone->InitialHalfspaces);
    };
  }
  *found = dd_TRUE;
  set_free(CandidateRows);
  if (cone->parent->InitBasisAtBottom==dd_TRUE) {
    cone->HalfspaceOrder=roworder_save;
  }
  if (cone->HalfspaceOrder==dd_MaxCutoff||
      cone->HalfspaceOrder==dd_MinCutoff||
      cone->HalfspaceOrder==dd_MixCutoff){
    cone->PreOrderedRun=dd_FALSE;
  } else cone->PreOrderedRun=dd_TRUE;
}

void dd_CheckEquality(dd_colrange d_size, dd_RayPtr*RP1, dd_RayPtr*RP2, dd_boolean *equal)
{
  long j;

  if (dd_debug)
    fprintf(stderr,"Check equality of two rays\n");
  *equal = dd_TRUE;
  j = 1;
  while (j <= d_size && *equal) {