src/bdd/cudd/cuddExact.c

Go to the documentation of this file.

#include "util_hack.h"
#include "cuddInt.h"

/*---------------------------------------------------------------------------*/
/* Constant declarations                                                     */
/*---------------------------------------------------------------------------*/


/*---------------------------------------------------------------------------*/
/* Stucture declarations                                                     */
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/* Type declarations                                                         */
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/* Variable declarations                                                     */
/*---------------------------------------------------------------------------*/

#ifndef lint
static char rcsid[] DD_UNUSED = "$Id: cuddExact.c,v 1.1.1.1 2003/02/24 22:23:52 wjiang Exp $";
#endif

#ifdef DD_STATS
static int ddTotalShuffles;
#endif

/*---------------------------------------------------------------------------*/
/* Macro declarations                                                        */
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/* Static function prototypes                                                */
/*---------------------------------------------------------------------------*/

static int getMaxBinomial ARGS((int n));
static int gcd ARGS((int x, int y));
static DdHalfWord ** getMatrix ARGS((int rows, int cols));
static void freeMatrix ARGS((DdHalfWord **matrix));
static int getLevelKeys ARGS((DdManager *table, int l));
static int ddShuffle ARGS((DdManager *table, DdHalfWord *permutation, int lower, int upper));
static int ddSiftUp ARGS((DdManager *table, int x, int xLow));
static int updateUB ARGS((DdManager *table, int oldBound, DdHalfWord *bestOrder, int lower, int upper));
static int ddCountRoots ARGS((DdManager *table, int lower, int upper));
static void ddClearGlobal ARGS((DdManager *table, int lower, int maxlevel));
static int computeLB ARGS((DdManager *table, DdHalfWord *order, int roots, int cost, int lower, int upper, int level));
static int updateEntry ARGS((DdManager *table, DdHalfWord *order, int level, int cost, DdHalfWord **orders, int *costs, int subsets, char *mask, int lower, int upper));
static void pushDown ARGS((DdHalfWord *order, int j, int level));
static DdHalfWord * initSymmInfo ARGS((DdManager *table, int lower, int upper));
static int checkSymmInfo ARGS((DdManager *table, DdHalfWord *symmInfo, int index, int level));

/*---------------------------------------------------------------------------*/
/* Definition of exported functions                                          */
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/* Definition of internal functions                                          */
/*---------------------------------------------------------------------------*/


int
cuddExact(
  DdManager * table,
  int  lower,
  int  upper)
{
    int k, i, j;
    int maxBinomial, oldSubsets, newSubsets;
    int subsetCost;
    int size;                   /* number of variables to be reordered */
    int unused, nvars, level, result;
    int upperBound, lowerBound, cost;
    int roots;
    char *mask = NULL;
    DdHalfWord  *symmInfo = NULL;
    DdHalfWord **newOrder = NULL;
    DdHalfWord **oldOrder = NULL;
    int *newCost = NULL;
    int *oldCost = NULL;
    DdHalfWord **tmpOrder;
    int *tmpCost;
    DdHalfWord *bestOrder = NULL;
    DdHalfWord *order;
#ifdef DD_STATS
    int  ddTotalSubsets;
#endif

    /* Restrict the range to be reordered by excluding unused variables
    ** at the two ends. */
    while (table->subtables[lower].keys == 1 &&
           table->vars[table->invperm[lower]]->ref == 1 &&
           lower < upper)
        lower++;
    while (table->subtables[upper].keys == 1 &&
           table->vars[table->invperm[upper]]->ref == 1 &&
           lower < upper)
        upper--;
    if (lower == upper) return(1); /* trivial problem */

    /* Apply symmetric sifting to get a good upper bound and to extract
    ** symmetry information. */
    result = cuddSymmSiftingConv(table,lower,upper);
    if (result == 0) goto cuddExactOutOfMem;

#ifdef DD_STATS
    (void) fprintf(table->out,"\n");
    ddTotalShuffles = 0;
    ddTotalSubsets = 0;
#endif

    /* Initialization. */
    nvars = table->size;
    size = upper - lower + 1;
    /* Count unused variable among those to be reordered.  This is only
    ** used to compute maxBinomial. */
    unused = 0;
    for (i = lower + 1; i < upper; i++) {
        if (table->subtables[i].keys == 1 &&
            table->vars[table->invperm[i]]->ref == 1)
            unused++;
    }

    /* Find the maximum number of subsets we may have to store. */
    maxBinomial = getMaxBinomial(size - unused);
    if (maxBinomial == -1) goto cuddExactOutOfMem;

    newOrder = getMatrix(maxBinomial, size);
    if (newOrder == NULL) goto cuddExactOutOfMem;

    newCost = ALLOC(int, maxBinomial);
    if (newCost == NULL) goto cuddExactOutOfMem;

    oldOrder = getMatrix(maxBinomial, size);
    if (oldOrder == NULL) goto cuddExactOutOfMem;

    oldCost = ALLOC(int, maxBinomial);
    if (oldCost == NULL) goto cuddExactOutOfMem;

    bestOrder = ALLOC(DdHalfWord, size);
    if (bestOrder == NULL) goto cuddExactOutOfMem;

    mask = ALLOC(char, nvars);
    if (mask == NULL) goto cuddExactOutOfMem;

    symmInfo = initSymmInfo(table, lower, upper);
    if (symmInfo == NULL) goto cuddExactOutOfMem;

    roots = ddCountRoots(table, lower, upper);

    /* Initialize the old order matrix for the empty subset and the best
    ** order to the current order. The cost for the empty subset includes
    ** the cost of the levels between upper and the constants. These levels
    ** are not going to change. Hence, we count them only once.
    */
    oldSubsets = 1;
    for (i = 0; i < size; i++) {
        oldOrder[0][i] = bestOrder[i] = (DdHalfWord) table->invperm[i+lower];
    }
    subsetCost = table->constants.keys;
    for (i = upper + 1; i < nvars; i++)
        subsetCost += getLevelKeys(table,i);
    oldCost[0] = subsetCost;
    /* The upper bound is initialized to the current size of the BDDs. */
    upperBound = table->keys - table->isolated;

    /* Now consider subsets of increasing size. */
    for (k = 1; k <= size; k++) {
#if DD_STATS
        (void) fprintf(table->out,"Processing subsets of size %d\n", k);
        fflush(table->out);
#endif
        newSubsets = 0;
        level = size - k;               /* offset of first bottom variable */

        for (i = 0; i < oldSubsets; i++) { /* for each subset of size k-1 */
            order = oldOrder[i];
            cost = oldCost[i];
            lowerBound = computeLB(table, order, roots, cost, lower, upper,
                                   level);
            if (lowerBound >= upperBound)
                continue;
            /* Impose new order. */
            result = ddShuffle(table, order, lower, upper);
            if (result == 0) goto cuddExactOutOfMem;
            upperBound = updateUB(table,upperBound,bestOrder,lower,upper);
            /* For each top bottom variable. */
            for (j = level; j >= 0; j--) {
                /* Skip unused variables. */
                if (table->subtables[j+lower-1].keys == 1 &&
                    table->vars[table->invperm[j+lower-1]]->ref == 1) continue;
                /* Find cost under this order. */
                subsetCost = cost + getLevelKeys(table, lower + level);
                newSubsets = updateEntry(table, order, level, subsetCost,
                                         newOrder, newCost, newSubsets, mask,
                                         lower, upper);
                if (j == 0)
                    break;
                if (checkSymmInfo(table, symmInfo, order[j-1], level) == 0)
                    continue;
                pushDown(order,j-1,level);
                /* Impose new order. */
                result = ddShuffle(table, order, lower, upper);
                if (result == 0) goto cuddExactOutOfMem;
                upperBound = updateUB(table,upperBound,bestOrder,lower,upper);
            } /* for each bottom variable */
        } /* for each subset of size k */

        /* New orders become old orders in preparation for next iteration. */
        tmpOrder = oldOrder; tmpCost = oldCost;
        oldOrder = newOrder; oldCost = newCost;
        newOrder = tmpOrder; newCost = tmpCost;
#ifdef DD_STATS
        ddTotalSubsets += newSubsets;
#endif
        oldSubsets = newSubsets;
    }
    result = ddShuffle(table, bestOrder, lower, upper);
    if (result == 0) goto cuddExactOutOfMem;
#ifdef DD_STATS
#ifdef DD_VERBOSE
    (void) fprintf(table->out,"\n");
#endif
    (void) fprintf(table->out,"#:S_EXACT   %8d: total subsets\n",
                   ddTotalSubsets);
    (void) fprintf(table->out,"#:H_EXACT   %8d: total shuffles",
                   ddTotalShuffles);
#endif

    freeMatrix(newOrder);
    freeMatrix(oldOrder);
    FREE(bestOrder);
    FREE(oldCost);
    FREE(newCost);
    FREE(symmInfo);
    FREE(mask);
    return(1);

cuddExactOutOfMem:

    if (newOrder != NULL) freeMatrix(newOrder);
    if (oldOrder != NULL) freeMatrix(oldOrder);
    if (bestOrder != NULL) FREE(bestOrder);
    if (oldCost != NULL) FREE(oldCost);
    if (newCost != NULL) FREE(newCost);
    if (symmInfo != NULL) FREE(symmInfo);
    if (mask != NULL) FREE(mask);
    table->errorCode = CUDD_MEMORY_OUT;
    return(0);

} /* end of cuddExact */


static int
getMaxBinomial(
  int  n)
{
    int *numerator;
    int i, j, k, y, g, result;

    k = (n & ~1) >> 1;

    numerator = ALLOC(int,k);
    if (numerator == NULL) return(-1);
    
    for (i = 0; i < k; i++)
        numerator[i] = n - i;
    
    for (i = k; i > 1; i--) {
        y = i;
        for (j = 0; j < k; j++) {
            if (numerator[j] == 1) continue;
            g = gcd(numerator[j], y);
            if (g != 1) {
                numerator[j] /= g;
                if (y == g) break;
                y /= g;
            }
        }
    }

    result = 1;
    for (i = 0; i < k; i++)
        result *= numerator[i];

    FREE(numerator);
    return(result);

}  /* end of getMaxBinomial */


static int
gcd(
  int  x,
  int  y)
{
    int a;
    int b;
    int lsbMask;

    /* GCD(n,0) = n. */
    if (x == 0) return(y);
    if (y == 0) return(x);

    a = x; b = y; lsbMask = 1;
    
    /* Here both a and b are != 0. The iteration maintains this invariant.
    ** Hence, we only need to check for when they become equal.
    */
    while (a != b) {
        if (a & lsbMask) {
            if (b & lsbMask) {  /* both odd */
                if (a < b) {
                    b = (b - a) >> 1;
                } else {
                    a = (a - b) >> 1;
                }
            } else {            /* a odd, b even */
                b >>= 1;
            }
        } else {
            if (b & lsbMask) {  /* a even, b odd */
                a >>= 1;
            } else {            /* both even */
                lsbMask <<= 1;
            }
        }
    }

    return(a);

} /* end of gcd */


static DdHalfWord **
getMatrix(
  int  rows /* number of rows */,
  int  cols /* number of columns */)
{
    DdHalfWord **matrix;
    int i;

    if (cols*rows == 0) return(NULL);
    matrix = ALLOC(DdHalfWord *, rows);
    if (matrix == NULL) return(NULL);
    matrix[0] = ALLOC(DdHalfWord, cols*rows);
    if (matrix[0] == NULL) return(NULL);
    for (i = 1; i < rows; i++) {
        matrix[i] = matrix[i-1] + cols;
    }
    return(matrix);

} /* end of getMatrix */


static void
freeMatrix(
  DdHalfWord ** matrix)
{
    FREE(matrix[0]);
    FREE(matrix);
    return;

} /* end of freeMatrix */


static int
getLevelKeys(
  DdManager * table,
  int  l)
{
    int isolated;
    int x;        /* x is an index */

    x = table->invperm[l];
    isolated = table->vars[x]->ref == 1;

    return(table->subtables[l].keys - isolated);

} /* end of getLevelKeys */


static int
ddShuffle(
  DdManager * table,
  DdHalfWord * permutation,
  int  lower,
  int  upper)
{
    DdHalfWord  index;
    int         level;
    int         position;
    int         numvars;
    int         result;
#ifdef DD_STATS
    long        localTime;
    int         initialSize;
#ifdef DD_VERBOSE
    int         finalSize;
#endif
    int         previousSize;
#endif

#ifdef DD_STATS
    localTime = util_cpu_time();
    initialSize = table->keys - table->isolated;
#endif

    numvars = table->size;

#if 0
    (void) fprintf(table->out,"%d:", ddTotalShuffles);
    for (level = 0; level < numvars; level++) {
        (void) fprintf(table->out," %d", table->invperm[level]);
    }
    (void) fprintf(table->out,"\n");
#endif

    for (level = 0; level <= upper - lower; level++) {
        index = permutation[level];
        position = table->perm[index];
#ifdef DD_STATS
        previousSize = table->keys - table->isolated;
#endif
        result = ddSiftUp(table,position,level+lower);
        if (!result) return(0);
    }

#ifdef DD_STATS
    ddTotalShuffles++;
#ifdef DD_VERBOSE
    finalSize = table->keys - table->isolated;
    if (finalSize < initialSize) {
        (void) fprintf(table->out,"-");
    } else if (finalSize > initialSize) {
        (void) fprintf(table->out,"+");
    } else {
        (void) fprintf(table->out,"=");
    }
    if ((ddTotalShuffles & 63) == 0) (void) fprintf(table->out,"\n");
    fflush(table->out);
#endif
#endif

    return(1);

} /* end of ddShuffle */


static int
ddSiftUp(
  DdManager * table,
  int  x,
  int  xLow)
{
    int        y;
    int        size;

    y = cuddNextLow(table,x);
    while (y >= xLow) {
        size = cuddSwapInPlace(table,y,x);
        if (size == 0) {
            return(0);
        }
        x = y;
        y = cuddNextLow(table,x);
    }
    return(1);

} /* end of ddSiftUp */


static int
updateUB(
  DdManager * table,
  int  oldBound,
  DdHalfWord * bestOrder,
  int  lower,
  int  upper)
{
    int i;
    int newBound = table->keys - table->isolated;

    if (newBound < oldBound) {
#ifdef DD_STATS
        (void) fprintf(table->out,"New upper bound = %d\n", newBound);
        fflush(table->out);
#endif
        for (i = lower; i <= upper; i++)
            bestOrder[i-lower] = (DdHalfWord) table->invperm[i];
        return(newBound);
    } else {
        return(oldBound);
    }

} /* end of updateUB */


static int
ddCountRoots(
  DdManager * table,
  int  lower,
  int  upper)
{
    int i,j;
    DdNode *f;
    DdNodePtr *nodelist;
    DdNode *sentinel = &(table->sentinel);
    int slots;
    int roots = 0;
    int maxlevel = lower;

    for (i = lower; i <= upper; i++) {
        nodelist = table->subtables[i].nodelist;
        slots = table->subtables[i].slots;
        for (j = 0; j < slots; j++) {
            f = nodelist[j];
            while (f != sentinel) {
                /* A node is a root of the DAG if it cannot be
                ** reached by nodes above it. If a node was never
                ** reached during the previous depth-first searches,
                ** then it is a root, and we start a new depth-first
                ** search from it.
                */
                if (!Cudd_IsComplement(f->next)) {
                    if (f != table->vars[f->index]) {
                        roots++;
                    }
                }
                if (!Cudd_IsConstant(cuddT(f))) {
                    cuddT(f)->next = Cudd_Complement(cuddT(f)->next);
                    if (table->perm[cuddT(f)->index] > maxlevel)
                        maxlevel = table->perm[cuddT(f)->index];
                }
                if (!Cudd_IsConstant(cuddE(f))) {
                    Cudd_Regular(cuddE(f))->next =
                        Cudd_Complement(Cudd_Regular(cuddE(f))->next);
                    if (table->perm[Cudd_Regular(cuddE(f))->index] > maxlevel)
                        maxlevel = table->perm[Cudd_Regular(cuddE(f))->index];
                }
                f = Cudd_Regular(f->next);
            }
        }
    }
    ddClearGlobal(table, lower, maxlevel);

    return(roots);

} /* end of ddCountRoots */


static void
ddClearGlobal(
  DdManager * table,
  int  lower,
  int  maxlevel)
{
    int i,j;
    DdNode *f;
    DdNodePtr *nodelist;
    DdNode *sentinel = &(table->sentinel);
    int slots;

    for (i = lower; i <= maxlevel; i++) {
        nodelist = table->subtables[i].nodelist;
        slots = table->subtables[i].slots;
        for (j = 0; j < slots; j++) {
            f = nodelist[j];
            while (f != sentinel) {
                f->next = Cudd_Regular(f->next);
                f = f->next;
            }
        }
    }

} /* end of ddClearGlobal */


static int
computeLB(
  DdManager * table             /* manager */,
  DdHalfWord * order            /* optimal order for the subset */,
  int  roots                    /* roots between lower and upper */,
  int  cost                     /* minimum cost for the subset */,
  int  lower                    /* lower level to be reordered */,
  int  upper                    /* upper level to be reordered */,
  int  level                    /* offset for the current top bottom var */
  )
{
    int i;
    int lb = cost;
    int lb1 = 0;
    int lb2;
    int support;
    DdHalfWord ref;

    /* The levels not involved in reordering are not going to change.
    ** Add their sizes to the lower bound.
    */
    for (i = 0; i < lower; i++) {
        lb += getLevelKeys(table,i);
    }
    /* If a variable is in the support, then there is going
    ** to be at least one node labeled by that variable.
    */
    for (i = lower; i <= lower+level; i++) {
        support = table->subtables[i].keys > 1 ||
            table->vars[order[i-lower]]->ref > 1;
        lb1 += support;
    }

    /* Estimate the number of nodes required to connect the roots to
    ** the nodes in the bottom part. */
    if (lower+level+1 < table->size) {
        if (lower+level < upper)
            ref = table->vars[order[level+1]]->ref;
        else
            ref = table->vars[table->invperm[upper+1]]->ref;
        lb2 = table->subtables[lower+level+1].keys -
            (ref > (DdHalfWord) 1) - roots;
    } else {
        lb2 = 0;
    }

    lb += lb1 > lb2 ? lb1 : lb2;

    return(lb);

} /* end of computeLB */


static int
updateEntry(
  DdManager * table,
  DdHalfWord * order,
  int  level,
  int  cost,
  DdHalfWord ** orders,
  int * costs,
  int  subsets,
  char * mask,
  int  lower,
  int  upper)
{
    int i, j;
    int size = upper - lower + 1;

    /* Build a mask that says what variables are in this subset. */
    for (i = lower; i <= upper; i++)
        mask[table->invperm[i]] = 0;
    for (i = level; i < size; i++)
        mask[order[i]] = 1;

    /* Check each subset until a match is found or all subsets are examined. */
    for (i = 0; i < subsets; i++) {
        DdHalfWord *subset = orders[i];
        for (j = level; j < size; j++) {
            if (mask[subset[j]] == 0)
                break;
        }
        if (j == size)          /* no mismatches: success */
            break;
    }
    if (i == subsets || cost < costs[i]) {              /* add or replace */
        for (j = 0; j < size; j++)
            orders[i][j] = order[j];
        costs[i] = cost;
        subsets += (i == subsets);
    }
    return(subsets);

} /* end of updateEntry */


static void
pushDown(
  DdHalfWord * order,
  int  j,
  int  level)
{
    int i;
    DdHalfWord tmp;

    tmp = order[j];
    for (i = j; i < level; i++) {
        order[i] = order[i+1];
    }
    order[level] = tmp;
    return;

} /* end of pushDown */


static DdHalfWord *
initSymmInfo(
  DdManager * table,
  int  lower,
  int  upper)
{
    int level, index, next, nextindex;
    DdHalfWord *symmInfo;

    symmInfo =  ALLOC(DdHalfWord, table->size);
    if (symmInfo == NULL) return(NULL);

    for (level = lower; level <= upper; level++) {
        index = table->invperm[level];
        next =  table->subtables[level].next;
        nextindex = table->invperm[next];
        symmInfo[index] = nextindex;
    }
    return(symmInfo);

} /* end of initSymmInfo */


static int
checkSymmInfo(
  DdManager * table,
  DdHalfWord * symmInfo,
  int  index,
  int  level)
{
    int i;

    i = symmInfo[index];
    while (i != index) {
        if (index < i && table->perm[i] <= level)
            return(0);
        i = symmInfo[i];
    }
    return(1);

} /* end of checkSymmInfo */

---