source: vis_dev/vis-2.3/src/fsm/fsmReach.c @ 45

/**CFile***********************************************************************

  FileName    [fsmReach.c]

  PackageName [fsm]

  Synopsis    [Routines to perform reachability on the FSM structure.]

  Author      [Rajeev K. Ranjan, In-Ho Moon, Kavita Ravi]

  Copyright   [Copyright (c) 1994-1996 The Regents of the Univ. of California.
  All rights reserved.

  Permission is hereby granted, without written agreement and without license
  or royalty fees, to use, copy, modify, and distribute this software and its
  documentation for any purpose, provided that the above copyright notice and
  the following two paragraphs appear in all copies of this software.

  IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY FOR
  DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
  OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
  CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

  THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
  FITNESS FOR A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN
  "AS IS" BASIS, AND THE UNIVERSITY OF CALIFORNIA HAS NO OBLIGATION TO PROVIDE
  MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.]

******************************************************************************/

#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                                                    */
/*---------------------------------------------------------------------------*/

/**AutomaticStart*************************************************************/

/*---------------------------------------------------------------------------*/
/* 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);

/**AutomaticEnd***************************************************************/


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

  Synopsis [Computes the set of reachable states of an FSM.]

  Description [Computes the set of reachable states of an FSM.  If
  this set has already been computed, then just returns the result of
  the previous computation. The recomputation is done anyway if onion
  rings are required and have not been previously computed. The
  calling application is responsible for freeing the returned MDD in
  all cases.  If verbosityLevel is greater than 1, then information is
  printed after every k steps of reachability, where k is defined by
  the -s option in the compute_reach command (the default for k is
  1). If reorderFlag is set, then dynamic reordering is invoked once
  when the size of the reachable state set reaches the threshold given
  by reorderThreshold.  This is independent of the
  dynamic_variable_reordering command. The shell flag indicates that
  the onion rings of frontier states should be stored. The default
  does not store the onion rings. The default value for approxFlag is
  Fsm_Rch_Default_c. See fsm.h for values.

  The procedure offers the following variations of reachability
  analysis: 1. Exact Reachability analysis by breadth-first search
  (default). 2. Reachability Analysis (HD) . 3. Reachability Analysis
  (upper bound by overapproximate FSM traversal). 4 ARDC-accelerated
  Reachability Analysis.  5. Invariant checking with early termination.

  A lower bound (specified by depthValue) on the reachable states may
  be obtained by performing a specified number of iterations
  (depthValue argument). It updates a field in the
  fsm.reachabilityInfo signalling that the reachable states are an
  under approximation. The default for the depthValue argument is 0
  which indicates that reachability analysis computation should
  proceed to a fixpoint. If a depthValue is specified after the
  computation of reachable states, the reachable states is returned,
  else depthValue number of computations are performed. Every non-zero
  depthValue is interpreted as the number of iterations (image
  computations) to be performed from the current state (not from the
  initial state).

  Reachability Analysis (HD) is requested as approxFlag =  Fsm_Rch_Hd_c.
  HD provides an exact set of reachable states with memory efficiency or
  a lower bound if reachability analysis does not complete.
  HD Principles: 1. compute image of small BDDs. 2. compute partial
  image if necessary 3. When no new states are generated (dead-end),
  compute new states from image of reached in a decomposed form.  Allow
  no partial image.  4. When there is a big jump in the size of
  reached, compute its entire image.  5. store states that will
  produce no new successors (interior states)

  The algorithm is implemented as follows
  from = init; reached = init;
  while (reached != 1 && from != 0) {
    Allow partial image;
    image = Image(from);
    Check if partial image;
    from = image - reached;
    if (!partial image) interior_states_candidate = from;
    if (from == 0 || slow growth of Reached for 5 iterations) {
      Dead-end computation(greedy): find new states from image of reached;
    }
    if (from != 0 and (reached+from) size has a big jump) {
      Dead-end computation(non-greedy): find the entire image of reached;
    }
    if (from != 0) {
      Compute dense subset of from with at least 1 new state;
      if subset contains from (i.e., all new states)
                       interior states += interior states candidate;
      reached = reached + from/subset;
    }
  }

  Legend:  fromLowerBound = from; initialStates = init;
           reachableStates = reached
           fromUpperBound = reachableStates (BFS)/fromLowerBound(HD)
           interiorStates = states with no new successors
           interiorStateCandidates = states that may be interior
                  states depending on whether all of them are added to
                  reachableStates.

  Currently, this does not support the incrementalFlag = 1 option.

  Reachability Analysis (upper bound by overapproximate FSM traversal)
  gives an upper bound of exact reachable states by overapproximation.
  It is requested as approxFlag = Fsm_Rch_Oa_c. To do this, first, it
  groups latch variables by using SSD(State Space Decomposition in the
  paper Cho et.al, TCAD96-DEC) algorithm, then computes the reachable
  states of each submachines until convergence, then the product of
  reachable states of all submachines is an overapproximate reachable
  states. In this case, depthValue and shellFlag and incrementalFlag
  are incompatible with this option.

  ARDC-accelerated Reachability Analysis may allow faster and smaller
  memory usage reachability analysis. First, it computes
  overapproximated reachable states (if already computed, just use it).
  Reachability Analysis is then performed (any value of approxFlag is
  applicable) with a transition relation minimized with the complement
  of overapproximate reachable states as don't cares (called
  ARDC). This option is requested as ardc = 1.

  invarStates is a set of invariant states that can be specified. The
  reachable states must satisfy all invariants. The default value is
  NULL. Any invalid invariant must also be set to NULL in the
  invarStaes array. This option is used by the check_invariant
  command. When an array of invariant states is provided, reachability
  analysis is done until an invariant is violated or when reachability
  finishes, whichever is earlier. If an approximation of the reachable
  states exists, then each invariant is checked for violation. If none
  are violated, the set of reachable states is returned. It is the
  caller's resposibility to remove the violating invariant from the
  invarStates array to continue checking the remaining
  invariants. The presence of invarStates takes precedence over a
  non-zero depthValue.

  In general, the reachable states may not be consistent with the
  onion rings i.e., sum of onion rings is always contained in the
  reachable states.  When the shell flag is specified and onion rings
  exist, the initial states are set to the sum of the onion ring
  states instead of the default initial states of the fsm or
  previously computed reachable states.

  Some random simulation of the reachable states can be done prior to
  reachability analysis. This is turned on using the rch_simulate
  flag. A number should be specifed in the command : set sim_reach
  <number>, which specifies the number of vectors to be simulated.
  The initial states are then set to the simulated set of states. If
  HD is specified, a subset may be taken.

  The recompute flag asks for an explicit recomputation. This means
  that even reachable states have been computed, they will be
  recomputed.

  On consecutive calls, recomputation is avoided as much as
  possible. The specific cases are when
  Case a. No shell Flag: Recomputation is avoided when reachability is
  done, when an underapproximation exists and it violates an
  invariant,  an overapproximation exists and all invariants pass.

  Case b. With shell flag: Recomputation is avoided when onion rings
  are up-to-date and reachability is done or when onion rings are
  up-to-date, an underapproximation exists and an invariant fails.

  Update of fsm.reachabiltiyInfo fields are also based on minimum
  recomputation. The reachability field is updated only if the current
  set of reachable states are larger than the previous set.

  printWarning flag turns is used to turn on/off the warning message
  printed when the previously computed reachable states are used to
  proceed with computation. The default is TRUE.  ]

  SideEffects [The result of the reachable state computation is stored
  with the FSM. i.e., depth, reachable set, underapprox flag. ]

  SeeAlso [Fsm_FsmComputeInitialStates
  Fsm_ArdcComputeOverApproximateReachableStates]

******************************************************************************/
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 */



/**Function********************************************************************

  Synopsis [Existentially quantifies variables from an implicit product.]

  Description [Existentially quantifies the smoothing variables from the
  product of the MDDs in mddArray, and returns the result. SmoothingVars is an
  array of MDD ids. For now, uses the simplistic approach of conjunction all
  the MDDs, and then quantifying the smoothing variables.]

  SideEffects []

******************************************************************************/
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;
}

/**Function********************************************************************

  Synopsis    [Prints the results of reachability analysis.]

  Description    [Prints the results of reachability analysis.]

  SideEffects []

  SeeAlso     []

******************************************************************************/
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                                          */
/*---------------------------------------------------------------------------*/

/**Function********************************************************************

  Synopsis    [Checks the validity of image]

  Description [In a properly formed image, there should not be any
  domain or quantify variables in its support. This function checks
  for that fact.]

  SideEffects []

******************************************************************************/
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;
}

/**Function********************************************************************

  Synopsis    [Counts the number of Bdd variables for the given Mdd id array].

  Description [Counts the number of Bdd variables for the given Mdd id array]


  SideEffects []

******************************************************************************/
/*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 */


/**Function********************************************************************

  Synopsis    [Add states.]

  Description [Add States.]

  SideEffects [ Free the original parts]

******************************************************************************/
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 */


/**Function********************************************************************

  Synopsis    [Simulates randomly specified number of states.]

  Description [Simulates randomly specified number of states.]

  SideEffects []

******************************************************************************/
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 */


/**Function********************************************************************

  Synopsis    [Prints stats per iteration]

  Description [Prints stats per iteration]

  SideEffects []

******************************************************************************/
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 */

/**Function********************************************************************

  Synopsis [Induce full dead-end if there is a big jump in Reached
  size]

  Description [Induce full dead-end if there is a big jump in Reached
  size and not enough increase in minterms of Reached. So throw away
  the new large-sized Reached and compute its entire image. Use the
  dead-end computation for this and use the non-greedy flag. This is a
  recovery mechanism to combat fragmentation of the state space.]

  SideEffects [Other quantities are adjusted accordingly. Now the Reached
  set before the dead-end computation will be added to the interior states.]

******************************************************************************/
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 */


/**Function********************************************************************

  Synopsis [Compute initial states.]

  Description [Compute initial states. Return fsm initial states if no
  previously computed states are found. If shell flag, compute union
  of rings. If not shell flag, return previously computed states.]

  SideEffects []

******************************************************************************/
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 */


/**Function********************************************************************

  Synopsis [Checks if the reachable states are contained in each
  element of an array of invariants.]

  Description [Checks if the reachable states are contained in each
  element of an array of invariants. Returns 0 as soon as a violation
  is found.]

  SideEffects []

******************************************************************************/
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;
}

/**Function********************************************************************

  Synopsis [Computes the reachable states, fromLowerBound and
  fromUpperBound for HD Traversal.]

  Description [Computes the reachable states, fromLowerBound and
  fromUpperBound for HD Traversal. Updates all other arguments such as
  interiorStates, deadEndResidue . There are 3 major operations in
  this procedure. 1. Induce (if slow growth in reached) and compute a
  dead-end if necessary. 2. Compute the image of Reached if there is a
  big jump in its size. 3. Compute a dense subset of states to be used
  as a frontier in the next iteration.

  1. A dead-end is the situation where no new states have been
  produced (fromLowerBound = 0), It is induced if the rate of growth
  of Reached is slow. If the image of the Entire Reached set has been
  computed in the previous iteration, then pretend this is a dead-end.

  2. If the jump in Reached is large, compute its entire image to
  recover from fragmentation.

  3. A dense subset of the frontier set is computed to be used as the
  frontier set in the next iteration. The update of reached is
  performed in this step (may not add scrap states if
  specified.). If a dead-end has just been performed, the new states
  are always added to Reached. Growth rate monitors are also
  updated. ]

  SideEffects []

******************************************************************************/
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 */

/**Function********************************************************************

  Synopsis [Initializes parameters for incremental computation.]

  Description [Initializes parameters for incremental computation.
  Instantiates arrayOfRoots and tableOfLeafs, and returns the number of latches
  in the sytstem.  Upon return, arrayOfRoots is an array of Ntk_Node_t *, and
  tableOfLeaves is an st_table of Ntk_Node_t *.]

  SideEffects []

******************************************************************************/
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);
}


/**Function********************************************************************

  Synopsis [Computes image for incremental computation.]

  Description [Computes image for  incremental computation.]

  SideEffects []

******************************************************************************/
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;
}

/**Function********************************************************************

  Synopsis [Prints information about onion rings]

  Description [Prints information about onion rings.]

  SideEffects []

******************************************************************************/
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 */

/**Function********************************************************************

  Synopsis [Reads in a sequence of integers separated by commas and
  stores them in an array.]

  Description [Reads in a sequence of integers separated by commas from the
  environment variable guided_search_sequence.  It stores them in an array that
  it returns.  Alpha characters and 0 integers are replaced by
  -1.  Pads the remaining entries corresponding to the guide array with -1.
  Every number indicates the number of times the corresponding (in sequence)
  hint should be used.  The value -1 means `to convergence'.]


  SideEffects []

******************************************************************************/
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;
}


/**Function********************************************************************

  Synopsis [Computes the reachable states, frontier and fromUpperBound.]

  Description [Computes the reachable states, frontier and
  fromUpperBound. For HD it also generates interior states, interior
  state candidates, dead-end residue, number of dead-ends, residue
  count, etc. For guided search, this procedure also generates a
  dead-end if the frontier is 0. Also this procedure keeps track of
  convergence and whether a dead-end is required at the end.]

  SideEffects [All arguments are updated.]

  SeeAlso [HdComputeReachabilityParameters FsmHdDeadEnd]

******************************************************************************/
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 */