src/fsm/fsmHD.c

Go to the documentation of this file.

#include "fsmInt.h"
#include "bdd.h"
#include "ntm.h"
#include "img.h"

static char rcsid[] UNUSED = "$Id: fsmHD.c,v 1.32 2005/05/16 06:23:29 fabio Exp $";

/*---------------------------------------------------------------------------*/
/* Constant declarations                                                     */
/*---------------------------------------------------------------------------*/
/* thresholds for computation */
#define FSM_HD_DEADEND_RESIDUE_LIMIT 5000/*size beyond which residue is freed*/
#define FSM_HD_FRONTIER_APPROX_THRESHOLD 3500 /* default frontier approximation
                                                 threshold */
#define FSM_HD_REACHED_THRESHOLD  2000   /* Threshold below which reached is
                                         not decomposed */ 
#define FSM_HD_DISJ_SIZE 5000            /* default Size of disjunct whose image
                            may be computed in dead end computation. */
#define FSM_HD_DEADEND_MAX_SIZE_FACTOR 5 /* Factor such that when multiplied with
        the approx threshold will allow partitioning of the reached disjunct. */
#define FSM_HD_FROM 0 /* Constant to indicate to subsetting method whether to use
                       * specified method or remap_ua */
#define FSM_HD_DEADEND 1 /* Constant to indicate to subsetting method whether to use
                       * specified method or remap_ua */
#define FSM_HD_MIN_SIZE_FROM 700 /* threshold to indicate not to subset if under this size */
#define FSM_HD_DONT_FREE 0 /* flag to indicate that mdds should not be freed */
#define FSM_HD_FREE 1 /* flag to indicate that mdds should be freed */
#define FSM_HD_SP_THRESHOLD 2000 /* Threshold for short_paths method */

/*---------------------------------------------------------------------------*/
/* Type declarations                                                         */
/*---------------------------------------------------------------------------*/
typedef struct FsmHdSizeStruct FsmHdSizeStruct_t;

/*---------------------------------------------------------------------------*/
/* Structure declarations                                                    */
/*---------------------------------------------------------------------------*/


struct FsmHdSizeStruct {
  mdd_t *bdd;
  int size;
  long rnd;
};
                  

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

static mdd_t * ComputeNewSeedsAtDeadEnd(mdd_manager *mddManager, Img_ImageInfo_t *imageInfo, mdd_t *reachableStates, FsmHdStruct_t *hdInfo, int greedy, int verbosity);
static mdd_t * ComputeImageOfDecomposedParts(mdd_manager *mddManager, Img_ImageInfo_t *imageInfo, mdd_t *reachedSet, mdd_t *reachableStates, FsmHdStruct_t *hdInfo, int frontierApproxThreshold, int nvars, int *failedPartition, int greedy, int verbosity);
static mdd_t * ProcessDisjunctsRecursive(mdd_manager *mddManager, Img_ImageInfo_t *imageInfo, mdd_t **reachedSet, mdd_t *reachableStates, FsmHdStruct_t *hdInfo, int frontierApproxThreshold, int nvars, int *failedPartition, int reachedSize, int greedy, int verbosity);
static mdd_t * ComputeSubsetBasedOnMethod(mdd_t *fromBetween, mdd_t *fromLowerBound, FsmHdTravOptions_t *options, int fromOrDeadEnd);
static mdd_t * AddStates(mdd_t *a, mdd_t *b, int freeA, int freeB);
static mdd_t * SubtractStates(mdd_t *a, mdd_t *b, int freeA, int freeB);
static int RandomCompare(const void *ptrX, const void *ptrY);

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

FsmHdTravOptions_t *
Fsm_FsmHdGetTravOptions(int nvars)
{
  char *frontierApproxMethodString, *frontierApproxThresholdString;
  char *scrapStatesString;
  char *deadEndString, *qualityString;
  char *newOnlyString, *partialImageString, *deadEndSubsetMethodString;
  FsmHdTravOptions_t *options;

  if (nvars < 0) return NIL(FsmHdTravOptions_t);
  options = ALLOC(FsmHdTravOptions_t, 1);
  if (options == NIL(FsmHdTravOptions_t)) {
    return NIL(FsmHdTravOptions_t);
  }

  options->nvars = nvars;
  /* set defaults */
  options->frontierApproxMethod = BDD_APPROX_RUA;
  options->frontierApproxThreshold = FSM_HD_FRONTIER_APPROX_THRESHOLD;
  options->quality = 1.0;
  options->deadEnd = Fsm_Hd_DE_Con_c;
  options->deadEndSubsetMethod = BDD_APPROX_RUA;
  options->scrapStates = TRUE;
  options->newOnly = FALSE;
  options->onlyPartialImage = FALSE;
  options->qualityBias = 0.5;
  
  /*  method for subsetting reached, from */
  frontierApproxMethodString = Cmd_FlagReadByName("hd_frontier_approx_method");
  if (frontierApproxMethodString != NIL(char)) {
    if (strcmp(frontierApproxMethodString, "heavy_branch") == 0) {
      options->frontierApproxMethod = BDD_APPROX_HB;
    } else if (strcmp(frontierApproxMethodString, "short_paths") == 0) {
      options->frontierApproxMethod = BDD_APPROX_SP;
    } else if (strcmp(frontierApproxMethodString, "under_approx") == 0) {
      options->frontierApproxMethod = BDD_APPROX_UA;
    } else if (strcmp(frontierApproxMethodString, "remap_ua") == 0) {
      options->frontierApproxMethod = BDD_APPROX_RUA;
    } else if (strcmp(frontierApproxMethodString, "compress") == 0) {
      options->frontierApproxMethod = BDD_APPROX_COMP;
    } else if (strcmp(frontierApproxMethodString, "biased_rua") == 0) {
      options->frontierApproxMethod = BDD_APPROX_BIASED_RUA;
    }
  }

  /* threshold value */
  frontierApproxThresholdString = Cmd_FlagReadByName("hd_frontier_approx_threshold");
  if (frontierApproxThresholdString != NIL(char)) {
    options->frontierApproxThreshold = strtol(frontierApproxThresholdString,
                                           NULL, 10);
  }
  /* adjust threshold if too small */
  if ((options->frontierApproxMethod == BDD_APPROX_HB) ||
      (options->frontierApproxMethod == BDD_APPROX_SP) ||
      (options->frontierApproxMethod == BDD_APPROX_COMP)) {
    if (options->frontierApproxThreshold < FSM_HD_SP_THRESHOLD) {
      fprintf(vis_stdout, "** hd warning: Threshold too low, Correcting Frontier Approx Threshold to %d\n", FSM_HD_SP_THRESHOLD);
      options->frontierApproxThreshold =  FSM_HD_SP_THRESHOLD;
    }
  } else if (options->frontierApproxThreshold < 0) {
    fprintf(vis_stdout, "** hd warning: Threshold negative, Correcting Frontier Approx Threshold to 0\n");
    options->frontierApproxThreshold = FSM_HD_FRONTIER_APPROX_THRESHOLD;
  }


  /* quality option for remap_ua */
  qualityString = Cmd_FlagReadByName("hd_approx_quality");
  if (qualityString != NIL(char)) {
    options->quality = strtod(qualityString, NULL);
    if (options->quality < 0.0) {
      options->quality = 1.0;
    }
  }

  /* quality option for biased_rua */
  qualityString = Cmd_FlagReadByName("hd_approx_bias_quality");
  if (qualityString != NIL(char)) {
    options->qualityBias = strtod(qualityString, NULL);
    if (options->qualityBias < 0.0) {
      options->qualityBias = 0.5;
    }
  }

  /* method for dead end computations */
  deadEndString = Cmd_FlagReadByName("hd_dead_end");
  if (deadEndString != NIL(char)) {
    if (strcmp(deadEndString, "brute_force") == 0) {
      options->deadEnd = Fsm_Hd_DE_BF_c;
    } else if (strcmp(deadEndString, "conjuncts") == 0) {
      options->deadEnd = Fsm_Hd_DE_Con_c;
    } else if (strcmp(deadEndString, "hybrid") == 0) {
      options->deadEnd = Fsm_Hd_DE_Hyb_c;
    }
  }
  /*  method for subsetting reached, from during deadEnds*/
  deadEndSubsetMethodString = Cmd_FlagReadByName("hd_dead_end_approx_method");
  if (deadEndSubsetMethodString != NIL(char)) {
    if (strcmp(deadEndSubsetMethodString, "heavy_branch") == 0) {
      options->deadEndSubsetMethod = BDD_APPROX_HB;
    } else if (strcmp(deadEndSubsetMethodString, "short_paths") == 0) {
      options->deadEndSubsetMethod = BDD_APPROX_SP;
    } else if (strcmp(deadEndSubsetMethodString, "under_approx") == 0) {
      options->deadEndSubsetMethod = BDD_APPROX_UA;
    } else if (strcmp(deadEndSubsetMethodString, "remap_ua") == 0) {
      options->deadEndSubsetMethod = BDD_APPROX_RUA;
    } else if (strcmp(deadEndSubsetMethodString, "compress") == 0) {
      options->deadEndSubsetMethod = BDD_APPROX_COMP;
    }
  }
  /* option to not add blocking states */
  scrapStatesString = Cmd_FlagReadByName("hd_no_scrap_states");
  if (scrapStatesString != NIL(char)) {
    options->scrapStates = FALSE;
  }

  /* option to consider only new states instead of fromBetween */
  newOnlyString = Cmd_FlagReadByName("hd_new_only");
  if (newOnlyString != NIL(char)) {
    options->newOnly = TRUE;
  }
    
  /* perform approximate traversal with only partial images */
  partialImageString = Cmd_FlagReadByName("hd_only_partial_image");
  if (partialImageString != NIL(char)) {
    options->onlyPartialImage = TRUE;
  }

  return (options);
} /* end of Fsm_FsmHdGetTraversalOptions */

void
Fsm_FsmHdFreeTravOptions(FsmHdTravOptions_t *options)
{
  if (options) FREE(options);
  return;
} /* end of Fsm_FsmHdFreeTravOptions */

int
Fsm_FsmHdOptionReadNumVars(FsmHdTravOptions_t *options)
{
  if (options) return(options->nvars);
  else return 0;
}

int 
Fsm_FsmHdOptionReadFrontierApproxThreshold(FsmHdTravOptions_t *options)
{
  if (options) return (options->frontierApproxThreshold);
  else return 0;
}

boolean
Fsm_FsmHdOptionReadOnlyPartialImage(FsmHdTravOptions_t *options)
{
  if (options) return (options->onlyPartialImage);
  else return FALSE;
}

boolean
Fsm_FsmHdOptionReadScrapStates(FsmHdTravOptions_t *options)
{
  if (options) return (options->scrapStates);
  else return TRUE;
}

double
Fsm_FsmHdOptionReadQuality(FsmHdTravOptions_t *options)
{
  if (options) return(options->quality);
  else return 1.0;
}

bdd_approx_type_t
Fsm_FsmHdOptionReadFrontierApproxMethod(FsmHdTravOptions_t *options)
{
  if (options) return(options->frontierApproxMethod);
  else return BDD_APPROX_RUA;
}

Fsm_HdDeadEndType_t
Fsm_FsmHdOptionReadDeadEnd(FsmHdTravOptions_t *options)
{
  if (options) return(options->deadEnd);
  else return Fsm_Hd_DE_Con_c;
}

boolean
Fsm_FsmHdOptionReadNewOnly(FsmHdTravOptions_t *options)
{
  if (options) return(options->newOnly);
  else return FALSE;
}

bdd_approx_type_t
Fsm_FsmHdOptionReadDeadEndSubsetMethod(FsmHdTravOptions_t *options)
{
  if (options) return(options->deadEndSubsetMethod);
  else return BDD_APPROX_RUA;
}

int 
Fsm_FsmHdOptionSetFrontierApproxThreshold(FsmHdTravOptions_t *options,
                                          int threshold)
{
  if ((options->frontierApproxMethod == BDD_APPROX_HB) ||
      (options->frontierApproxMethod == BDD_APPROX_SP) ||
      (options->frontierApproxMethod == BDD_APPROX_COMP)) {
    if  (threshold < FSM_HD_SP_THRESHOLD) {
      fprintf(vis_stdout, "** hd warning: Threshold too low, Correcting Frontier Approx Threshold to %d\n", FSM_HD_SP_THRESHOLD);
      options->frontierApproxThreshold =  FSM_HD_SP_THRESHOLD;
      return 0;
    } else {
      options->frontierApproxThreshold = threshold;
      return 1;
    }
  } else if (threshold < 0) {
    fprintf(vis_stdout, "** hd warning: Threshold negative, Correcting Frontier Approx Threshold to 0\n");
    options->frontierApproxThreshold = FSM_HD_FRONTIER_APPROX_THRESHOLD;
    return 0;
  } else {
    options->frontierApproxThreshold = threshold;
    return 1;
  }
}

void
Fsm_FsmHdOptionSetOnlyPartialImage(FsmHdTravOptions_t *options,
                                   boolean onlyPartial)
{
  if (options) options->onlyPartialImage = onlyPartial;
  return;
}

void
Fsm_FsmHdOptionSetScrapStates(FsmHdTravOptions_t *options, boolean scrapStates)
{
  if (options) options->scrapStates = scrapStates;
  return;
}

void
Fsm_FsmHdOptionSetQuality(FsmHdTravOptions_t *options, double quality)
{
  if (quality > 0.0) 
    if (options) options->quality = quality;
  return;
}

void
Fsm_FsmHdOptionSetFrontierApproxMethod(FsmHdTravOptions_t *options,
                                       bdd_approx_type_t approxMethod)
{
  if (options) options->frontierApproxMethod = approxMethod;
  if ((options->frontierApproxMethod == BDD_APPROX_HB) ||
      (options->frontierApproxMethod == BDD_APPROX_SP) ||
      (options->frontierApproxMethod == BDD_APPROX_COMP)) {
    if  (options->frontierApproxThreshold < FSM_HD_SP_THRESHOLD) {
      fprintf(vis_stdout, "** hd warning: Threshold too low, Correcting Frontier Approx Threshold to %d\n", FSM_HD_SP_THRESHOLD);
      options->frontierApproxThreshold =  FSM_HD_SP_THRESHOLD;
      return;
    }
  } else if (options->frontierApproxThreshold < 0) {
    fprintf(vis_stdout, "** hd warning: Threshold negative, Correcting Frontier Approx Threshold to 0\n");
    options->frontierApproxThreshold = FSM_HD_FRONTIER_APPROX_THRESHOLD;
    return;
  }

  return;
}

void
Fsm_FsmHdOptionSetDeadEnd(FsmHdTravOptions_t *options,
                          Fsm_HdDeadEndType_t deadEnd)
{
  if (options) options->deadEnd = deadEnd;
  return;
}

void
Fsm_FsmHdOptionSetNewOnly(FsmHdTravOptions_t *options, boolean newOnly)
{
  if (options) options->newOnly = newOnly;
  return;
}

void
Fsm_FsmHdOptionSetDeadEndSubsetMethod(FsmHdTravOptions_t *options,
                                      bdd_approx_type_t approxMethod)
{
  if (options) options->deadEndSubsetMethod = approxMethod;
  return;
}


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

FsmHdMintSizeStats_t *
FsmHdStatsStructAlloc(void)
{
  FsmHdMintSizeStats_t *stats;

  stats = ALLOC(FsmHdMintSizeStats_t, 1);
  if (stats == NULL) return NULL;

  stats->sizeTo = 0;
  stats->mintermTo = 0.0;
  stats->sizeFrom = 0;
  stats->mintermFrom = 0.0;
  stats->sizeFromSubset = 0;
  stats->mintermFromSubset = 0.0;
  stats->sizeReached = 0;
  stats->mintermReached = 0.0;
  
  return stats;
} /* end of FsmHdStatsStructAlloc */


void
FsmHdStatsStructFree(FsmHdMintSizeStats_t *stats)
{
  FREE(stats);
  return;
} /* end of FsmHdStatsStructFree */



void
FsmHdStatsComputeSizeAndMinterms(mdd_manager *mddManager,
                       mdd_t *f,
                       array_t *varArray,
                       int nvars,
                       Fsm_Hd_Quant_t field,
                       FsmHdMintSizeStats_t *stats)
{
  int size;
  double minterms = 0.0;
  
  size = mdd_size(f);
  if (nvars <= 1023)
    minterms = mdd_count_onset(mddManager, f, varArray);
  
  switch (field) {
  case Fsm_Hd_From: {
    stats->sizeFrom = size;
    stats->mintermFrom = minterms; 
    break;
  }
  case Fsm_Hd_To: {
    stats->sizeTo = size;
    stats->mintermTo = minterms;
    break;
  }
  case Fsm_Hd_FromSubset: {
    stats->sizeFromSubset = size;
    stats->mintermFromSubset = minterms;
    break;
  }
  case Fsm_Hd_Reached: {
    stats->sizeReached = size;
    stats->mintermReached = minterms;
    break;
  }
  default: break;
  }
  return;
} /* end of FsmHdStatsComputeSizeAndMinterms */

void
FsmHdStatsReset(FsmHdMintSizeStats_t *stats,
                Fsm_Hd_Quant_t field)
{

  switch (field) {
  case Fsm_Hd_To:
    stats->sizeTo = 0;
    stats->mintermTo = 0.0;
    break;
  case Fsm_Hd_From:
    stats->sizeFrom = 0;
    stats->mintermFrom = 0.0;
    break;
  case Fsm_Hd_FromSubset:
    stats->sizeFromSubset = 0;
    stats->mintermFromSubset = 0.0;
    break;
  case Fsm_Hd_Reached:
    stats->sizeReached = 0;
    stats->mintermReached = 0.0;
    break;
  }
  return;
} /* end of FsmHdStatsReset */


void
FsmHdPrintOptions(void)
{
  FsmHdTravOptions_t *options;
  char *dummy;
  options = Fsm_FsmHdGetTravOptions( 0);
  if (options == NIL(FsmHdTravOptions_t)) {
    fprintf(vis_stderr, "HD: Unalble to get options\n");
    return;
  }

  dummy = ALLOC(char, 20);
  if (dummy == NULL) return;
  switch (options->frontierApproxMethod) {
  case BDD_APPROX_HB: sprintf(dummy, "%s", "heavy_branch"); break;
  case BDD_APPROX_SP: sprintf(dummy, "%s", "short_paths"); break;

  case BDD_APPROX_UA: sprintf(dummy, "%s", "under_approx"); break;
  case BDD_APPROX_RUA: sprintf(dummy, "%s", "remap_ua"); break;
  case BDD_APPROX_BIASED_RUA: sprintf(dummy, "%s", "biased_rua"); break;

  case BDD_APPROX_COMP: sprintf(dummy, "%s", "compress"); break;
  default: sprintf(dummy, "%s", "remap_ua"); break;
  }
    
  fprintf(vis_stdout, "HD: Approx method for frontier subsets = %s\n", dummy);
  fprintf(vis_stdout, "HD: Threshold for frontier approximation = %d\n", options->frontierApproxThreshold);
  fprintf(vis_stdout, "HD: Quality option for approx methods = %g\n", options->quality);
  fprintf(vis_stdout, "HD: Bias quality option for the biased approx method = %g\n", options->qualityBias);

  switch (options->deadEnd) {
  case Fsm_Hd_DE_BF_c: sprintf(dummy, "%s", "brute_force"); break;
  case Fsm_Hd_DE_Con_c: sprintf(dummy, "%s", "conjuncts"); break;
  case Fsm_Hd_DE_Hyb_c: sprintf(dummy, "%s", "hybrid"); break;
  }
    
  fprintf(vis_stdout, "HD: Dead-End method = %s\n", dummy);

  switch (options->deadEndSubsetMethod) {
  case BDD_APPROX_HB: sprintf(dummy, "%s", "heavy_branch"); break;
  case BDD_APPROX_SP: sprintf(dummy, "%s", "short_paths"); break;

  case BDD_APPROX_UA: sprintf(dummy, "%s", "under_approx"); break;
  case BDD_APPROX_RUA: sprintf(dummy, "%s", "remap_ua"); break;

  case BDD_APPROX_COMP: sprintf(dummy, "%s", "compress"); break;
  default: sprintf(dummy, "%s", "remap_ua"); break;
  }
    
  fprintf(vis_stdout, "HD: Dead-end approx method = %s\n", dummy);
  FREE(dummy);
  
  fprintf(vis_stdout, "HD: Add scrap states option = ");
  if (options->scrapStates == TRUE) fprintf(vis_stdout, "yes\n");
  else fprintf(vis_stdout, "no\n");

  fprintf(vis_stdout, "HD: Only new states in frontier option = ");
  if (options->newOnly == TRUE) fprintf(vis_stdout, "yes\n");
  else fprintf(vis_stdout, "no\n");

  fprintf(vis_stdout, "HD: Only partial image option = ");
  if (options->onlyPartialImage == TRUE) fprintf(vis_stdout, "yes\n");
  else fprintf(vis_stdout, "no\n");

  Fsm_FsmHdFreeTravOptions(options);
  Img_ImagePrintPartialImageOptions();
  
  fprintf(vis_stdout, "See help page on print_hd_options for setting these options\n");
  return;

} /* end of FsmHdPrintOptions */

mdd_t *
FsmHdDeadEnd(mdd_manager *mddManager,
             Img_ImageInfo_t *imageInfo,
             mdd_t *reachableStates,
             FsmHdStruct_t *hdInfo, 
             array_t *varArray,
             int greedy,
             int verbosity)
{
  mdd_t *fromLowerBound;
  
  /* if residue is around to long, throw it away */
  if (hdInfo->residueCount >= 4) {
    if (hdInfo->deadEndResidue != NULL) {
      mdd_free(hdInfo->deadEndResidue);
      hdInfo->deadEndResidue = NULL;
    }
  }
      
  if ((hdInfo->deadEndResidue != NULL) && (greedy)) {
    /* dont compute new seeds out of dead end
     * if some states are left over from previous iterations.
     */
    fromLowerBound = hdInfo->deadEndResidue;
    hdInfo->deadEndResidue = NULL;
    assert(!mdd_is_tautology(fromLowerBound, 0));
    if (hdInfo->interiorStates != NIL(mdd_t))
      assert(mdd_lequal(fromLowerBound, hdInfo->interiorStates, 1, 0));
    if (verbosity >= 2) fprintf(vis_stdout, "Residue had new seeds\n");
    (hdInfo->residueCount)++;
  } else {
    /* compute new seeds by decomposing reachable states */
    if (hdInfo->deadEndResidue != NULL) {
      mdd_free(hdInfo->deadEndResidue);
      hdInfo->deadEndResidue = NULL;
    }
    hdInfo->residueCount = 0;
    fromLowerBound = ComputeNewSeedsAtDeadEnd(mddManager,
                                              imageInfo,
                                              reachableStates,
                                              hdInfo,
                                              greedy,
                                              verbosity);

  }
            
  /* mark that this is a dead end */
  hdInfo->deadEndComputation = TRUE;

  return (fromLowerBound);
} /* end of FsmHdDeadEnd */

void
FsmHdFromComputeDenseSubset(mdd_manager *mddManager,
                            mdd_t **fromLowerBound,
                            mdd_t **fromUpperBound,
                            mdd_t **reachableStates,
                            mdd_t **image,
                            mdd_t *initialStates,
                            FsmHdStruct_t *hdInfo,
                            int shellFlag,
                            array_t *onionRings,
                            int sizeOfOnionRings,
                            array_t *varArray)
{

  int thisTimeSubset;
  mdd_t *fromSubset, *fromBetween;
  mdd_t *temp, *temp1, *nonInteriorStates;
  float densityFrom = 0.0, densityFromSubset = 0.0; /* initialize to pacify */
  int numReorders;
  double mintResidue;
  int sizeResidue;
  int frontierApproxThreshold;
  int nvars;
  
  numReorders = bdd_read_reorderings(mddManager);
  frontierApproxThreshold = Fsm_FsmHdOptionReadFrontierApproxThreshold(hdInfo->options);
  nvars =  Fsm_FsmHdOptionReadNumVars(hdInfo->options);

   /* if image is bigger and smaller in size, switch. Add initial states.
   * Reset dead states since they are no longer valid.
   */

#if 0 /* taken out for counterexamples  */
  if ((FsmHdStatsReadMintermReached(hdInfo->stats) <
       FsmHdStatsReadMintermTo(hdInfo->stats)) &&
      (FsmHdStatsReadSizeReached(hdInfo->stats) > FsmHdStatsReadSizeTo(hdInfo->stats))) {
    if (hdInfo->interiorStates != NIL(mdd_t)) mdd_free(hdInfo->interiorStates);
    hdInfo->interiorStates = bdd_zero(mddManager);
    bdd_free(*reachableStates);
    mdd_free(*fromUpperBound);
    *reachableStates = mdd_or(initialStates, *image, 1, 1);
    *reachableStates = AddStates(*fromLowerBound, *reachableStates,
                                 FSM_HD_DONT_FREE, FSM_HD_FREE);
    FsmHdStatsComputeSizeAndMinterms(mddManager, *reachableStates,
                             varArray, nvars, Fsm_Hd_Reached, hdInfo->stats);
    *fromUpperBound = mdd_dup(*reachableStates);
    if (shellFlag && *onionRings) {
      int onionRingCount;
      mdd_t *onionRing;
      fprintf(vis_stdout, "Onion Rings no longer valid. \n");
      arrayForEachItem(mdd_t *, *onionRings, onionRingCount, onionRing) {
        mdd_free(onionRing);
      }
      *onionRings = NIL(array_t);
    }
  }
#endif  

  /* states which may produce new states as successors */
  nonInteriorStates = SubtractStates(*reachableStates, hdInfo->interiorStates,
                                     FSM_HD_DONT_FREE, FSM_HD_DONT_FREE);

  /* try and reinject the residue from a dead end into the frontier */
  if (hdInfo->deadEndResidue != NULL) {
    temp = mdd_or(*fromLowerBound, hdInfo->deadEndResidue, 1, 1);
  } else {
    temp = mdd_dup(*fromLowerBound);
  }
  /* calculate the actual "from" going to be used */
  if ( Fsm_FsmHdOptionReadNewOnly(hdInfo->options) == FALSE) {
    fromBetween = bdd_squeeze(temp, nonInteriorStates);
    mdd_free(temp);
  } else {
    fromBetween = temp;
  }
  mdd_free(nonInteriorStates);

  if (nvars <= 1023) {
    FsmHdStatsComputeSizeAndMinterms(mddManager, fromBetween,
                                     varArray, nvars,
                                     Fsm_Hd_From, hdInfo->stats);
  }
    
  /* if from is bigger in size than image, use image as from*/
  if (( Fsm_FsmHdOptionReadNewOnly(hdInfo->options) == FALSE) && (*image != NULL) &&
      (FsmHdStatsReadSizeTo(hdInfo->stats) < FsmHdStatsReadSizeFrom(hdInfo->stats))) {
    mdd_free(fromBetween);
    FsmHdStatsSetSizeFrom(hdInfo->stats, FsmHdStatsReadSizeTo(hdInfo->stats));
    FsmHdStatsSetMintermFrom(hdInfo->stats, FsmHdStatsReadMintermTo(hdInfo->stats));
    fromBetween = *image;
    *image = NULL;
  }

  /* free image since no longer required */
  if (*image != NULL) mdd_free(*image);

  assert(FsmHdStatsReadMintermFromSubset(hdInfo->stats) == 0.0);
  assert(FsmHdStatsReadSizeFromSubset(hdInfo->stats) == 0);

  /* flag to indicate whether a subset is computed this time */
  thisTimeSubset = 0;

  /* if size of from exceeds threshold, subset */
  if ((FsmHdStatsReadSizeFrom(hdInfo->stats) >  frontierApproxThreshold) &&
      (FsmHdStatsReadSizeFrom(hdInfo->stats) > FSM_HD_MIN_SIZE_FROM)) {

    /* get subset based on the different methods */
    fromSubset = ComputeSubsetBasedOnMethod(fromBetween,
                                            *fromLowerBound, hdInfo->options,
                                            FSM_HD_FROM);
    if (nvars <= 1023) {
      /* record stats */
      FsmHdStatsComputeSizeAndMinterms(mddManager, fromSubset,
                             varArray, nvars, 
                           Fsm_Hd_FromSubset, hdInfo->stats);
      /* compute density */
      densityFromSubset = FsmHdStatsReadDensityFromSubset(hdInfo->stats);
      densityFrom = FsmHdStatsReadDensityFrom(hdInfo->stats);
    }
                
    /* check density of subset, discard the subset if not required*/
    if (((nvars <= 1023) &&
         (((densityFrom < densityFromSubset) ||
            (FsmHdStatsReadSizeFrom(hdInfo->stats) > 2*frontierApproxThreshold)))) ||
        ((nvars > 1023) &&
         ((bdd_apa_compare_ratios(nvars, fromSubset, fromBetween,
                                  FsmHdStatsReadSizeFromSubset(hdInfo->stats),
                                  FsmHdStatsReadSizeFrom(hdInfo->stats)) ||
            (FsmHdStatsReadSizeFrom(hdInfo->stats) > 2*frontierApproxThreshold))))) {

      /* subset chosen this time */
      thisTimeSubset = 1;
      
      mdd_free(fromBetween);

      /* discard the old bounds only if not required later */
      if (hdInfo->deadEndComputation || hdInfo->firstSubset) {
        /* store the residue since subset has been taken. Use them in
           subsequent iterations*/
        hdInfo->deadEndResidue = SubtractStates(*fromLowerBound, fromSubset,
                                         FSM_HD_DONT_FREE, FSM_HD_DONT_FREE);
        if(mdd_is_tautology(hdInfo->deadEndResidue, 0)) {
          mdd_free(hdInfo->deadEndResidue);
          hdInfo->deadEndResidue = NULL;
        } else if (mdd_size(hdInfo->deadEndResidue) > FSM_HD_DEADEND_RESIDUE_LIMIT) {
          mdd_free(hdInfo->deadEndResidue);
          hdInfo->deadEndResidue = NULL;
        } else {
          hdInfo->firstSubset = 0;
        }
                
      }
      
      /* set from lower bound for next image computation */
      mdd_free(*fromLowerBound);
      *fromLowerBound = fromSubset;
                    
    }  else {
      /* discard the old bounds; since no subset is discarded */
      mdd_free(*fromLowerBound);
      /* Undo subset computation */
      mdd_free(fromSubset);
      *fromLowerBound = fromBetween;
    }
  } else {
    /* No subsetting here */
    mdd_free(*fromLowerBound);
    *fromLowerBound = fromBetween;
  }

  /* remove frontier from dead-end residue */
  if ((hdInfo->deadEndResidue != NULL) && (!hdInfo->deadEndComputation)) {
    hdInfo->deadEndResidue = SubtractStates(hdInfo->deadEndResidue, *fromLowerBound,
                                     FSM_HD_FREE, FSM_HD_DONT_FREE);
    if (mdd_is_tautology(hdInfo->deadEndResidue, 0)) {
      mdd_free(hdInfo->deadEndResidue);
      hdInfo->deadEndResidue = NULL;
    } else if (mdd_size(hdInfo->deadEndResidue) > FSM_HD_DISJ_SIZE) {
      mdd_free(hdInfo->deadEndResidue);
      hdInfo->deadEndResidue = NULL;
    }
  }

  /*
   * check size of reached if there were any reorderings since
   * the last time the size was measured
   */
  if (bdd_read_reorderings(mddManager)> numReorders) {
    FsmHdStatsSetSizeReached(hdInfo->stats, mdd_size(*reachableStates));
  }
  /* Modify Reachable states if necessary. The conditions are if
   * blocking states are not to be included in Reached. The exceptions
   * are if no subset was computed or the size of Reached is less than
   * 30000 or if this were a DEAD-END computation (very important:
   * this rule isnt violated). That is ALL states computed out of a
   * DEAD_END must ALWAYS be added to reached.
   */
  if (thisTimeSubset && (Fsm_FsmHdOptionReadScrapStates(hdInfo->options) == FALSE) &&
      (!hdInfo->deadEndComputation) && 
      (FsmHdStatsReadSizeReached(hdInfo->stats) > 30000)) {
    temp = AddStates(*fromUpperBound, *fromLowerBound,
                     FSM_HD_FREE, FSM_HD_DONT_FREE);
    temp1 = bdd_squeeze(temp, *reachableStates);
    mdd_free(*reachableStates);
    mdd_free(temp);
    *reachableStates = temp1;
    FsmHdStatsComputeSizeAndMinterms(mddManager, *reachableStates, varArray,
                           nvars, Fsm_Hd_Reached, hdInfo->stats);
    if (shellFlag && onionRings) { /* restrict the onion rings to the modified reachable set */
      int onionRingCount;
      for (onionRingCount = sizeOfOnionRings; onionRingCount < array_n(onionRings);
           onionRingCount++) {
        mdd_t *onionRing;
        onionRing = array_fetch(mdd_t *, onionRings, onionRingCount);
        temp = mdd_and(onionRing, *reachableStates, 1, 1);
        mdd_free(onionRing);
        array_insert(mdd_t *, onionRings, onionRingCount, temp);
      }
    }

  } else {
    mdd_free(*fromUpperBound);
  }

  /* Throw away dead-end residue when it gets to be a pain.
   * i.e., when its density is larger than Reached or when the
   * fraction of minterms retained in it are too small.
   */
  if (hdInfo->deadEndResidue != NULL) {
    mintResidue = mdd_count_onset(mddManager, hdInfo->deadEndResidue, varArray);
    sizeResidue = mdd_size(hdInfo->deadEndResidue);
    if ((mintResidue <  0.001*FsmHdStatsReadMintermReached(hdInfo->stats)) ||
        (mintResidue/sizeResidue > FsmHdStatsReadDensityReached(hdInfo->stats))) {
      mdd_free(hdInfo->deadEndResidue);
      hdInfo->deadEndResidue = NULL;
    }
  }

  /* set fromUpperBound for the next iteration */
  *fromUpperBound = *fromLowerBound;
  return;
} /* end of FsmHdFromComputeDenseSubset */

 
/*---------------------------------------------------------------------------*/
/* Definition of static functions                                          */
/*---------------------------------------------------------------------------*/

static mdd_t *
ComputeNewSeedsAtDeadEnd(mdd_manager *mddManager,
                         Img_ImageInfo_t *imageInfo,
                         mdd_t *reachableStates,
                         FsmHdStruct_t *hdInfo, 
                         int greedy,
                         int verbosity)
{
  mdd_t *toCareSet, *newSeeds = NULL, *possibleSeeds, *temp;
  mdd_t *reachedSubset, *reachedSet;
  int failedPartition = 0;
  int frontierApproxThreshold =  Fsm_FsmHdOptionReadFrontierApproxThreshold(hdInfo->options);
  int deadEnd =  Fsm_FsmHdOptionReadDeadEnd(hdInfo->options);
  int nvars =  Fsm_FsmHdOptionReadNumVars(hdInfo->options);
  
  /* in brute force method, compute the image of the entire reached */
  if (deadEnd == Fsm_Hd_DE_BF_c) {
    toCareSet = mdd_not(reachableStates);
    possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                          reachableStates,
                                                          reachableStates,
                                                          toCareSet);
    mdd_free(toCareSet);
    newSeeds = SubtractStates(possibleSeeds, reachableStates,
                              FSM_HD_FREE, FSM_HD_DONT_FREE);
    return (newSeeds);
  }

  /* in HYBRID method, first try with  a subset of the reached set */
  reachedSet = SubtractStates(reachableStates, hdInfo->interiorStates,
                            FSM_HD_DONT_FREE, FSM_HD_DONT_FREE);
  if (deadEnd ==  Fsm_Hd_DE_Hyb_c) {
    /* compute image of a subset of reached */
    reachedSubset = ComputeSubsetBasedOnMethod(reachedSet,
                                               NIL(mdd_t),
                                               hdInfo->options,
                                               FSM_HD_DEADEND);
                                               
    if (bdd_size(reachedSubset) < 2*frontierApproxThreshold) {
      toCareSet = mdd_not(reachableStates);
      possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                            reachedSubset,
                                                            reachableStates,
                                                            toCareSet);
      mdd_free(toCareSet);
      /* compute new states */
      newSeeds = SubtractStates(possibleSeeds, reachableStates,
                              FSM_HD_FREE, FSM_HD_DONT_FREE);
      /* add to interior states */
      hdInfo->interiorStates = AddStates(hdInfo->interiorStates, reachedSubset,
                                  FSM_HD_FREE, FSM_HD_DONT_FREE);
      if (verbosity >= 2) {
        fprintf(vis_stdout, "Interior states added, Size = %d\n",
                mdd_size(hdInfo->interiorStates));
      }

      if ((!mdd_is_tautology(newSeeds, 0)) && (greedy))  {
        mdd_free(reachedSet);
        mdd_free(reachedSubset);
        return (newSeeds);
      } else {
        /* reset newSeeds value for further computations */
        if (mdd_is_tautology(newSeeds, 0)) {
          mdd_free(newSeeds);
          newSeeds = NULL;
        }
        /* if new seeds not found or iter_decomp unsuccessful, set
         * reachedSet to the other part (which may be a part of the
         * reachable states (former) or the reachable states itself
         * (latter)).
         */
        reachedSet = SubtractStates(reachedSet, reachedSubset,
                                      FSM_HD_FREE, FSM_HD_FREE);
      }
    } else {
      mdd_free(reachedSubset);
    }
  }
  
  /* compute the decomposed image of reachedSet */
  if ((deadEnd == Fsm_Hd_DE_Con_c) ||
      (deadEnd == Fsm_Hd_DE_Hyb_c)) {
    /* this flag is set when the image of some decomposed part is
     * not computed due to its size.
     */
    failedPartition = 0;
    /* save any previously computed new states */
    temp = newSeeds;
    /* compute new seeds by recursively decomposing reachedSet */
    newSeeds = ComputeImageOfDecomposedParts(mddManager,
                                             imageInfo,
                                             reachedSet,
                                             reachableStates,
                                             hdInfo,
                                             frontierApproxThreshold,
                                             nvars,
                                             &failedPartition,
                                             greedy,
                                             verbosity);
    /* combine new states */
    newSeeds = AddStates(newSeeds, temp, FSM_HD_FREE, FSM_HD_FREE);

    /* if no remains or partition didnt fail, return new Seeds */
    /* if greedy and new states found, return */
    if ((mdd_lequal(reachedSet, hdInfo->interiorStates, 1, 1)) ||
        ((newSeeds != NULL) && (greedy))) {
      /* assert that the entire image was computed */
      if (mdd_lequal(reachedSet, hdInfo->interiorStates, 1, 1)) assert(!failedPartition);
      mdd_free(reachedSet);
      if (newSeeds == NULL)
        newSeeds = bdd_zero(mddManager);
      else
        assert(!mdd_is_tautology(newSeeds, 0));
      return(newSeeds);
    }

    /* continue if either no new states have been found or the
     * computation is not greedy */
    /* new set whose image is to be computed  */
    reachedSet = SubtractStates(reachedSet, hdInfo->interiorStates,
                                FSM_HD_FREE, FSM_HD_DONT_FREE);
    
    
    /* one last shot with a subsetting the remains */
    /* reachedSubset is the subset of reached remains */
    reachedSubset = ComputeSubsetBasedOnMethod(reachedSet,
                                               NIL(mdd_t),
                                               hdInfo->options,
                                               FSM_HD_DEADEND);
        
    /* compute image of the subset above */
    toCareSet = mdd_not(reachableStates);
    possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                          reachedSubset,
                                                          reachableStates,
                                                          toCareSet);
    mdd_free(toCareSet);
    /* add to interior states */
    hdInfo->interiorStates = AddStates(hdInfo->interiorStates, reachedSubset,
                                FSM_HD_FREE, FSM_HD_FREE);
    if (verbosity >= 2) {
      fprintf(vis_stdout, "Interior states added, Size = %d\n",
              mdd_size(hdInfo->interiorStates));
    }
    /* combine all new states */
    temp = newSeeds;
    newSeeds = SubtractStates(possibleSeeds, reachableStates,
                              FSM_HD_FREE, FSM_HD_DONT_FREE);
    newSeeds = AddStates(newSeeds, temp, FSM_HD_FREE, FSM_HD_FREE);

    /* if no new seeds, compute the image of (reachedSet - reachedSubset) */
    if ((mdd_is_tautology(newSeeds, 0) != 1) && (greedy)) {
      mdd_free(reachedSet);
      mdd_free(reachedSubset);
      return newSeeds;
    } else {
      /* reset newSeeds for further computation */
      if (mdd_is_tautology(newSeeds, 0)) {
        mdd_free(newSeeds);
        newSeeds = NULL;
      }
      /* set reached = reachedSet - reachedSubset */
      reachedSet  = SubtractStates(reachedSet, reachedSubset,
                                   FSM_HD_FREE, FSM_HD_FREE);
    }

    /* save previously computed new states */
    temp = newSeeds;
    failedPartition = 0;
    /* try decomposition again */
    newSeeds = ComputeImageOfDecomposedParts(mddManager,
                                             imageInfo,
                                             reachedSet,
                                             reachableStates, hdInfo, 
                                             frontierApproxThreshold,
                                             nvars,
                                             &failedPartition,
                                             greedy,
                                             verbosity);
    /* combine all new states */
    newSeeds = AddStates(newSeeds, temp, FSM_HD_FREE, FSM_HD_FREE);
    /* if no remains or partition didnt fail, return new Seeds */
    /* if greedy and new states found, return */
    if ((mdd_lequal(reachedSet, hdInfo->interiorStates, 1, 1)) ||
        ((newSeeds != NULL) && (greedy))) {
      /* assert that the entire image was computed */
      if (mdd_lequal(reachedSet, hdInfo->interiorStates, 1, 1)) assert(!failedPartition);
      mdd_free(reachedSet);
      if (newSeeds == NULL)
        newSeeds = bdd_zero(mddManager);
      else
        assert(!mdd_is_tautology(newSeeds, 0));
      return(newSeeds);
    }

    /* continue if either no new states have been found or the
     * computation is not greedy */
    /* new set whose image is to be computed  */
    reachedSet = SubtractStates(reachedSet, hdInfo->interiorStates,
                                FSM_HD_FREE, FSM_HD_DONT_FREE);
    
    /* now bite the bullet and compute the image of the whole
     * remainder.
     */
    toCareSet = mdd_not(reachableStates);
    possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                          reachedSet,
                                                          reachableStates,
                                                          toCareSet);
    mdd_free(toCareSet);
    /* add to interios states */
    hdInfo->interiorStates = AddStates(hdInfo->interiorStates, reachedSet,
                                FSM_HD_FREE, FSM_HD_FREE);
    if (verbosity >= 2) {
      fprintf(vis_stdout, "Interior states added, Size = %d\n",
              mdd_size(hdInfo->interiorStates));
    }
    /* combine all new states */
    temp = newSeeds;
    newSeeds = SubtractStates(possibleSeeds, reachableStates,
                              FSM_HD_FREE, FSM_HD_DONT_FREE);
    newSeeds = AddStates(newSeeds, temp, FSM_HD_FREE, FSM_HD_FREE);

  }
  return newSeeds;
}


static mdd_t *
ComputeImageOfDecomposedParts(mdd_manager *mddManager,
                              Img_ImageInfo_t *imageInfo,
                              mdd_t *reachedSet,
                              mdd_t *reachableStates,
                              FsmHdStruct_t *hdInfo, 
                              int frontierApproxThreshold,
                              int nvars,
                              int *failedPartition,
                              int greedy,
                              int verbosity)
{

  mdd_t *temp, *newSeeds, *possibleSeeds;
  mdd_t *toCareSet;
  int reachedSize, threshold;
    
  /* initialize some variables */
  newSeeds = NULL;
  reachedSize = mdd_size(reachedSet);
  threshold = (FSM_HD_REACHED_THRESHOLD > frontierApproxThreshold) ? FSM_HD_REACHED_THRESHOLD : frontierApproxThreshold;

  /* compute the image of reached if it is small */
  if (reachedSize < threshold) {
    toCareSet = mdd_not(reachableStates);
    possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                          reachedSet,
                                                          reachableStates,
                                                          toCareSet);
    mdd_free(toCareSet);
    newSeeds = SubtractStates(possibleSeeds, reachableStates,
                              FSM_HD_FREE, FSM_HD_DONT_FREE);
    if (mdd_is_tautology(newSeeds, 0)) {
      mdd_free(newSeeds);
      newSeeds = NULL;
    }
    hdInfo->interiorStates = AddStates(hdInfo->interiorStates, reachedSet,
                                FSM_HD_FREE, FSM_HD_DONT_FREE);
    return (newSeeds);
  }

  /* create conjuncts of the complement, disjuncts of the
   * reachable set
   */
  temp = mdd_dup(reachedSet);
  /* set threshold to some non-trivial value */
  if (frontierApproxThreshold <= 0)
    frontierApproxThreshold = FSM_HD_REACHED_THRESHOLD/
      FSM_HD_DEADEND_MAX_SIZE_FACTOR;
  /* compute the image of decomposed parts. Return when some new seeds are
   * found.
   */
  newSeeds = ProcessDisjunctsRecursive(mddManager,
                                       imageInfo,
                                       &temp,
                                       reachableStates,
                                       hdInfo, 
                                       frontierApproxThreshold,
                                       nvars,
                                       failedPartition,
                                       reachedSize,
                                       greedy,
                                       verbosity);
                                         
  return (newSeeds);
    
}

static mdd_t *
ProcessDisjunctsRecursive(mdd_manager *mddManager,
                          Img_ImageInfo_t *imageInfo,
                          mdd_t **reachedSet,
                          mdd_t *reachableStates,
                          FsmHdStruct_t *hdInfo, 
                          int frontierApproxThreshold,
                          int nvars,
                          int *failedPartition,
                          int reachedSize,
                          int greedy,
                          int verbosity)
{
  mdd_t *toCareSet, *newSeeds, *possibleSeeds;
  mdd_t *temp, *subset;
  int  tempSize;
  int hitConstant;
  mdd_t **disjArray;
  char *decompString;
  int i;
  int numParts;
  FsmHdSizeStruct_t sizeStruct;
  array_t *sizeArray;
  
  hitConstant = 0;

  assert((!mdd_is_tautology(*reachedSet, 1)) &&
         (!mdd_is_tautology(*reachedSet, 0)));
  if (verbosity >= 2) {
    fprintf(vis_stdout, "Orig Size = %d, ",mdd_size(*reachedSet));
  }

  /* decompose reached set into disjuntive parts */
  decompString = Cmd_FlagReadByName("hd_decomp");
  if (decompString == NIL(char)) {
    numParts = bdd_var_decomp(*reachedSet,
                                        (bdd_partition_type_t)BDD_DISJUNCTS,
                                        &disjArray);
  } else if (strcmp(decompString, "var") == 0) {
    numParts = bdd_var_decomp(*reachedSet,
                                        (bdd_partition_type_t)BDD_DISJUNCTS,
                                        &disjArray);
  } else if (strcmp(decompString, "gen") == 0) {
    numParts = bdd_gen_decomp(*reachedSet,
                                        (bdd_partition_type_t)BDD_DISJUNCTS,
                                        &disjArray);
  } else if (strcmp(decompString, "approx") == 0) {
    numParts = bdd_approx_decomp(*reachedSet,
                                           (bdd_partition_type_t)BDD_DISJUNCTS,
                                           &disjArray);
  } else if (strcmp(decompString, "iter") == 0) {
    numParts = bdd_iter_decomp(*reachedSet,
                                         (bdd_partition_type_t)BDD_DISJUNCTS,
                                         &disjArray);
  } else {
    numParts = bdd_var_decomp(*reachedSet,
                                        (bdd_partition_type_t)BDD_DISJUNCTS,
                                        &disjArray);
  }
  mdd_free(*reachedSet);

  /* decomposition failed */
  if (numParts == 0) return NULL;

  /* decomposition failed, stop recurring */
  if (numParts == 1) hitConstant = 1;
  
  /* allocate array with bdds and sizes */
  sizeArray = array_alloc(FsmHdSizeStruct_t,  0);  
  if (sizeArray == NIL(array_t)) {
    for (i = 0; i < numParts; i++) {
      mdd_free(disjArray[i]);
    }
    FREE(disjArray);
    return NULL;
  }
  
  /* free all constants and partitions that are larger than the
   * original reached size
   */
  for (i = 0; i < numParts; i++) {
    assert(!mdd_is_tautology(disjArray[i], 1));
    tempSize = bdd_size(disjArray[i]);
    if ((tempSize >= reachedSize) ||
        (tempSize == 1)) {
      if (tempSize == 1) {
        hitConstant = 1;
        assert (mdd_is_tautology(disjArray[i], 0));
      } else {
        *failedPartition = 1;
      }
      mdd_free(disjArray[i]);
    } else {
      sizeStruct.bdd = disjArray[i];
      sizeStruct.size = tempSize;
      sizeStruct.rnd = util_random();
      array_insert_last(FsmHdSizeStruct_t, sizeArray, sizeStruct);
    }
  }
  FREE(disjArray);

  if (array_n(sizeArray) == 0) {
    array_free(sizeArray);
    return (NULL);
  }
  
  /* shuffle parts randomly. */
  if (array_n(sizeArray) > 1) {
    array_sort(sizeArray, RandomCompare);
  } 

  if (verbosity >= 2) {
    arrayForEachItem(FsmHdSizeStruct_t, sizeArray, i, sizeStruct) {
      fprintf(vis_stdout , "disjSize[%d] = %d ", i, sizeStruct.size);
    }
    fprintf(vis_stdout, "\n");
  }
  
  /* now compute image of parts or recur*/
  newSeeds = NULL;
  arrayForEachItem(FsmHdSizeStruct_t, sizeArray, i, sizeStruct) {
    temp = NULL;
    if ((sizeStruct.size <
         frontierApproxThreshold*FSM_HD_DEADEND_MAX_SIZE_FACTOR) ||
        (sizeStruct.size < FSM_HD_DISJ_SIZE)) {
      /* if this partition is small enough, compute its image */
      /* compute image */
      toCareSet = mdd_not(reachableStates);
      possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                            sizeStruct.bdd,
                                                            reachableStates,
                                                            toCareSet); 
      mdd_free(toCareSet);
      /* update interior states */
      hdInfo->interiorStates = AddStates(hdInfo->interiorStates, sizeStruct.bdd,
                                  FSM_HD_FREE, FSM_HD_FREE);
      if (verbosity >= 2) {
        fprintf(vis_stdout, "Interior states added, Size = %d\n",
                mdd_size(hdInfo->interiorStates));
      }
      /* store new states if any were already found */
      temp = newSeeds;
      /* compute new states of this image */
      newSeeds = SubtractStates(possibleSeeds, reachableStates,
                              FSM_HD_FREE, FSM_HD_DONT_FREE);
      /* add up states or set newSeeds to NULL if no new states */
      if (temp != NULL) {
        assert(!mdd_is_tautology(temp, 0));
        newSeeds = AddStates(newSeeds, temp, FSM_HD_FREE, FSM_HD_FREE);
      } else if (mdd_is_tautology(newSeeds, 0) == 1) {
        mdd_free(newSeeds);
        newSeeds = NULL;
      } 
    } else if ((!hitConstant) && ((newSeeds == NULL) ||
                                  (!greedy) ||
                                  (sizeStruct.size < 3* FSM_HD_DISJ_SIZE))) {
      /* recur if size is small enough, or non-greedy computation or
       * no new states have been found yet. But not if a constant is hit. */
      /* store any new states already computed */
      temp = newSeeds;
      newSeeds = ProcessDisjunctsRecursive(mddManager,
                                           imageInfo,
                                           &(sizeStruct.bdd),
                                           reachableStates, hdInfo, 
                                           frontierApproxThreshold,
                                           nvars,
                                           failedPartition,
                                           reachedSize,
                                           greedy,
                                           verbosity);
      /* if new seeds not NULL, then not zero either */
      if (newSeeds != NULL) {
        assert(!mdd_is_tautology(newSeeds, 0));
        /* add states */
        if (temp != NULL) {
          newSeeds = AddStates(temp, newSeeds, FSM_HD_FREE, FSM_HD_FREE);
        }
      } else {
        newSeeds = temp;
      }
    } else {
      /* if this partition is too large, compute the image of its
       *  subset if the subset is small enough */
      /* set failed partition as image is not being computed of this
       *  whole part */
      *failedPartition = 1;
      /* compute subset */
      subset = bdd_approx_remap_ua(sizeStruct.bdd,
                                   (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                   nvars, 0, 1.0);
      mdd_free(sizeStruct.bdd);
      tempSize = bdd_size(subset);
      /* compute image of subset if small enough */
      if ((tempSize < frontierApproxThreshold*FSM_HD_DEADEND_MAX_SIZE_FACTOR) ||
          (tempSize < FSM_HD_DISJ_SIZE)) {
        /* compute image */
        toCareSet = mdd_not(reachableStates);
        possibleSeeds = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                              subset,
                                                              reachableStates,
                                                              toCareSet);
        mdd_free(toCareSet);
        /* update dead states */
        hdInfo->interiorStates = AddStates(hdInfo->interiorStates, subset,
                                    FSM_HD_FREE, FSM_HD_FREE);
        if (verbosity >= 2) {
          fprintf(vis_stdout, "Interior states added, Size = %d\n",
                  mdd_size(hdInfo->interiorStates));
        }
        /* store new states if any were already found */
        temp = newSeeds;
        /* compute new states of this image */
        newSeeds = SubtractStates(possibleSeeds, reachableStates,
                                FSM_HD_FREE, FSM_HD_DONT_FREE);
        /* add up states or set newSeeds to NULL if no new states */
        if (temp != NULL) {
          assert(!mdd_is_tautology(temp, 0));
          newSeeds = AddStates(newSeeds, temp, FSM_HD_FREE, FSM_HD_FREE);
        } else if (mdd_is_tautology(newSeeds, 0) == 1) {
          mdd_free(newSeeds);
          newSeeds = NULL;
        }       
      } else {
        mdd_free(subset);
      }
    }
  } /* end of array for each */

  array_free(sizeArray);
                
  return (newSeeds);

} /* end of ProcessDisjunctsRecursive */



static mdd_t *
ComputeSubsetBasedOnMethod(mdd_t *fromBetween,
                           mdd_t *fromLowerBound,
                           FsmHdTravOptions_t *options,
                           int fromOrDeadEnd)
{
  mdd_t *fromSubset;
  int frontierApproxThreshold =  Fsm_FsmHdOptionReadFrontierApproxThreshold(options);
  int frontierApproxMethod;
  double quality =  Fsm_FsmHdOptionReadQuality(options);
  int nvars =  Fsm_FsmHdOptionReadNumVars(options);
                             
  
  frontierApproxMethod = (fromOrDeadEnd == FSM_HD_FROM) ?
     Fsm_FsmHdOptionReadFrontierApproxMethod(options) :
    ((fromOrDeadEnd == FSM_HD_DEADEND) ?
      Fsm_FsmHdOptionReadDeadEndSubsetMethod(options) :
      Fsm_FsmHdOptionReadFrontierApproxMethod(options));
  
  /* based on approximation method specified, compute subset */
  switch (frontierApproxMethod) {
  case BDD_APPROX_HB:
    frontierApproxThreshold = (frontierApproxThreshold < FSM_HD_SP_THRESHOLD) ? FSM_HD_SP_THRESHOLD: frontierApproxThreshold;
    fromSubset = bdd_approx_hb(fromBetween,
                               (bdd_approx_dir_t)BDD_UNDER_APPROX,
                               nvars, frontierApproxThreshold);
    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_hb(fromLowerBound,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, frontierApproxThreshold);
        }
      }
    }
                    
    break;
  case BDD_APPROX_SP:
    frontierApproxThreshold = (frontierApproxThreshold < FSM_HD_SP_THRESHOLD) ? FSM_HD_SP_THRESHOLD: frontierApproxThreshold;
    fromSubset = bdd_approx_sp(fromBetween,
                               (bdd_approx_dir_t)BDD_UNDER_APPROX,
                               nvars, frontierApproxThreshold,
                               0);
    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_sp(fromLowerBound,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, frontierApproxThreshold,
                                     0);
        }
      }
    }
                    
    break;
  case BDD_APPROX_UA:
    fromSubset = bdd_approx_ua(fromBetween,
                               (bdd_approx_dir_t)BDD_UNDER_APPROX,
                               nvars, frontierApproxThreshold,
                               0, quality);
    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_ua(fromLowerBound,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, frontierApproxThreshold,
                                     0, quality);
        }
      }
    }
    break;
  case BDD_APPROX_RUA:
    fromSubset = bdd_approx_remap_ua(fromBetween,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, frontierApproxThreshold,
                                     quality );
    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_remap_ua(fromLowerBound,
                                           (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                           nvars, frontierApproxThreshold,
                                           quality);
        }
      }
    }
                    
    break;
  case BDD_APPROX_BIASED_RUA:
    if (fromLowerBound == NIL(mdd_t)) {
      fromLowerBound = fromBetween;
    }
    fromSubset = bdd_approx_biased_rua(fromBetween,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                       fromLowerBound,
                                       nvars, frontierApproxThreshold,
                                     quality, options->qualityBias );
    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_remap_ua(fromLowerBound,
                                           (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                           nvars, frontierApproxThreshold,
                                           quality);
        }
      }
    }
    break;
  case BDD_APPROX_COMP:
    frontierApproxThreshold = (frontierApproxThreshold < FSM_HD_SP_THRESHOLD) ? FSM_HD_SP_THRESHOLD: frontierApproxThreshold;
    fromSubset = bdd_approx_compress(fromBetween,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, frontierApproxThreshold);

    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_compress(fromLowerBound,
                                           (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                           nvars, frontierApproxThreshold);
        }
      }
    }

    break;
  default:
    fromSubset = bdd_approx_remap_ua(fromBetween,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, frontierApproxThreshold,
                                     quality );
    if (fromLowerBound != NIL(mdd_t)) {
      if (!bdd_equal(fromBetween,fromLowerBound))  {
        if (mdd_lequal(fromSubset, fromLowerBound, 1, 0)) {
          mdd_free(fromSubset);
          fromSubset = bdd_approx_remap_ua(fromLowerBound,
                                           (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                           nvars, frontierApproxThreshold,
                                           quality);
        }
      }
    }
    break;
  }

  return (fromSubset);
} /* end of ComputeSubsetBasedOnMethod */


static mdd_t *
AddStates(mdd_t *a, mdd_t *b, int freeA, int freeB)
{
  mdd_t *result;
  if ((a != NULL) && (b != NULL)) {
    result = mdd_or(a, b, 1, 1);
  } else if ((a == NULL) && (b == NULL)) {
    result = a;
  } else if (a == NULL) {
    result = mdd_dup(b);
  } else {
    result =  mdd_dup(a);
  }
  if ((freeA == FSM_HD_FREE) && (a != NULL))  mdd_free(a);
  if ((freeB == FSM_HD_FREE) && (b != NULL)) mdd_free(b);
  return result;
} /* end of AddStates */



static mdd_t *
SubtractStates(mdd_t *a, mdd_t *b, int freeA, int freeB)
{
  mdd_t *result;
  if ((a != NULL) && (b != NULL)) {
    result = mdd_and(a, b, 1, 0);
  } else if ((a == NULL) && (b == NULL)) {
    result = a;
  } else if (b == NULL) {
    result =  mdd_dup(a);
  } else { /* a = NULL , b != NULL */
    result = NULL;
  }
  if ((freeA == FSM_HD_FREE) && (a != NULL))  mdd_free(a);
  if ((freeB == FSM_HD_FREE) && (b != NULL)) mdd_free(b);
  return result;
} /* end of SubtractStates */


static int
RandomCompare(const void *ptrX,
              const void *ptrY)
{
  FsmHdSizeStruct_t *sizeStruct1 = (FsmHdSizeStruct_t *)ptrX;
  FsmHdSizeStruct_t *sizeStruct2 = (FsmHdSizeStruct_t *)ptrY;

  if (sizeStruct1->rnd > sizeStruct2->rnd) return 1;
  else if (sizeStruct1->rnd < sizeStruct2->rnd) return -1;
  else return 0;
} /* end of BddSizeCompare */