src/fsm/fsmReach.c

Go to the documentation of this file.

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

static char rcsid[] UNUSED = "$Id: fsmReach.c,v 1.102 2009/04/11 01:40:54 fabio Exp $";

/*---------------------------------------------------------------------------*/
/* Constant declarations                                                     */
/*---------------------------------------------------------------------------*/
/* HD constants */
#define FSM_HD_NONGREEDY 0 /* flag to indicate non-greedy computation at
                              dead-ends */
#define FSM_HD_GREEDY 1 /* flag to indicate greedy computation at dead-ends */
#define FSM_HD_LARGE_SIZE 100000 /* size considered large for reached */
#define FSM_HD_MID_SIZE 50000    /* size considered medium for reached*/
#define FSM_HD_MINT_GROWTH 0.5   /* factor measuring minterm jump of reached*/
#define FSM_HD_GROWTH_RATE 0.04 /* Growth rate for reached in measuring slow
                                 * growth. */
#define FSM_MDD_DONT_FREE 0  /* flag to indicate that bdd should not be freed. */
#define FSM_MDD_FREE 1       /* flag to indicate that bdd should be freed.*/
#define FSM_HD_NUM_SLOW_GROWTHS 5 /* number of iterations of allowed
                                   * slow growth of reached. */
/*---------------------------------------------------------------------------*/
/* Structure declarations                                                    */
/*---------------------------------------------------------------------------*/

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

static int CheckImageValidity(mdd_manager *mddManager, mdd_t *image, array_t *domainVarMddIdArray, array_t *quantifyVarMddIdArray);
static int ComputeNumberOfBinaryStateVariables(mdd_manager *mddManager, array_t *mddIdArray);
static mdd_t * AddStates(mdd_t *a, mdd_t *b, int freeA, int freeB);
static void RandomSimulation(int simNVec, Fsm_Fsm_t *fsm, Fsm_RchType_t approxFlag, mdd_t *initialStates, mdd_t **reachableStates, mdd_t **fromLowerBound, FsmHdStruct_t *hdInfo);
static void PrintStatsPerIteration(Fsm_RchType_t approxFlag, int nvars, int depth, FsmHdStruct_t *hdInfo, mdd_manager *mddManager, mdd_t *reachableStates, mdd_t *fromLowerBound, array_t *mintermVarArray);
static void HdInduceFullDeadEndIfNecessary(mdd_manager *mddManager, Img_ImageInfo_t *imageInfo, mdd_t **fromLowerBound, mdd_t **fromUpperBound, mdd_t **reachableStates, mdd_t **image, FsmHdStruct_t *hdInfo, int *depth, int shellFlag, array_t *onionRings, array_t *mintermVarArray, int oldSize, double oldMint, int verbosityLevel);
static mdd_t * ComputeInitialStates(Fsm_Fsm_t *fsm, int shellFlag, boolean printWarning);
static int CheckStatesContainedInInvariant(mdd_t *reachableStates, array_t *invarStates);
static void HdComputeReachabilityParameters(mdd_manager *mddManager, Img_ImageInfo_t *imageInfo, mdd_t **reachableStates, mdd_t **fromUpperBound, mdd_t **fromLowerBound, mdd_t **image, FsmHdStruct_t *hdInfo, mdd_t *initialStates, array_t *mintermVarArray, int *depth, int printStep, int shellFlag, array_t *onionRings, int verbosityLevel);
static int InitializeIncrementalParameters(Fsm_Fsm_t *fsm, Ntk_Network_t *network, array_t **arrayOfRoots, st_table **tableOfLeaves);
static mdd_t * IncrementalImageCompute(Ntk_Network_t *network, Fsm_Fsm_t *fsm, mdd_manager *mddManager, mdd_t *fromLowerBound, mdd_t *fromUpperBound, mdd_t *toCareSet, array_t *arrayOfRoots, st_table *tableOfLeaves, array_t *smoothVarArray, array_t *relationArray, int numLatch);
static void PrintOnionRings(Fsm_Fsm_t *fsm, int printStep, Fsm_RchType_t approxFlag, int nvars);
static array_t * GenerateGuidedSearchSequenceArray(array_t *guideArray);
static void ComputeReachabilityParameters(mdd_manager *mddManager, Img_ImageInfo_t *imageInfo, Fsm_RchType_t approxFlag, mdd_t **reachableStates, mdd_t **fromUpperBound, mdd_t **fromLowerBound, mdd_t **image, FsmHdStruct_t *hdInfo, mdd_t *initialStates, array_t *mintermVarArray, int *depth, int printStep, int shellFlag, array_t *onionRings, int verbosityLevel, boolean guidedSearchMode, mdd_t *hint, int *hintDepth, boolean *notConverged);

/*---------------------------------------------------------------------------*/
/* Definition of exported functions                                          */
/*---------------------------------------------------------------------------*/
mdd_t *
Fsm_FsmComputeReachableStates(
  Fsm_Fsm_t *fsm,
  int incrementalFlag,
  int verbosityLevel,
  int printStep,
  int depthValue,
  int shellFlag,
  int reorderFlag,
  int reorderThreshold,
  Fsm_RchType_t approxFlag,
  int ardc,
  int recompute,
  array_t *invarStates,
  boolean printWarning,
  array_t *guideArray)
{
  /* BFS variables */
  Img_ImageInfo_t *imageInfo = NIL(Img_ImageInfo_t); /* image info structure*/
  Img_ImageInfo_t *oldImageInfo = NIL(Img_ImageInfo_t);
  mdd_t           *reachableStates;/* reachable states */
  mdd_t           *unreachableStates;
  mdd_t           *fromLowerBound; /* new states in each iteration */
  mdd_t           *fromUpperBound; /* set to reachable states */
  mdd_t           *image, *newImage;/* image computed at each iteration */
  mdd_t           *initialStates;  /* initial states */
  mdd_t           *toCareSet;      /* the complement of the reached set */
  int              depth = 0;      /* depth upto which the computation is
                                    * performed.
                                    */
  graph_t         *partition, *oldPartition;
  mdd_manager     *mddManager = Fsm_FsmReadMddManager(fsm); /* mdd manager */
  Ntk_Network_t   *network = fsm->network; /* network */
  bdd_reorder_type_t currentMethod =
    Ntk_NetworkReadDynamicVarOrderingMethod(network); /* current reordering
                                                       * method */
  int firstTimeFlag = TRUE, firstReorderFlag = FALSE;
                         /* a flag to indicate the start of the computation */
  Fsm_RchStatus_t rchStatus = Fsm_Rch_Under_c;
                         /* under approximation of the reached set computed */
  int reachableStateSetSize; /* bdd size of the reached set */
  array_t *mintermVarArray; /* array of present state variables */
  array_t *onionRings = NIL(array_t);
                            /* onionringarray for shellFlag computation */

  /* Incremental flag  variables */
  array_t *arrayOfRoots = NIL(array_t);
  st_table *tableOfLeaves = NIL(st_table);
  int numLatch = 0;
  array_t *relationArray = NIL(array_t), *smoothVarArray = NIL(array_t);

  boolean notConverged = FALSE; /* flag that says that the fixpoint procedure
                                   did not converge. */
  /* HD variables */
  FsmHdStruct_t *hdInfo = NULL;
  int nvars;       /* number of present state variables */


  /* Simulation variables */
  int simNVec = 0;                 /* number of simulation vectors. */
  char *simString = Cmd_FlagReadByName("rch_simulate");

  /* Guided Search */
  array_t *hintDepthArray = NIL(array_t);
  int hintDepth = -1;
  mdd_t   *hint = NIL(mdd_t); /* iterates over guide array */
  int      hintnr;            /* idem                      */
  boolean guidedSearchMode = FALSE;
  boolean guideArrayAllocated = FALSE;

  boolean invariantFailed = FALSE;

  /* initializations */
  /* if fsm is a submachine, takes realPsVars to count minterm correctly. */
  mintermVarArray = Fsm_FsmReadPresentStateVars(fsm);

  /* compute number of present state variables */
  nvars = ComputeNumberOfBinaryStateVariables(mddManager, mintermVarArray);
  assert(nvars > 0);

  /* initializations */
  if (recompute) {
    if (incrementalFlag &&
        ((approxFlag == Fsm_Rch_Hd_c) || (approxFlag == Fsm_Rch_Oa_c))) {
      fprintf(vis_stdout,
              "** rch error: Incremental flag and HD computation are not compatible\n");
      return NIL(mdd_t);
    }
    FsmResetReachabilityFields(fsm, approxFlag);
  } else if ((!shellFlag) || Fsm_FsmTestReachabilityOnionRingsUpToDate(fsm)) {
    /* If already computed and exact , just return a copy. */
    reachableStates = Fsm_FsmReadReachableStates(fsm);
    if (reachableStates != NIL(mdd_t)){
      if (printWarning) {
        fprintf(vis_stdout, "** rch info: Reachable states have been previously computed.\n");
        fprintf(vis_stdout, "** rch info: Not recomputing reachable states.\n");
        fprintf(vis_stdout, "** rch info: Use compute_reach -F to recompute if necessary.\n");
      }
      if (!shellFlag) {
        if (printWarning) {
          if (FsmReadReachabilityComputationMode(fsm) == Fsm_Rch_Bfs_c) {
            fprintf(vis_stdout, "** rch info: All previous computations done using BFS (-A 0 option).\n");
          } else {
            fprintf(vis_stdout, "** rch info: Some previous computations done using BFS/DFS mode (-A 1 or -g option).\n");
          }
        }
        return (mdd_dup(reachableStates));
      } else if (Fsm_FsmReadReachabilityOnionRings(fsm)) {
        /* if onion rings exist, since they are up-to-date, return the
         * reachable states
         */
        if (printWarning) {
          if (Fsm_FsmReadReachabilityOnionRingsMode(fsm) == Fsm_Rch_Bfs_c) {
            fprintf(vis_stdout, "** rch info: Onion rings have been computed using BFS (-A 0 option)\n");
          } else {
            fprintf(vis_stdout, "** rch info: Some onion rings have been computed using BFS/DFS (-A 1 or -g option)\n");
          }
        }
        if (printStep)
          PrintOnionRings(fsm, printStep, approxFlag, nvars);
        return (mdd_dup(reachableStates));
      }
    }

    /* if some computed states exist and they violate an invariant,
       return the current set of reachable states */
    reachableStates = Fsm_FsmReadCurrentReachableStates(fsm);
    if (reachableStates != NIL(mdd_t)) {
      if (!shellFlag) {
        /* if an underApprox  exists, check invariant, if any */
        if (invarStates != NIL(array_t)) {
          if (!CheckStatesContainedInInvariant(reachableStates, invarStates)) {
            if (printWarning) {
              fprintf(vis_stdout, "** rch info: Invariant violation using previously computed states\n");
            }
            /* return violating reachable states */
            return (mdd_dup(reachableStates));
          }
        }
      } else if (Fsm_FsmReadReachabilityOnionRings(fsm)) {
        /* an invariant fails. */
        if (invarStates != NIL(array_t)) {
          if (!CheckStatesContainedInInvariant(reachableStates, invarStates)) {
            if (printWarning) {
              fprintf(vis_stdout, "** rch info: Invariant violation using previously computed states\n");
              if (Fsm_FsmReadReachabilityOnionRingsMode(fsm) == Fsm_Rch_Bfs_c) {
                fprintf(vis_stdout, "** rch info: Previous onion rings computed in BFS (-A 0 option).\n");
              } else {
                fprintf(vis_stdout, "** rch info: Some set of onion rings have been computed using BFS/DFS (-A 1 or -g option)\n");
              }
              fprintf(vis_stdout, "** rch info: Use compute_reach -F to recompute if necessary.\n");
            }
            if (printStep) {
              PrintOnionRings(fsm, printStep, approxFlag, nvars);
            }
            return (mdd_dup(reachableStates));
          }
        }
      }
    }
    /* if an over approx exists, check invariant, if any */
    if ((invarStates != NIL(array_t)) && (!shellFlag) && (!depthValue) &&
        (Fsm_FsmReadOverApproximateReachableStates(fsm)
         != NIL(array_t))) {
      if (FsmArdcCheckInvariant(fsm, invarStates)) {
        fprintf(vis_stdout, "** rch info: Over approximation exists, using over approximation for invariant checking.\n");
         return(Fsm_ArdcGetMddOfOverApproximateReachableStates(fsm));
      }
    }

    if (incrementalFlag &&
        ((approxFlag == Fsm_Rch_Hd_c) || (approxFlag == Fsm_Rch_Oa_c))) {
      fprintf(vis_stdout,
              "** rch error: Incremental flag and HD computation are not compatible\n");

      return NIL(mdd_t);
    }
  }
  /* onion rings make no sense with Over approx */
  if (shellFlag && (approxFlag == Fsm_Rch_Oa_c)) {
    fprintf(vis_stdout, "** rch error: Can't generate onion rings with over approx\n");
    return NIL(mdd_t);
  }

  /* depth value makes no sense with over approximation */
  if (depthValue && (approxFlag == Fsm_Rch_Oa_c)) {
    fprintf(vis_stdout, "** rch error: Can't generate over approx with depth Value\n");
    return NIL(mdd_t);
  }

  /* guided search cannot be done with Over approximation */
  if ((guideArray != NIL(array_t)) && (approxFlag == Fsm_Rch_Oa_c)) {
    fprintf(vis_stdout, "Guided search is not implemented with Over approximations\n");
    return NIL(mdd_t);
  }

  /* initializations */
  reachableStateSetSize = 0;
  if ((approxFlag == Fsm_Rch_Hd_c) || (guideArray != NIL(array_t))) {
    /* this structure is needed for HD and guidedSearch */
    hdInfo = FsmHdStructAlloc(nvars);
  }

  simNVec = 0;
  if (simString != NIL(char)) {
    simNVec = strtol(simString, NULL, 10);
    if (simNVec <= 0)  simNVec = 0;
  }
  /* initialize onion ring array */
  if (shellFlag) {
    if (Fsm_FsmReadReachabilityOnionRings(fsm) == NIL(array_t)) {
      onionRings = array_alloc(mdd_t *, 0);
    } else {
      onionRings = Fsm_FsmReadReachabilityOnionRings(fsm);
    }
  }

  /* Computes ARDCs, and  overapproximation by SSD traversal, if required. */
  if (approxFlag == Fsm_Rch_Oa_c || ardc) {
    long initialTime = 0;
    long finalTime;

    assert(!incrementalFlag);

    if (verbosityLevel > 0)
      initialTime = util_cpu_time();

    (void)Fsm_FsmComputeOverApproximateReachableStates(fsm, incrementalFlag,
        verbosityLevel, printStep, depthValue, shellFlag, reorderFlag,
        reorderThreshold, recompute);

    if (approxFlag == Fsm_Rch_Oa_c) {
      if (hdInfo != NIL(FsmHdStruct_t)) FsmHdStructFree(hdInfo);
      return(Fsm_ArdcGetMddOfOverApproximateReachableStates(fsm));
    }

    if (verbosityLevel > 0) {
      finalTime = util_cpu_time();
      Fsm_ArdcPrintReachabilityResults(fsm, finalTime-initialTime);
    }
  }

  if(Cmd_FlagReadByName("linear_compute_range")) {
    partition = Part_NetworkCreatePartition(network, 0, "temp", (lsList)0,
        (lsList)0, NIL(mdd_t), Part_Fine_c, 0, 0, 0, 0);
    oldPartition = fsm->partition;
    oldImageInfo = fsm->imageInfo;
    fsm->imageInfo = 0;
    imageInfo = Fsm_FsmReadOrCreateImageInfoForComputingRange(fsm, 0, 1);
    unreachableStates = Img_ImageGetUnreachableStates(imageInfo);
    fprintf(vis_stdout, "%-20s%10g\n", "reachable states =",
                     mdd_count_onset(mddManager, unreachableStates,
                                     fsm->fsmData.presentStateVars));
    fflush(vis_stdout);
    /* get unreachable states */
    Img_ImageInfoFree(imageInfo);
    fsm->imageInfo = oldImageInfo;
    fsm->partition = oldPartition;
    Part_PartitionFree(partition);
    exit(0);
  }
  /* Compute the set of initial states. Start with the
   * underapproximation of Reached states if it exists or sum of
   * onion rings if shell flag is specified or the initial states of
   * the fsm.
   */
  initialStates = ComputeInitialStates(fsm, shellFlag, printWarning);
  if (initialStates == NIL(mdd_t)) {
    fprintf(vis_stdout, "** rch error: No initial states computed.");
    fprintf(vis_stdout, " No reachability will be performed.\n");
    if (hdInfo != NIL(FsmHdStruct_t)) FsmHdStructFree(hdInfo);
    return NIL(mdd_t);
  }

  if (incrementalFlag){
    /* Note: if incrementalflag is set, no imageinfo is created! */
    numLatch = InitializeIncrementalParameters(fsm, network, &arrayOfRoots,
                                               &tableOfLeaves);
    smoothVarArray = array_join(fsm->fsmData.presentStateVars,
                                fsm->fsmData.inputVars);
    relationArray = array_alloc(mdd_t *, numLatch+1);
  } else {
    /* Create an imageInfo for image computations, if not already created. */
    imageInfo = Fsm_FsmReadOrCreateImageInfo(fsm, 0, 1);
    if (ardc) {
      Fsm_ArdcMinimizeTransitionRelation(fsm, Img_Forward_c);
    }
  }

  /* Get approximate traversal options */
  if ((!incrementalFlag) && (approxFlag == Fsm_Rch_Hd_c)) {
    /* only partial image allowed. No approximation of frontier set */
    hdInfo->onlyPartialImage =  Fsm_FsmHdOptionReadOnlyPartialImage(hdInfo->options);
    /* states which can produce no new successors. */
    hdInfo->interiorStates = bdd_zero(mddManager);
  }

  /* pertaining to simulation */
  if (simNVec > 0) {
    RandomSimulation(simNVec, fsm, approxFlag, initialStates,
                     &reachableStates, &fromLowerBound, hdInfo);
  } else {
    /* start reached set with initial states */
    reachableStates = mdd_dup(initialStates);
    fromLowerBound = mdd_dup(reachableStates);
  }

  /* initialize variables to print */
  if ((!incrementalFlag) && (approxFlag == Fsm_Rch_Hd_c)) {

    FsmHdStatsComputeSizeAndMinterms(mddManager, reachableStates,
                           mintermVarArray, nvars,
                           Fsm_Hd_Reached, hdInfo->stats);
    FsmHdStatsComputeSizeAndMinterms(mddManager, fromLowerBound,
                           mintermVarArray, nvars,
                           Fsm_Hd_From, hdInfo->stats);
    FsmHdStatsSetMintermFromSubset(hdInfo->stats, FsmHdStatsReadMintermFrom(hdInfo->stats));
    FsmHdStatsSetSizeFromSubset(hdInfo->stats, FsmHdStatsReadSizeFrom(hdInfo->stats));
  }

  if (printStep && ((depth % printStep) == 0) && (approxFlag == Fsm_Rch_Hd_c)){
    fprintf(vis_stdout, "\nIndex: R = Reached; I = Image, F = Frontier, FA = Approximation of Frontier\n");
    fprintf(vis_stdout, "       f = Minterms in f; |f| = Bdd nodes in f\n\n");
  }
  /* initialize for the first iteration */
   fromUpperBound = reachableStates;

  /* hint array */
  if (guideArray == NIL(array_t)) {
    guideArray = array_alloc(mdd_t *, 0);
    guideArrayAllocated = TRUE;
  }
  /* pad with 1 when no hints are provided */
  if (array_n(guideArray) == 0) {
    array_insert_last(mdd_t *, guideArray, bdd_one(mddManager));
  }
  /* depth sequence for hints. 1-to-1 correspondence with guidearray */
  /* -1 implies apply hint to convergence. */
  hintDepthArray = GenerateGuidedSearchSequenceArray(guideArray);

  assert(array_n(hintDepthArray) == array_n(guideArray));

  if (array_n(guideArray) == 0) hint = mdd_one(mddManager);

  arrayForEachItem(mdd_t *, guideArray, hintnr, hint){

    /* check the depth for which this hint is to be applied. */
    hintDepth = array_fetch(int, hintDepthArray, hintnr);

    if(hintnr == 0 && mdd_is_tautology(hint, 1) && array_n(guideArray) == 1){
      assert(hintDepth == -1);
    }  else {
      assert(!incrementalFlag);
      guidedSearchMode  = TRUE;
      Fsm_InstantiateHint(fsm, hint, 1, 0, Ctlp_Underapprox_c,
                          verbosityLevel > 1);
      if (hintnr < (array_n(guideArray)-1))
        fprintf(vis_stdout, "**GS info: Instantiating  hint number %d\n", hintnr+1);
      else
        fprintf(vis_stdout, "**GS info: Restoring original transition relation\n");

      if (approxFlag == Fsm_Rch_Hd_c) {
        /* clean up all the invalid state information */
        if (hdInfo->interiorStates != NIL(mdd_t))
          mdd_free(hdInfo->interiorStates);
        hdInfo->interiorStates = NIL(mdd_t); /* for every new hint, this is invalid. */
        if (hdInfo->interiorStateCandidate != NIL(mdd_t))
          mdd_free(hdInfo->interiorStateCandidate);
        hdInfo->interiorStateCandidate = NIL(mdd_t); /* for every new hint, this is invalid */
        hdInfo->imageOfReachedComputed = FALSE; /* since it doesn't matter for every new hint */
        if (hdInfo->deadEndResidue != NIL(mdd_t))
          mdd_free(hdInfo->deadEndResidue);
        hdInfo->deadEndResidue = NIL(mdd_t);
      }
      /* generate a frontier by doing a dead-end computation if the
         frontier is empty. This is non-greedy for BFS and greedy for
         HD. */
      if (mdd_is_tautology(fromLowerBound, 0)) {
        ComputeReachabilityParameters(mddManager, imageInfo, approxFlag,
                                      &reachableStates, &fromUpperBound,
                                      &fromLowerBound, &image, hdInfo,
                                      initialStates,
                                      mintermVarArray,
                                      &depth, printStep,
                                      shellFlag, onionRings,
                                      verbosityLevel,
                                      guidedSearchMode, hint, &hintDepth,
                                      &notConverged);
      }
    }
    /* initialize */
    notConverged = FALSE;
    if (!incrementalFlag)
      Img_ImageInfoResetUseOptimizedRelationFlag(imageInfo);
    /* Main loop of reachability computation. */
    /* Iterate until fromLowerBound is empty or all states are reached. */
    while (!mdd_is_tautology(fromLowerBound, 0)) {
      if(depth > 0 && !incrementalFlag)
        Img_ImageInfoSetUseOptimizedRelationFlag(imageInfo);
      /* fromLowerBound is the "onion shell" of states just reached. */
      if ((shellFlag && (approxFlag != Fsm_Rch_Hd_c) &&
          (!firstTimeFlag ||
           (Fsm_FsmReadReachabilityOnionRings(fsm) == NIL(array_t)))) ||
          (shellFlag && (approxFlag == Fsm_Rch_Hd_c) &&
          (Fsm_FsmReadReachabilityOnionRings(fsm) == NIL(array_t)) &&
          firstTimeFlag)) {
        array_insert_last(mdd_t*, onionRings, mdd_dup(fromLowerBound));
      }

      /* if it is the print step, print out the iteration and the */
      if (printStep && ((depth % printStep) == 0)) {
        int step;
        /* for BFS, preserve earlier format of printing */
        step = (shellFlag) ? array_n(onionRings) : Fsm_FsmGetReachableDepth(fsm)+depth;
        PrintStatsPerIteration(approxFlag,
                               nvars,
                               step,
                               hdInfo,
                               mddManager,
                               reachableStates,
                               fromLowerBound,
                               mintermVarArray);
      }

      /* Check if invariant contains the new states.  In case of HD
       * computation, check the containment of the entire reached set.
       */
      if (invarStates != NIL(array_t)) {
        mdd_t *temp;
        temp = (approxFlag == Fsm_Rch_Hd_c) ? reachableStates : fromLowerBound;
        if (!CheckStatesContainedInInvariant(temp, invarStates)) {
          notConverged = TRUE;
          invariantFailed = TRUE; /* flag used for early termination */
          break;
        }
      }

      if (mdd_is_tautology(reachableStates, 1)) {
        break;
      }

      if ((depthValue && (depth == depthValue)) ||
          /* limited depth hint is implemented this way */
          (hintDepth == 0)) {
        /* No more steps */
        notConverged = TRUE;
        break;
      }

      if (reorderFlag && !firstReorderFlag &&
          (reachableStateSetSize >= reorderThreshold )){
        firstReorderFlag = TRUE;
        (void) fprintf(vis_stdout, "Dynamic variable ordering forced ");
        (void) fprintf(vis_stdout, "with method sift\n");
        Ntk_NetworkSetDynamicVarOrderingMethod(network,
                                               BDD_REORDER_SIFT,
                                               BDD_REORDER_VERBOSITY_DEFAULT);
        bdd_reorder(Ntk_NetworkReadMddManager(network));
        Ntk_NetworkSetDynamicVarOrderingMethod(network, currentMethod,
                                               BDD_REORDER_VERBOSITY_DEFAULT);
      }

      firstTimeFlag = FALSE;
      /* careSet for image computation is the set of all states not reached
       * so far.
       */
      toCareSet = mdd_not(reachableStates);

      /*
       * Image of some set between lower and upper, where we care
       * about the image only on unreached states.
       */
      if (incrementalFlag){
        image = IncrementalImageCompute(network, fsm, mddManager,
                                        fromLowerBound,
                                        fromUpperBound, toCareSet,
                                        arrayOfRoots, tableOfLeaves,
                                        smoothVarArray, relationArray,
                                        numLatch);
      } else{
        assert(!incrementalFlag);

        if (approxFlag == Fsm_Rch_Hd_c) {
          /* allow partial image computation */
          Img_ImageAllowPartialImage(imageInfo, TRUE);
        }

        /* compute the image of this iteration */
        image = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                      fromLowerBound,
                                                      fromUpperBound,
                                                      toCareSet);



        /* record if partial image or not */
        if (approxFlag == Fsm_Rch_Hd_c) {
          /* Check if partial image computation, disallow partial image
           * computation.
           */
          hdInfo->partialImage = Img_ImageWasPartial(imageInfo);
          /* record this to prevent dead-end computation. */
          hdInfo->imageOfReachedComputed = (!hdInfo->partialImage) &&
            mdd_equal(reachableStates, fromLowerBound);
          /* potential candidates for interior states (only if not a
           * partial image)
           */
          if (!hdInfo->partialImage) hdInfo->interiorStateCandidate = mdd_dup(fromLowerBound);
        }
      }

      /* free the used bdds */
      mdd_free(toCareSet);

      /* New lower bound is the states just reached. */
      mdd_free(fromLowerBound);
      fromLowerBound = mdd_and(image, reachableStates, 1, 0);
      depth++;
      hintDepth--;


      if ((!incrementalFlag) && (approxFlag == Fsm_Rch_Hd_c)) {
        FsmHdStatsComputeSizeAndMinterms(mddManager, image,
                                         mintermVarArray, nvars,
                                         Fsm_Hd_To, hdInfo->stats);
      }

      /* need image only if HD and not dead-end */
      if (approxFlag != Fsm_Rch_Hd_c) {
          mdd_free(image);
          image = NULL;
      } else {
        if (hdInfo->imageOfReachedComputed ||   hdInfo->onlyPartialImage ||
          (hdInfo->slowGrowth >= FSM_HD_NUM_SLOW_GROWTHS) ||
          mdd_is_tautology(fromLowerBound, 0)) {
          mdd_free(image);
          image = NULL;
        }
      }

      /* computes reachable states, fromlowerBound and fromupperbound
       * for BFS.  For HD it computes the above as well as interior
       * statess, dead-end residue, etc. All the other parameteres are
       * updated. For guided search, this call does not perform a
       * non-greedy dead-end computation, even if the frontier is 0. */
       ComputeReachabilityParameters(mddManager, imageInfo, approxFlag,
                                    &reachableStates,
                                    &fromUpperBound,
                                    &fromLowerBound,
                                    &image,
                                     hdInfo,
                                    initialStates,
                                    mintermVarArray,
                                     &depth,  printStep,
                                     shellFlag, onionRings,
                                    verbosityLevel, FALSE, hint,
                                    &hintDepth,
                                    &notConverged);


       if (mdd_is_tautology(fromLowerBound, 0) &&
           (approxFlag == Fsm_Rch_Bfs_c) &&
           (!guidedSearchMode) && (!depthValue))
         depth--;

    } /* end of main while loop (while frontier != 0) */

    /* stop iterating over hint if invariants are violated or all
       states are reachable or limited depth value has been
       reached. */
    if (invariantFailed) break;
    if (mdd_is_tautology(reachableStates, 1)) {
      break;
    }
    if (depthValue && (depth == depthValue)) {
      /* No more steps */
      break;
    }

  }/* for each hint */


  if(!incrementalFlag &&
     Img_ImageInfoObtainMethodType(imageInfo) == Img_Linear_c &&
     Img_ImageInfoObtainOptimizeType(imageInfo) == Opt_NS &&
     Img_IsTransitionRelationOptimized(imageInfo)) {
    Img_ImageInfoResetUseOptimizedRelationFlag(imageInfo);
    fprintf(vis_stdout, "%-20s%10g\n", "reachable states =",
                 mdd_count_onset(mddManager, reachableStates,
                                 fsm->fsmData.presentStateVars));
    fromLowerBound = mdd_dup(reachableStates);
    fromUpperBound = reachableStates;
    while (1) {
      toCareSet = mdd_dup(reachableStates);
      image = Img_ImageInfoComputeFwdWithDomainVars(imageInfo,
                                                  fromLowerBound,
                                                  fromUpperBound,
                                                  toCareSet);
      newImage = bdd_or(image, initialStates, 1, 1);
      bdd_free(image);
      image = newImage;
      if(mdd_equal(fromLowerBound, image)) {
        mdd_free(toCareSet);
        mdd_free(fromLowerBound);
        mdd_free(reachableStates);
        fromLowerBound = image;
        reachableStates = mdd_dup(fromLowerBound);
        fromUpperBound = reachableStates;
        fprintf(vis_stdout, "%-20s%10g\n", "reachable states =",
                 mdd_count_onset(mddManager, reachableStates,
                                 fsm->fsmData.presentStateVars));
        break;
      }
      mdd_free(toCareSet);
      mdd_free(fromLowerBound);
      mdd_free(reachableStates);
      fromLowerBound = image;
      reachableStates = mdd_dup(fromLowerBound);
      fromUpperBound = reachableStates;

      fprintf(vis_stdout, "%-20s%10g\n", "reachable states =",
                 mdd_count_onset(mddManager, reachableStates,
                                 fsm->fsmData.presentStateVars));
    }
  }


  /* clean up for Guided search mainly, only if needed */
  /* we rely on the fact that the variable hint is set to the last hint from
   * the previous loop.
   */
  if (guidedSearchMode && !mdd_is_tautology(hint, 1))
    Fsm_CleanUpHints(fsm, 1, 0, printStep);
  array_free(hintDepthArray);

  /* If (1) all states are reachable, or (2) not limited depth computation
     and not a converge with a hint, or (3) limited depth and not
     guided search mode and BFS and frontier = 0, or (4) in guided
     search, if limited depth and frontier = 0 then one full dead-end
     should have occured */
  if(mdd_is_tautology(reachableStates, 1) ||
     ((notConverged == FALSE) && mdd_is_tautology(hint, 1)) ||
     (!guidedSearchMode && (notConverged == TRUE) &&
      mdd_is_tautology(fromLowerBound, 0) && (approxFlag == Fsm_Rch_Bfs_c)) ||
     (guidedSearchMode && hdInfo->deadEndWithOriginalTR &&
      mdd_is_tautology(fromLowerBound, 0) && (approxFlag == Fsm_Rch_Bfs_c))) {

    rchStatus = Fsm_Rch_Exact_c;
    if (mdd_is_tautology(reachableStates, 1) &&
        !mdd_is_tautology(fromLowerBound, 0) &&
        printStep && ((depth % printStep) == 0) ) {
      int step;
      /* for BFS, preserve earlier format of printing */
      step = (shellFlag) ?
        array_n(onionRings) :
        Fsm_FsmGetReachableDepth(fsm)+depth;
      PrintStatsPerIteration(approxFlag, nvars, step, hdInfo,
                             mddManager, reachableStates,
                             fromLowerBound, mintermVarArray);
    }
  }

  /* clean up the hints if allocated here */
  if (guideArrayAllocated) mdd_array_free(guideArray);

  if (hdInfo != NIL(FsmHdStruct_t)) {
    FsmHdStructFree(hdInfo);
  }

  mdd_free(fromLowerBound);
  mdd_free(initialStates);

  /* Update FSM reahcability fields */

  /* indicate whether the shells are consistent with the reachable states */
  if (!shellFlag) {
    /* if any states were added from a previous computation, shells
       are not up-to-date */
    if (Fsm_FsmReadCurrentReachableStates(fsm) != NIL(mdd_t)) {
      if (!mdd_lequal(reachableStates, Fsm_FsmReadCurrentReachableStates(fsm), 1, 1)) {
        FsmSetReachabilityOnionRingsUpToDateFlag(fsm, FALSE);
      } /* otherwise the previous iteration was computed with onion
      rings, in which case it is still up-to-date or the previous
      iteration was computed without onion rings and it was set to not
      up-to-date in the previous iteration */
    } else {
      assert(Fsm_FsmReadReachabilityOnionRings(fsm) == NIL(array_t));
      /* not up-to-date because no onion rings */
      FsmSetReachabilityOnionRingsUpToDateFlag(fsm, FALSE);
    }

  } else {
    /* shells are up to date only if current set adds upto previously
       computed states or more */
    if (Fsm_FsmReadCurrentReachableStates(fsm)) {
      if (!mdd_lequal(Fsm_FsmReadCurrentReachableStates(fsm),
                    reachableStates, 1, 1)) {
        FsmSetReachabilityOnionRingsUpToDateFlag(fsm, FALSE);
      } else {
        FsmSetReachabilityOnionRingsUpToDateFlag(fsm, TRUE);
      }

    } else FsmSetReachabilityOnionRingsUpToDateFlag(fsm, TRUE);

    /* Onion rings are not BFS, i.e., HD or Guided search. */
    if ((approxFlag == Fsm_Rch_Hd_c) ||
        (Fsm_FsmReadReachabilityOnionRingsMode(fsm) == Fsm_Rch_Hd_c) ||
        guidedSearchMode)
      FsmSetReachabilityOnionRingsMode(fsm, Fsm_Rch_Hd_c);

    if (array_n(onionRings) == 0) {
      array_free(onionRings);
      onionRings = NIL(array_t);
    }
    FsmSetReachabilityOnionRings(fsm, onionRings);
  }

  /* update fsm structure */
  /* record depth of traversal, depth is consistent with reachable
     states, hence may not be consistent with the onion rings */
  /* replace reachable states if only more states have been computed now */
  if (Fsm_FsmReadCurrentReachableStates(fsm) != NIL(mdd_t)) {
    if (!mdd_lequal(reachableStates, Fsm_FsmReadCurrentReachableStates(fsm), 1, 1)) {
      Fsm_FsmFreeReachableStates(fsm);
      fsm->reachabilityInfo.reachableStates = mdd_dup(reachableStates);
      if (!shellFlag) fsm->reachabilityInfo.depth = Fsm_FsmGetReachableDepth(fsm) + depth;
      else if (Fsm_FsmReadReachabilityOnionRings(fsm))
        fsm->reachabilityInfo.depth = array_n(Fsm_FsmReadReachabilityOnionRings(fsm));
      /* set computation mode as HD if HD or Guided search, basically to indicate not BFS */
      if ((approxFlag == Fsm_Rch_Hd_c) ||
          (FsmReadReachabilityComputationMode(fsm)  == Fsm_Rch_Hd_c) ||
          guidedSearchMode) {
        FsmSetReachabilityComputationMode(fsm, Fsm_Rch_Hd_c);
      }
    }
  } else {
    fsm->reachabilityInfo.reachableStates = mdd_dup(reachableStates);
    if (!shellFlag) fsm->reachabilityInfo.depth = depth;
    else if (Fsm_FsmReadReachabilityOnionRings(fsm))
      fsm->reachabilityInfo.depth = array_n(Fsm_FsmReadReachabilityOnionRings(fsm));
    /* set computation mode as HD if HD or Guided search, basically to indicate not BFS */
      if ((approxFlag == Fsm_Rch_Hd_c) || guidedSearchMode) {
        FsmSetReachabilityComputationMode(fsm, Fsm_Rch_Hd_c);
      } else {
        FsmSetReachabilityComputationMode(fsm, Fsm_Rch_Bfs_c);
      }

  }
  /* indicate whether reachability analysis is completed or not */
  FsmSetReachabilityApproxComputationStatus(fsm, rchStatus);

  if (incrementalFlag){
    array_free(relationArray);
    array_free(arrayOfRoots);
    st_free_table(tableOfLeaves);
    array_free(smoothVarArray);
    /* Need to put in next state function data for the partition */
  }

  return reachableStates;
} /* end of Fsm_FsmComputeReachableStates */



mdd_t *
Fsm_MddMultiwayAndSmooth(mdd_manager *mddManager, array_t *mddArray,
                         array_t *smoothingVars)
{
  int i;
  mdd_t *resultMdd;
  mdd_t *product = mdd_one(mddManager);

  for (i = 0; i<array_n(mddArray); i++){
    mdd_t *mdd         = array_fetch(mdd_t*, mddArray, i);
    mdd_t *tempProduct = mdd_and(product, mdd, 1, 1);
    mdd_free(product);
    product = tempProduct;
  }

  if (array_n(smoothingVars) == 0) {
    resultMdd = mdd_dup(product);
  }
  else{
    resultMdd =  mdd_smooth(mddManager, product, smoothingVars);
  }

  mdd_free(product);
  return resultMdd;
}

void
Fsm_FsmReachabilityPrintResults(
  Fsm_Fsm_t *fsm,
  long elapseTime,
  int approxFlag)
{
  mdd_manager   *mddManager = Fsm_FsmReadMddManager(fsm);
  mdd_t         *reachableStates = fsm->reachabilityInfo.reachableStates;
  array_t       *psVarsArray;
  int           nvars;

  if (!reachableStates)
    return;

  psVarsArray = Fsm_FsmReadPresentStateVars(fsm);
  /* compute number of present state variables */
  nvars = ComputeNumberOfBinaryStateVariables(mddManager, psVarsArray);

  (void) fprintf(vis_stdout,"********************************\n");
  (void) fprintf(vis_stdout,"Reachability analysis results:\n");
  if ((FsmReadReachabilityComputationMode(fsm) == Fsm_Rch_Bfs_c) &&
      (Fsm_FsmReadReachableStates(fsm) != NIL(mdd_t))) {
    (void) fprintf(vis_stdout, "%-20s%10d\n", "FSM depth =",
                   fsm->reachabilityInfo.depth);
  } else {
    (void) fprintf(vis_stdout, "%-20s%10d\n", "computation depth =",
                   fsm->reachabilityInfo.depth);
  }
  if (nvars <= 1023) {
    (void) fprintf(vis_stdout, "%-20s%10g\n", "reachable states =",
                   mdd_count_onset(mddManager, reachableStates,
                                   fsm->fsmData.presentStateVars));
  } else {
    /*
    (void) fprintf(vis_stdout, "%-20s", "reachable states = ");
    bdd_print_apa_minterm(vis_stdout, reachableStates, nvars, (long)6);
    */
    EpDouble    epd;
    char        strValue[20];

    mdd_epd_count_onset(mddManager, reachableStates,
                        fsm->fsmData.presentStateVars, &epd);
    EpdGetString(&epd, strValue);
    (void) fprintf(vis_stdout, "%-20s%10s\n", "reachable states =", strValue);
  }
  (void) fprintf(vis_stdout, "%-20s%10d\n", "BDD size =",
                 mdd_size(reachableStates));
  (void) fprintf(vis_stdout, "%-20s%10g\n", "analysis time =",
                 (double)elapseTime/1000.0);
  if (Fsm_FsmReadReachableStates(fsm) != NIL(mdd_t)) {
    (void) fprintf(vis_stdout, "%-20s%10s\n", "reachability analysis = ", "complete");
  } else {
    (void) fprintf(vis_stdout, "%-20s%10s\n", "reachability analysis = ", "incomplete");
  }
}


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

static int
CheckImageValidity(mdd_manager *mddManager, mdd_t *image, array_t
                   *domainVarMddIdArray, array_t *quantifyVarMddIdArray)
{
  int i;
  array_t *imageSupportArray = mdd_get_support(mddManager, image);
  st_table *imageSupportTable = st_init_table(st_numcmp, st_numhash);
  for (i=0; i<array_n(imageSupportArray); i++){
    int mddId = array_fetch(int, imageSupportArray, i);
    (void) st_insert(imageSupportTable, (char *) (long) mddId, NIL(char));
  }
  for (i=0; i<array_n(domainVarMddIdArray); i++){
    int domainVarId;
    domainVarId = array_fetch(int, domainVarMddIdArray, i);
    assert(st_is_member(imageSupportTable, (char *)(long)domainVarId) == 0);
  }
  for (i=0; i<array_n(quantifyVarMddIdArray); i++){
    int quantifyVarId;
    quantifyVarId = array_fetch(int, quantifyVarMddIdArray, i);
    assert(st_is_member(imageSupportTable, (char *)(long)quantifyVarId) == 0);
  }
  st_free_table(imageSupportTable);
  array_free(imageSupportArray);
  return 1;
}

static int
ComputeNumberOfBinaryStateVariables(mdd_manager *mddManager,
                               array_t *mddIdArray)
{
  int mddId, i;
  mvar_type mddVar;
  int numEncodingBits, count = 0;
  array_t *mvar_list = mdd_ret_mvar_list(mddManager);

  arrayForEachItem(int, mddIdArray, i, mddId) {
    mddVar = array_fetch(mvar_type, mvar_list, mddId);
    numEncodingBits = mddVar.encode_length;
    count += numEncodingBits;
  }
  return count;
} /* end of ComputeNumberOfBinaryStateVariables */


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_MDD_FREE) && (a != NULL)) mdd_free(a);
  if ((freeB == FSM_MDD_FREE) && (b != NULL)) mdd_free(b);
  return result;
} /* end of AddStates */


static void
RandomSimulation(int simNVec,
                 Fsm_Fsm_t *fsm,
                 Fsm_RchType_t approxFlag,
                 mdd_t *initialStates,
                 mdd_t **reachableStates,
                 mdd_t **fromLowerBound,
                 FsmHdStruct_t *hdInfo)
{
  Ntk_Network_t   *network = fsm->network;

  /* Compute specified number of simulated states. */
  *reachableStates = Sim_RandomSimulate(network, simNVec, 0);
  *reachableStates = AddStates(*reachableStates, initialStates, FSM_MDD_FREE, FSM_MDD_DONT_FREE);
  *fromLowerBound = mdd_dup(*reachableStates);

  if (approxFlag == Fsm_Rch_Hd_c) {
    if (mdd_size(*fromLowerBound) >
        Fsm_FsmHdOptionReadFrontierApproxThreshold(hdInfo->options)) {
      /* Reset frontier if too large. */
      mdd_free(*fromLowerBound);
      *fromLowerBound = bdd_between(initialStates, *reachableStates);
      if (!mdd_equal(*fromLowerBound, *reachableStates)) {
        hdInfo->firstSubset = FALSE;
        hdInfo->deadEndResidue = mdd_and(*reachableStates,
                                         *fromLowerBound, 1, 0);
      }
    }
  }
  return;
} /* end of RandomSimulation */


static void
PrintStatsPerIteration(Fsm_RchType_t approxFlag,
                       int nvars,
                       int depth,
                       FsmHdStruct_t *hdInfo,
                       mdd_manager *mddManager,
                       mdd_t *reachableStates,
                       mdd_t *fromLowerBound,
                       array_t *mintermVarArray)
{
  if (approxFlag != Fsm_Rch_Hd_c) {
    int reachableStateSetSize = mdd_size(reachableStates);
    if (nvars <= 1023) {
      (void)fprintf(vis_stdout, "BFS step = %d\t|states| = %8g\tBdd size = %8ld\n",
                    depth, mdd_count_onset(mddManager, reachableStates,
                                           mintermVarArray),
                    (long)reachableStateSetSize);
      /*                            bdd_print_minterm(reachableStates); */
    } else {
      (void)fprintf(vis_stdout, "BFS step = %d\tBdd size = %8d\t|states| = ",
                    depth, reachableStateSetSize);
      (void) bdd_print_apa_minterm(vis_stdout, reachableStates, nvars, (long)6);
    }

  } else { /* HD */

    if (nvars <= 1023) {
      (void)fprintf(vis_stdout, "Step = %d, R = %8g, |R| = %8d,  I = %8g,  |I| = %8d\n",
                    depth,  FsmHdStatsReadMintermReached(hdInfo->stats),
                    FsmHdStatsReadSizeReached(hdInfo->stats), FsmHdStatsReadMintermTo(hdInfo->stats) ,
                    FsmHdStatsReadSizeTo(hdInfo->stats));
      (void)fprintf(vis_stdout, "Completed iteration %d at %g seconds\n",
                    depth, (util_cpu_time()/1000.0));

      (void)fprintf(vis_stdout, "Step = %d, F = %8g, |F| = %8d, FA = %8g, |FA| = %8d\n",
                    depth+1, FsmHdStatsReadMintermFrom(hdInfo->stats),
                    FsmHdStatsReadSizeFrom(hdInfo->stats),
                    FsmHdStatsReadMintermFromSubset(hdInfo->stats),
                    FsmHdStatsReadSizeFromSubset(hdInfo->stats));
    } else {
      int status;
      (void)fprintf(vis_stdout, "Step = %d, |R| = %8d, R = ",
                    depth, FsmHdStatsReadSizeReached(hdInfo->stats));
      status = bdd_print_apa_minterm(vis_stdout, reachableStates, nvars, (long)6);
      if (!status) {
        error_append("** rch error: Error in printing arbitrary precision minterm count of Reached");
      }

      (void)fprintf(vis_stdout, "Completed iteration %d at %g seconds\n",
                    depth, (util_cpu_time()/1000.0));
      (void)fprintf(vis_stdout, "Step = %d, |F| = %8d, F = %8g, |FA| = %8d, FA = ",
                    depth, FsmHdStatsReadSizeFrom(hdInfo->stats), 0.0,
                    FsmHdStatsReadSizeFromSubset(hdInfo->stats) );
      status = bdd_print_apa_minterm(vis_stdout, fromLowerBound, nvars, 6);
      if (!status) {
        error_append("** rch error: Error in printing arbitrary precision minterm count of From");
      }
    }
  }
  return;
} /* end of PrintStatsPerIteration */

static void
HdInduceFullDeadEndIfNecessary(mdd_manager *mddManager,
                               Img_ImageInfo_t *imageInfo,
                               mdd_t **fromLowerBound,
                               mdd_t **fromUpperBound,
                               mdd_t **reachableStates,
                               mdd_t **image,
                               FsmHdStruct_t *hdInfo,
                               int *depth,
                               int shellFlag,
                               array_t *onionRings,
                               array_t *mintermVarArray,
                               int oldSize,
                               double oldMint,
                               int verbosityLevel)
{
  int nvars =  Fsm_FsmHdOptionReadNumVars(hdInfo->options);
  /*
   * Compute a full (non-greedy) dead-end when a dead-end computation
   * hasn't occurred, and the jump in minterms is lower than MINT_GROWTH and
   * either the old reached size is between MID_SIZE and LARGE_SIZE and the
   * jump is MID_SIZE or the new size is larger than LARGE_SIZE and the jump
   * is LARGE_SIZE.
   */
  if ((hdInfo->deadEndComputation == FALSE) &&
      ((FsmHdStatsReadMintermReached(hdInfo->stats) - oldMint)/oldMint < FSM_HD_MINT_GROWTH) &&
      (((FsmHdStatsReadSizeReached(hdInfo->stats)  > FSM_HD_MID_SIZE) &&
        (FsmHdStatsReadSizeReached(hdInfo->stats) - oldSize > FSM_HD_MID_SIZE) &&
        (oldSize > FSM_HD_MID_SIZE) && (oldSize < FSM_HD_LARGE_SIZE)) ||
       ((FsmHdStatsReadSizeReached(hdInfo->stats) > FSM_HD_LARGE_SIZE) &&
        (FsmHdStatsReadSizeReached(hdInfo->stats) - oldSize > FSM_HD_LARGE_SIZE)))) {
    if (verbosityLevel >= 2) {
      fprintf(vis_stdout, "Full Dead End: New Reached = %d, Minterms = %g\n",
              FsmHdStatsReadSizeReached(hdInfo->stats), FsmHdStatsReadMintermReached(hdInfo->stats));
    }

    if (*image != NULL) {
      mdd_free(*image);
      *image = NULL;
    }
    /* throw away the reached with the jump in size */
    mdd_free(*reachableStates);
    if (shellFlag && onionRings) {
      mdd_free(array_fetch(mdd_t *, onionRings, array_n(onionRings)-1));
      array_insert(mdd_t *, onionRings, array_n(onionRings)-1, NIL(mdd_t));
    }

    *reachableStates = *fromUpperBound;
    /* throw away interior state candidates since we just threw away
     * the new states
     */
    if (hdInfo->interiorStateCandidate != NIL(mdd_t)) {
      mdd_free(hdInfo->interiorStateCandidate);
      hdInfo->interiorStateCandidate = NIL(mdd_t);
    }

    if (verbosityLevel >= 2) {
      /* record number of dead ends */
      (hdInfo->numDeadEnds)++;
      fprintf(vis_stdout, "Dead-End %d at %g seconds - FULL\n",
              hdInfo->numDeadEnds, (util_cpu_time()/1000.0));
    }

    /* perform a full non-greedy dead-end computation */
    *fromLowerBound = FsmHdDeadEnd(mddManager, imageInfo, *reachableStates,
                                  hdInfo,
                                  mintermVarArray,
                                  FSM_HD_NONGREEDY, verbosityLevel);
    (*depth)++;
    /* update reached since a new set of seeds have been found */
    assert (!mdd_is_tautology(*fromLowerBound, 0));
    if (shellFlag && onionRings) {
      mdd_t *temp = mdd_and(*fromLowerBound, *reachableStates, 1, 0);
      if (!mdd_is_tautology(temp, 0))
        array_insert(mdd_t *, onionRings, array_n(onionRings)-1, temp);
      else
        mdd_free(temp);
    }
    *fromUpperBound = *reachableStates;
    *reachableStates = mdd_or(*fromLowerBound, *fromUpperBound, 1, 1);
    FsmHdStatsComputeSizeAndMinterms(mddManager, *reachableStates,
                                     mintermVarArray,  nvars,
                                     Fsm_Hd_Reached, hdInfo->stats);
  } /* end of if jump in Reached */
} /* end of  HdInduceFullDeadEndIfNecessary */


static mdd_t *
ComputeInitialStates(Fsm_Fsm_t *fsm, int shellFlag, boolean printWarning)
{
  mdd_t *initialStates = NIL(mdd_t);
  boolean rings = FALSE;
  array_t *mintermVarArray = Fsm_FsmReadPresentStateVars(fsm);

  if (Fsm_FsmReadCurrentReachableStates(fsm) == NIL(mdd_t)) {
    initialStates = Fsm_FsmComputeInitialStates(fsm);
    return (initialStates);
  } else if (shellFlag) {
    (void) Fsm_FsmReachabilityOnionRingsStates(fsm, &initialStates);
    if (initialStates == NIL(mdd_t)) {
      initialStates = Fsm_FsmComputeInitialStates(fsm);
      return (initialStates);
    } else if (mdd_is_tautology(initialStates, 0)) {
      initialStates = Fsm_FsmComputeInitialStates(fsm);
      return (initialStates);
    }
    /* some states are found in the onion rings */
    rings = TRUE;
  } else {
    /* use previously computed states */
    initialStates = mdd_dup(Fsm_FsmReadCurrentReachableStates(fsm));
  }
  if (printWarning) {
    fprintf(vis_stdout, "** rch warning: Starting reachability analysis from previously computed states.\n");
    fprintf(vis_stdout, "** rch warning: Use compute_reach -F to recompute if necessary.\n");
    fprintf(vis_stdout, "** rch warning: Previously computed states = %g, Size = %d, computation depth = 0-%d\n", mdd_count_onset(Fsm_FsmReadMddManager(fsm), initialStates, mintermVarArray), mdd_size(initialStates), fsm->reachabilityInfo.depth);
    if (!rings) {
      if (FsmReadReachabilityComputationMode(fsm) == Fsm_Rch_Bfs_c) {
        fprintf(vis_stdout, "** rch warning: Previous computation done with BFS (-A 0 option)\n");
      } else if (FsmReadReachabilityComputationMode(fsm) == Fsm_Rch_Hd_c) {
        fprintf(vis_stdout, "** rch warning: Some previous computations done using BFS/DFS mode (-A 1 or -g option)\n");
      }
    } else {
      if (Fsm_FsmReadReachabilityOnionRingsMode(fsm) == Fsm_Rch_Bfs_c) {
        fprintf(vis_stdout, "** rch warning: Previous onion rings computed with BFS (-A 0 option)\n");
      } else if (Fsm_FsmReadReachabilityOnionRingsMode(fsm) == Fsm_Rch_Hd_c) {
        fprintf(vis_stdout, "** rch warning: Some set of onion rings have been computed using BFS/DFS (-A 1 or -g option)\n");
      }
    }
  }
  return (initialStates);
} /* end of ComputeInitialStates */


static int
CheckStatesContainedInInvariant(mdd_t *reachableStates,
                                array_t *invarStates)
{
  int i;
  mdd_t *invariant;

  arrayForEachItem(mdd_t *, invarStates, i, invariant) {
    if (invariant == NIL(mdd_t)) continue;
    if (!mdd_lequal(reachableStates, invariant, 1, 1)) {
      return 0;
    }
  }
  return 1;
}

static void
HdComputeReachabilityParameters(mdd_manager *mddManager,
                                Img_ImageInfo_t *imageInfo,
                                mdd_t **reachableStates,
                                mdd_t **fromUpperBound,
                                mdd_t **fromLowerBound,
                                mdd_t **image,
                                FsmHdStruct_t *hdInfo,
                                mdd_t *initialStates,
                                array_t *mintermVarArray,
                                int *depth,
                                int printStep,
                                int shellFlag,
                                array_t *onionRings,
                                int verbosityLevel)
{
  float growth;
  boolean scrapStates;
  int nvars = Fsm_FsmHdOptionReadNumVars(hdInfo->options);
  int oldSize;
  double oldMint;
  int sizeOfOnionRings = 0;
  if (shellFlag && onionRings) {
    sizeOfOnionRings = array_n(onionRings);
  }

  verbosityLevel = (printStep && (((*depth) % printStep) == 0)) ? verbosityLevel : 0;

  /* keep image if small, record size of the image */
  oldSize = FsmHdStatsReadSizeReached(hdInfo->stats);
  oldMint = FsmHdStatsReadMintermReached(hdInfo->stats);
  FsmHdStatsReset(hdInfo->stats, Fsm_Hd_Reached);
  FsmHdStatsReset(hdInfo->stats, Fsm_Hd_FromSubset);
  FsmHdStatsReset(hdInfo->stats, Fsm_Hd_From);

  scrapStates =  Fsm_FsmHdOptionReadScrapStates(hdInfo->options);
  hdInfo->onlyPartialImage =  Fsm_FsmHdOptionReadOnlyPartialImage(hdInfo->options);

  /* Step 1: Check comments at the top of this function */
  /* Check if this is a dead end */
  /* if growth rate of reached is really small, induce a dead-end */
  hdInfo->deadEndComputation = FALSE;
  if ((hdInfo->imageOfReachedComputed) ||
      (mdd_is_tautology(*fromLowerBound, 0)) ||
      (hdInfo->slowGrowth >= 5)) {
    /*
     * if image of reached is computed, add to the dead states.
     * Set dead-end computation flag and add all new states to Reached.
     */

    /* add to interior states */
    hdInfo->interiorStates = AddStates(hdInfo->interiorStates,
                                       hdInfo->interiorStateCandidate,
                                       FSM_MDD_FREE, FSM_MDD_FREE);
    hdInfo->interiorStateCandidate = NIL(mdd_t);
    if ((verbosityLevel >= 2) && (hdInfo->interiorStates != NIL(mdd_t))) {
      fprintf(vis_stdout, "Interior states added, Size = %d, Minterms = %g\n",
              mdd_size(hdInfo->interiorStates), mdd_count_onset(mddManager,
                                        hdInfo->interiorStates, mintermVarArray));
    }

    hdInfo->deadEndComputation = TRUE;
    hdInfo->firstSubset = FALSE;
    hdInfo->slowGrowth = 0;
    hdInfo->lastIter = *depth;

    if (!hdInfo->imageOfReachedComputed) {
      /* Either slow growth or real dead-end */
      if (verbosityLevel >= 2) {
        if (!mdd_is_tautology(*fromLowerBound, 0))
          fprintf(vis_stdout, "Dead-End triggered by slow growth\n");
        /* record number of dead ends */
        (hdInfo->numDeadEnds)++;
        fprintf(vis_stdout, "Dead-End %d at %g seconds\n", hdInfo->numDeadEnds,
                (util_cpu_time()/1000.0));
      }

      /* add another onion ring. */
      if (shellFlag && onionRings && !mdd_is_tautology(*fromLowerBound, 0)) {
        mdd_t *temp = mdd_and(*fromLowerBound, *reachableStates, 1, 0);
        if (!mdd_is_tautology(temp, 0))
          array_insert_last(mdd_t *, onionRings, temp);
        else
          mdd_free(temp);
      }
      *reachableStates = AddStates(*fromLowerBound, *reachableStates,
                                  FSM_MDD_FREE, FSM_MDD_FREE);
      *fromUpperBound = *reachableStates;

          /* perform a dead-end computation */
      *fromLowerBound = FsmHdDeadEnd(mddManager, imageInfo, *reachableStates,
                                    hdInfo,
                                    mintermVarArray,
                                    FSM_HD_GREEDY, verbosityLevel);

      (*depth)++;
    }
  } /* end of if dead-end computation */

  /* Step 2: Check comments at the top of this function */
  if (!mdd_is_tautology(*fromLowerBound, 0)) {
    /* add another onion ring. */
    if (shellFlag && onionRings) {
      mdd_t *temp = mdd_and(*fromLowerBound, *reachableStates, 1, 0);
      if (!mdd_is_tautology(temp, 0))
        array_insert_last(mdd_t *, onionRings, temp);
      else
        mdd_free(temp);
    }
    /* compute reachable states with new states */
    *fromUpperBound = *reachableStates;
    *reachableStates = mdd_or(*fromLowerBound, *fromUpperBound, 1, 1);
    FsmHdStatsComputeSizeAndMinterms(mddManager, *reachableStates,
                                     mintermVarArray, nvars,
                                     Fsm_Hd_Reached, hdInfo->stats);


    /*
         * induce a full dead-end computation if jump in size of Reached is
         * large. But when the image of Reached has been computed or a dead-end
         * computation was just done, this is not done because deadEndComputation
         * has been set to 1.
         */
    HdInduceFullDeadEndIfNecessary(mddManager, imageInfo, fromLowerBound,
                                   fromUpperBound, reachableStates, image,
                                   hdInfo, depth, shellFlag, onionRings,
                                   mintermVarArray, oldSize,
                                   oldMint, verbosityLevel);
  } /* end of  if (!mdd_is_tautology(fromLowerBound, 0)) */

  /* Step 3: Check comments at the top of this function */
  /* Create subset if necessary */
  if (!mdd_is_tautology(*fromLowerBound, 0)) {
    if (!hdInfo->onlyPartialImage) {
      mdd_t *newStates;
      /*
       * Now fromUpperBound has the old reached states, reachableStates
       * has the new set of states.
       */
      newStates = mdd_dup(*fromLowerBound);
      FsmHdFromComputeDenseSubset(mddManager, fromLowerBound,
                                  fromUpperBound, reachableStates,
                                  image, initialStates, hdInfo,
                                  shellFlag, onionRings, sizeOfOnionRings,
                                  mintermVarArray);

      /* if all new states are chosen from the previous iteration and
       * all these states have been added then add to interior states.
       * All new states are sometimes not added when
       * hdInfo->options->scrapStates is 0. Dead-end computations need not be
       * considered here since the interior states were already added then.
       */
      if ((hdInfo->interiorStateCandidate != NIL(mdd_t)) &&
          ((mdd_lequal(newStates, *fromLowerBound, 1, 1)) ||
           (scrapStates == TRUE))) {
        hdInfo->interiorStates = AddStates(hdInfo->interiorStates, hdInfo->interiorStateCandidate,
                                   FSM_MDD_FREE, FSM_MDD_FREE);
        hdInfo->interiorStateCandidate = NIL(mdd_t);
        if (verbosityLevel >= 2) {
          fprintf(vis_stdout,
                  "Interior states added, Size = %d, Minterms = %g\n",
                  mdd_size(hdInfo->interiorStates),
                  mdd_count_onset(mddManager, hdInfo->interiorStates,
                                  mintermVarArray));
        }
      }
      mdd_free(newStates);

      /* end of if fromLowerbound == 0 && subset Flag*/
    } else  {
      mdd_free(*fromUpperBound);
      if (printStep && ((*depth) % printStep == 0)) {
        FsmHdStatsComputeSizeAndMinterms(mddManager, *fromLowerBound,
                                         mintermVarArray, nvars, Fsm_Hd_From, hdInfo->stats);
        FsmHdStatsSetMintermFromSubset(hdInfo->stats, 0.0);
        FsmHdStatsSetSizeFromSubset(hdInfo->stats, 0);
      }
      *fromUpperBound = *reachableStates;
    }
  } /* end of if fromLowerBound == 0 */
  if (hdInfo->interiorStateCandidate != NIL(mdd_t)) {
    mdd_free(hdInfo->interiorStateCandidate);
    hdInfo->interiorStateCandidate = NIL(mdd_t);
  }

  /* monitor growth rate of reached, record if slow growth in consecutive
   * iterations
   */
  if (!mdd_is_tautology(*fromLowerBound, 0)) {
    FsmHdStatsComputeSizeAndMinterms(mddManager, *reachableStates,
                                     mintermVarArray, nvars, Fsm_Hd_Reached, hdInfo->stats);
    growth = (float)(FsmHdStatsReadMintermReached(hdInfo->stats) - oldMint)/
      (float)oldMint;
    if ((growth < FSM_HD_GROWTH_RATE) &&
        (((hdInfo->slowGrowth) && ((*depth) == (hdInfo->lastIter) + 1)) || (!(hdInfo->slowGrowth)))) {
      if (verbosityLevel >= 2) {
        fprintf(vis_stdout, "Growth rate = %f\n", growth);
      }
      (hdInfo->slowGrowth)++;
      hdInfo->lastIter = *depth;
    } else {
      /* reset value */
      hdInfo->slowGrowth = 0;
    }
  }
  return;
} /* end of HdComputeReachabilityParameters */

static int
InitializeIncrementalParameters(
  Fsm_Fsm_t *fsm,
  Ntk_Network_t *network,
  array_t **arrayOfRoots,
  st_table **tableOfLeaves)
{
  int numLatch;
  Ntk_Node_t *node;
  lsGen gen;
  int i, mddId;  /* iterate over the nextStateVars of the fsm */

  numLatch = Ntk_NetworkReadNumLatches(network);

  /* This changed to fix a bug.  Instead of reading the roots using
     Ntk_NetworkForEachCombOutput, we now get them from
     fsm->fsdData.nextStateVars, which means that the order is consistent with
     the latter field, as required later on in the incremental reachability
     code

     old code:
     *arrayOfRoots = array_alloc(Ntk_Node_t *, numLatch);
     Ntk_NetworkForEachCombOutput(network, gen, node) {
     if(Ntk_NodeTestIsLatchDataInput(node)){
     array_insert_last(Ntk_Node_t *, *arrayOfRoots, node);
    }
  }
  */

  assert(numLatch == array_n(fsm->fsmData.nextStateVars));
  *arrayOfRoots = array_alloc(Ntk_Node_t *, numLatch);

  arrayForEachItem(int, fsm->fsmData.nextStateVars, i, mddId){
    Ntk_Node_t *shadow;
    Ntk_Node_t *latch;
    Ntk_Node_t *dataInput;

    shadow = Ntk_NetworkFindNodeByMddId(network, mddId);
    assert(Ntk_NodeTestIsShadow(shadow));
    latch = Ntk_ShadowReadOrigin(shadow);
    assert(Ntk_NodeTestIsLatch(latch));
    dataInput = Ntk_LatchReadDataInput(latch);
    assert(dataInput != NIL(Ntk_Node_t));

    array_insert(Ntk_Node_t *, *arrayOfRoots, i, dataInput);
  }

  *tableOfLeaves = st_init_table(st_ptrcmp, st_ptrhash);
  Ntk_NetworkForEachCombInput(network, gen, node) {
    st_insert(*tableOfLeaves, (char *)node, (char *) (long) (-1) );
  }
  return (numLatch);
}


static mdd_t *
IncrementalImageCompute(Ntk_Network_t *network,
                        Fsm_Fsm_t *fsm,
                        mdd_manager *mddManager,
                        mdd_t *fromLowerBound,
                        mdd_t *fromUpperBound,
                        mdd_t *toCareSet,
                        array_t *arrayOfRoots,
                        st_table *tableOfLeaves,
                        array_t *smoothVarArray,
                        array_t *relationArray,
                        int numLatch)
{
  Mvf_Function_t *function;
  mdd_t *relation, *optimizedRelation, *toCareSetRV, *imageRV;
  mdd_t  *subOptimizedRelation, *image;
  int mddId, i;
  array_t *arrayOfMvf;
  mdd_t *frontierSet;

  frontierSet = bdd_between(fromLowerBound, fromUpperBound);
  /* Create the array of transition relations */
  toCareSetRV = mdd_substitute(mddManager, toCareSet,
                               fsm->fsmData.presentStateVars,
                               fsm->fsmData.nextStateVars);

  arrayOfMvf = Ntm_NetworkBuildMvfs(network, arrayOfRoots,
                                    tableOfLeaves, frontierSet);
  for (i=0; i < numLatch; i++){
    function = array_fetch(Mvf_Function_t *, arrayOfMvf, i);
    mddId = array_fetch(int, fsm->fsmData.nextStateVars, i);
    /* Since both arrayOfMvf and fsmData have been obtained from */
    /* the same generator */
    relation = Mvf_FunctionBuildRelationWithVariable(function,
                                                     mddId);
    subOptimizedRelation = bdd_cofactor(relation, toCareSetRV);
    mdd_free(relation);
    optimizedRelation = bdd_cofactor(subOptimizedRelation, frontierSet);
    mdd_free(subOptimizedRelation);
    array_insert(mdd_t *, relationArray, i, optimizedRelation);
  }
  Mvf_FunctionArrayFree(arrayOfMvf);
  mdd_free(toCareSetRV);
  array_insert(mdd_t *, relationArray, numLatch, frontierSet);
  imageRV = Fsm_MddMultiwayAndSmooth(mddManager, relationArray,
                                     smoothVarArray);
  /*
  ** The above call can be substituted by more sophisticated
  ** Img_MultiwayLinearAndSmooth. However, that will have its
  ** associated overhead, and could offset any advantage. We
  ** expect that individual transition relation relations should
  ** be small and monolithic way to handle them would be ok.
  ** However, a good strategy would be find the quantification
  ** schedule which does not change with the computation of
  ** incremental transition relations.
  */
  for (i = 0; i <= numLatch; i++){
    mdd_free(array_fetch(mdd_t*, relationArray, i));
  }
  assert(CheckImageValidity(mddManager, imageRV,
                            fsm->fsmData.presentStateVars,
                            fsm->fsmData.inputVars));
  image = mdd_substitute(mddManager, imageRV,
                         fsm->fsmData.nextStateVars,
                         fsm->fsmData.presentStateVars);
  mdd_free(imageRV);
  assert(CheckImageValidity(mddManager, image,
                            fsm->fsmData.nextStateVars,
                            fsm->fsmData.inputVars));
  return image;
}

static void
PrintOnionRings(Fsm_Fsm_t *fsm, int printStep,
                Fsm_RchType_t approxFlag, int nvars)
{
  int i;
  mdd_t *reachableStates;
  mdd_manager *mddManager = Fsm_FsmReadMddManager(fsm);
  array_t *mintermVarArray = Fsm_FsmReadPresentStateVars(fsm);

  fprintf(vis_stdout, "** rch warning: Not possible to compute size of reachable states at every iteration\n");
  if (Fsm_FsmReadReachabilityOnionRings(fsm) != NIL(array_t)) {
    arrayForEachItem(mdd_t *, Fsm_FsmReadReachabilityOnionRings(fsm),
                     i, reachableStates) {
      if (printStep && (i % printStep == 0)) {
        if (approxFlag != Fsm_Rch_Hd_c) {
          if (nvars <= 1023) {
            (void)fprintf(vis_stdout, "BFS step = %d\t|onion ring states| = %8g",
                          i, mdd_count_onset(mddManager, reachableStates,
                                             mintermVarArray));
            if (i != array_n(Fsm_FsmReadReachabilityOnionRings(fsm))-1) {
              fprintf(vis_stdout, "\n");
            } else {
              fprintf(vis_stdout, "\tSize = %d\n", mdd_size(reachableStates));
            }

          } else {
            (void)fprintf(vis_stdout, "BFS step = %d\t|onion ring states| = ",
                          i);
            (void) bdd_print_apa_minterm(vis_stdout, reachableStates, nvars, (long)6);
            if (i != array_n(Fsm_FsmReadReachabilityOnionRings(fsm))-1) {
              fprintf(vis_stdout, "Final Size = %d\n", mdd_size(reachableStates));
            }
          }
        } else {

          if (nvars <= 1023) {
            (void)fprintf(vis_stdout, "Step = %d\t|onion ring states| = %8g",
                          i, mdd_count_onset(mddManager, reachableStates,
                                             mintermVarArray));
            if (i != array_n(Fsm_FsmReadReachabilityOnionRings(fsm))-1) {
              fprintf(vis_stdout, "\n");
            } else {
              fprintf(vis_stdout, "\tSize = %d\n", mdd_size(reachableStates));
            }
          } else {
            (void)fprintf(vis_stdout, "Step = %d\t|onion ring states| = ",
                          i);
            (void) bdd_print_apa_minterm(vis_stdout, reachableStates, nvars, (long)6);
            if (i != array_n(Fsm_FsmReadReachabilityOnionRings(fsm))-1) {
              fprintf(vis_stdout, "Final Size = %d\n", mdd_size(reachableStates));
            }
          }
        }
      }
    }
  }
  return;
}/* printonionrings */

static array_t *
GenerateGuidedSearchSequenceArray(array_t *guideArray)
{
  char *seqStr = Cmd_FlagReadByName("guided_search_sequence");
  char *preStr;
  int intEntry;
  array_t *depthArray = array_alloc(int, 0);
  int i;

  /* Skip parsing if we have no guided search sequence  */
  if (seqStr != NIL(char)  && strcmp(seqStr, "") != 0){
    preStr = strtok(seqStr, ",");

    while(preStr != NIL(char)){
      intEntry = strtol(preStr, (char **) NULL, 10);
      array_insert_last(int, depthArray, (intEntry == 0 ? -1 : intEntry));
      preStr = strtok(NIL(char), ",");
    }
  }

  if(array_n(depthArray) > array_n(guideArray)){
    fprintf(vis_stdout, "** rch warning: guided_search_sequence has more entries\n");
    fprintf(vis_stdout,"than there are hints.  Extra entires will be ignored\n");
  }

  /* pad with -1s */
  for (i = array_n(depthArray); i < array_n(guideArray); i++) {
    array_insert_last(int, depthArray, -1);
  }

  assert(array_n(depthArray) >= array_n(guideArray));

  return depthArray;
}


static void
ComputeReachabilityParameters(mdd_manager *mddManager,
                              Img_ImageInfo_t *imageInfo,
                              Fsm_RchType_t approxFlag,
                              mdd_t **reachableStates,
                              mdd_t **fromUpperBound,
                              mdd_t **fromLowerBound,
                              mdd_t **image,
                              FsmHdStruct_t *hdInfo,
                              mdd_t *initialStates,
                              array_t *mintermVarArray,
                              int *depth,
                              int printStep,
                              int shellFlag,
                              array_t *onionRings,
                              int verbosityLevel,
                              boolean guidedSearchMode,
                              mdd_t *hint,
                              int *hintDepth,
                              boolean *notConverged)
{
  if (approxFlag == Fsm_Rch_Bfs_c) {
    /* New set of reachable states is old set plus new states in image. */
    *reachableStates = AddStates(*fromLowerBound, *reachableStates,
                                 FSM_MDD_DONT_FREE, FSM_MDD_FREE);
    *fromUpperBound = *reachableStates;

    /* in guided search, if the mode is to convergence or the
       limited depth has not been reached, and the frontier is empty,
       then generate a frontier in BFS mode (that is a full dead-end).  */
    if (guidedSearchMode &&
        mdd_is_tautology(*fromLowerBound, 0) /* no frontier */ &&
        (hdInfo->deadEndWithOriginalTR == FALSE) /* and guided search */ ) {
      assert(imageInfo != NIL(Img_ImageInfo_t));
      if (*hintDepth != 0) {
        /* add another onion ring. */
        if (shellFlag && onionRings) {
          mdd_t *temp = mdd_and(*fromLowerBound, *reachableStates, 1, 0);
          if (!mdd_is_tautology(temp, 0))
            array_insert_last(mdd_t *, onionRings, temp);
          else
            mdd_free(temp);
        }
        /* do a full dead-end (non-greedy), compute all the frontier states */
        mdd_free(*fromLowerBound);
        *fromLowerBound = FsmHdDeadEnd(mddManager, imageInfo, *reachableStates,
                                       hdInfo,
                                       mintermVarArray,
                                       FSM_HD_NONGREEDY, verbosityLevel);
        mdd_free(hdInfo->interiorStates);
        hdInfo->interiorStates = NIL(mdd_t);
        *reachableStates = AddStates(*fromLowerBound, *reachableStates,
                                     FSM_MDD_DONT_FREE, FSM_MDD_FREE);
        *fromUpperBound = *reachableStates;
        (*hintDepth)--;
        (*depth)++;
        if (mdd_is_tautology(hint, 1)) hdInfo->deadEndWithOriginalTR = TRUE;
      } else {
        *notConverged = TRUE;
        /* set this flag for correct detection of convergence
           Corner case: when limited depth and frontier = 0 */
      }
    }
  } else if (approxFlag == Fsm_Rch_Hd_c) {
    assert(imageInfo != NIL(Img_ImageInfo_t));
      if (!guidedSearchMode || (*hintDepth != 0)) {
        /* computes reachable states, fromLowerBound and fromUpperBound.
         * Updates all other quantities such as interiorStates,
         * deadEndResidue. */
        HdComputeReachabilityParameters(mddManager, imageInfo,
                                        reachableStates,
                                        fromUpperBound,
                                        fromLowerBound,
                                        image,
                                        hdInfo,
                                        initialStates,
                                        mintermVarArray,
                                        depth, printStep,
                                        shellFlag, onionRings,
                                        verbosityLevel);
        if (guidedSearchMode) (*hintDepth)--;

      } else { /* guidedSearchMode and hint depth is 0 */
        if (mdd_is_tautology(*fromLowerBound, 0)) {
          *notConverged = TRUE;
        /* set this flag for correct detection of convergence
           Corner case: when limited depth and frontier = 0 */
        }
        /* add another onion ring. */
        if (shellFlag && onionRings) {
          mdd_t *temp = mdd_and(*fromLowerBound, *reachableStates, 1, 0);
          if (!mdd_is_tautology(temp, 0))
            array_insert_last(mdd_t *, onionRings, temp);
          else
            mdd_free(temp);
        }

      }
  }
  return;
} /* End of ComputeReachabilityParameters */