source: vis_dev/vis-2.3/src/img/imgIwls95.c @ 17

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

  FileName    [imgIwls95.c]

  PackageName [img]

  Synopsis    [Routines for image computation using component transition
               relation approach described in the proceedings of IWLS'95,
               (henceforth referred to Iwls95 technique).]

  Description [<p>The initialization process performs following steps:
  <UL>
  <LI> Create the relation of the roots at the bit level
  in terms of the quantify and domain variables.
  <LI> Order the bit level relations.
  <LI> Group the relations of bits together, making a cluster
  whenever the BDD size reaches a threshold.
  <LI> For each cluster, quantify out the quantify variables which
  are local to that particular cluster.
  <LI> Order the clusters using the algorithm given in
  "Efficient BDD Algorithms for FSM Synthesis and
  Verification", by R. K. Ranjan et. al. in the proceedings of
  IWLS'95{1}.
  <LI> The orders of the clusters for forward and backward image are
  calculated and stored. Also stored is the schedule of
  variables for early quantification.
  </UL>
  For forward and backward image computation the corresponding
  routines are called with appropriate ordering of clusters and
  early quantification schedule. <p>
  Note that, a cubic implementation of the above algorithm is
  straightforward. However, to obtain quadratic complexity requires
  significant amount of book-keeping. This is the reason for the complexity
  of the code in this file.]

  SeeAlso     []

  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 "imgInt.h"
#include "fsm.h"

static char rcsid[] UNUSED = "$Id: imgIwls95.c,v 1.127 2005/05/18 18:11:57 jinh Exp $";

/*---------------------------------------------------------------------------*/
/* Constant declarations                                                     */
/*---------------------------------------------------------------------------*/
#define domainVar_c 0
#define rangeVar_c 1
#define quantifyVar_c 2
#define IMG_IWLS95_DEBUG

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

typedef struct Iwls95InfoStruct Iwls95Info_t;
typedef struct CtrInfoStruct CtrInfo_t;
typedef struct VarInfoStruct VarInfo_t;
typedef struct CtrItemStruct CtrItem_t;
typedef struct VarItemStruct VarItem_t;

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

  Synopsis    [This structure stores the information needed for iterative image
  computation by IWLS'95 method.]

  Description [A step of image computation is computed by taking the product of
  "From" set with the BDD's of fwdClusteredRelationArray in order.
  "fwdArraySmoothVarBddArray[i]" is an array of Bdd variables, which represents
  the set of variables which can be quantified out once
  fwdClusteredRelationArray[i] is multiplied in the product. Similar analysis
  holds good for bwdClusteredRelationArray.

  <p>The "option" struct contains option values for computing the above
  mentioned arrays and also controls the verbosity.

  <p>The bit relation array contains the transition relation in a bit relation
  form.

  <p>The quantifyVarMddIdArray contains functionData->quantifyVars and any
  intermediate variables that need to be quantified.

  <p> The *SmoothVarCubeArrays hold the same information as the
  *ArraySmoothVarBddArrays, but in the form of an array of cubes, one cube for
  each entry of *ArraySmoothVarBddArrays, because the quantification functions
  expect a cube. 

  <p>AllowPartialImage is an option for high-density analysis, if it is true,
  image is allowed to return a subset instead of the correct result.  If it
  does so, it will set wasPartialImage.  The PIOptions field keeps the options
  for partial image computation. (See the definition of the type)

  <p> The bwdClusteredCofactoredRelationArray is the
  clusteredCofactoredRelationArray minimized wrt careSetArray.  Minimization
  only takes place when the careSetArray that is passed is changed.  They can
  both be NIL, which means that they need to be recomputed.
  
  <p> RB What are savedArray*, method, eliminateDepend, nFoundDependVars,
  averageFOundDependVars, nPrev*, and lazySiftFlag????]

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

struct Iwls95InfoStruct {
  array_t *bitRelationArray; 

  array_t *quantifyVarMddIdArray;

  array_t *fwdOptClusteredRelationArray;
  array_t *fwdClusteredRelationArray;
  array_t *fwdOriClusteredRelationArray;
  array_t *bwdClusteredRelationArray;
  array_t *bwdOriClusteredRelationArray;
  array_t *bwdClusteredCofactoredRelationArray;
  array_t *careSetArray;

  array_t *fwdOptArraySmoothVarBddArray;
  array_t *fwdArraySmoothVarBddArray;
  array_t *bwdArraySmoothVarBddArray;
  array_t *fwdOptSmoothVarCubeArray;
  array_t *fwdSmoothVarCubeArray;
  array_t *bwdSmoothVarCubeArray;

  boolean allowPartialImage;
  boolean wasPartialImage;
  ImgPartialImageOption_t *PIoption;

  array_t *savedArraySmoothVarBddArray;
  array_t *savedSmoothVarCubeArray;

  Img_MethodType method;
  ImgTrmOption_t *option;

  int eliminateDepend;
        /* 0 : do not perform this
           1 : eliminates dependent variables from
                constraint and modifies function vector by
                composition at every image computation.
           2 : once the number of eliminated variables
                becomes zero, do not perform.
        */
  int nFoundDependVars;
  float averageFoundDependVars;

  int nPrevEliminatedFwd;
  int linearComputeRange;

  boolean lazySiftFlag;
};


/**Struct**********************************************************************

  Synopsis    [This structure contains the information about a particular
               cluster.]

  Description [Corresponding to each cluster one ctrInfoStruct is created, The
               various fields of this structure contain the information for
               ordering the clusters for image computation
               purposes. "varItemList" is the linked list of "varItemStruct"
               corresponding to all the support variables of the cluster.]

******************************************************************************/
struct CtrInfoStruct {
  int ctrId; /* Unique id of the cluster */
  int numLocalSmoothVars; /*
                           * Number of smooth variables in the support of the
                           * cluster which can be quantified out if the cluster
                           * is multiplied in the current product.
                           */
  int numSmoothVars; /* Number of smooth variables in the support */
  int numIntroducedVars; /*
                          * Number of range variables introduced if the cluster
                          * is multiplied in the current product.
                          */
  int maxSmoothVarIndex; /* Max. index of a smooth variable in the support */
  int rankMaxSmoothVarIndex; /* Rank of maximum index amongst all ctr's */
  lsList varItemList; /* List of varItemStruct of the support variables */
  lsHandle ctrInfoListHandle; /*
                               * Handle to the list of info struct of ctr's
                               * which are yet to be ordered.
                               */
};



/**Struct**********************************************************************

  Synopsis    [This structure contains information about a particular variable
               (domain, range or quantify).]

  Description [Corresponding to each variable a structure is created. Like the
               varItemList field of ctrInfoStruct, ctrItemList contains the
               list of ctrInfoStruct corresponding to various clusters which
               depend on this variable.]

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

struct VarInfoStruct {
  int bddId;  /* BDD id of the variable */
  int numCtr; /* Number of components which depend on this variable. */
  int varType; /* Domain, range or quantify variable ?*/
  lsList ctrItemList; /* List of ctrItemStruct corresponding to the clusters
                         that depend on it */
};



/**Struct**********************************************************************

  Synopsis    [This structure is the wrapper around a particular variable which
               is in the support of a particular cluster.]

  Description [Every cluster has a linked list of varItemStruct's. Each of
               these structs is a wrapper around varInfo. The "ctrItemStruct"
               corresponding to the cluster itself is present in ctrItemList,
               the linked list of varInfo. "ctrItemHandle" is the handle to
               the ctrItemStruct in "ctrItemList" of varInfo.]

******************************************************************************/
struct VarItemStruct {
  VarInfo_t *varInfo; /* The pointer to the corresponding variable info
                       * structure */
  lsHandle ctrItemHandle; /*
                           * The list handle to the ctrItemStruct
                           * corresponding to the ctrInfo (the one which
                           * contains this varItem in the varItemList) in the
                           * ctrItemList field of the varInfo.
                           */
};


/**Struct**********************************************************************

  Synopsis    [This structure is the wrapper around a particular ctr which
               depend on a particular variable.]

  Description [Every variable has a linked list of ctrItemStruct's. Each of
               these structs is a wrapper around ctrInfo. The "varItemStruct"
               corresponding to the variable itself is present in varItemList,
               the linked list of ctrInfo. "varItemHandle" is the handle to
               the varItemStruct in "varItemList" of ctrInfo.]

******************************************************************************/
struct CtrItemStruct {
  CtrInfo_t *ctrInfo;     /* Pointer to the corresponding cluster
                           * info structure. */
  lsHandle varItemHandle; /*
                           * The list handle to the varItemStruct
                           * corresponding to the varInfo (the one which
                           * contains this ctrItem in the ctrItemList) in the
                           * varItemList field of the ctrInfo.
                           */
};



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

static  double  *signatures;

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


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

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

static void ResetClusteredCofactoredRelationArray(mdd_manager *mddManager, Iwls95Info_t *info);
static void FreeClusteredCofactoredRelationArray(Iwls95Info_t *info);
static ImgTrmOption_t * TrmGetOptions(void);
static void CreateBitRelationArray(Iwls95Info_t *info, ImgFunctionData_t *functionData);
static void FreeBitRelationArray(Iwls95Info_t *info);
static void OrderClusterOrder(mdd_manager *mddManager, array_t *relationArray, array_t *fromVarBddArray, array_t *toVarBddArray, array_t *quantifyVarBddArray, array_t **orderedClusteredRelationArrayPtr, array_t **smoothVarBddArrayPtr, ImgTrmOption_t *option, boolean freeRelationArray);
static array_t * CreateClusters(bdd_manager *bddManager, array_t *relationArray, ImgTrmOption_t *option);
static void OrderRelationArray(mdd_manager *mddManager, array_t *relationArray, array_t *domainVarBddArray, array_t *quantifyVarBddArray, array_t *rangeVarBddArray, ImgTrmOption_t *option, array_t **orderedRelationArrayPtr, array_t **arraySmoothVarBddArrayPtr);
static array_t * RelationArraySmoothLocalVars(array_t *relationArray, array_t *ctrInfoArray, array_t *varInfoArray, st_table *bddIdToBddTable);
static void OrderRelationArrayAux(array_t *relationArray, lsList remainingCtrInfoList, array_t *ctrInfoArray, array_t *varInfoArray, int *sortedMaxIndexVector, int numSmoothVarsRemaining, int numIntroducedVarsRemaining, st_table *bddIdToBddTable, ImgTrmOption_t *option, array_t *domainAndQuantifyVarBddArray, array_t **orderedRelationArrayPtr, array_t **arraySmoothVarBddArrayPtr, array_t *arrayDomainQuantifyVarsWithZeroNumCtr);
static array_t * UpdateInfoArrays(CtrInfo_t *ctrInfo, st_table *bddIdToBddTable, int *numSmoothVarsRemainingPtr, int *numIntroducedVarsRemainingPtr);
static float CalculateBenefit(CtrInfo_t *ctrInfo, int maxNumLocalSmoothVars, int maxNumSmoothVars, int maxIndex, int maxNumIntroducedVars, ImgTrmOption_t *option);
static void PrintOption(Img_MethodType method, ImgTrmOption_t *option, FILE *fp);
static Iwls95Info_t * Iwls95InfoStructAlloc(Img_MethodType method);
static CtrInfo_t * CtrInfoStructAlloc(void);
static void CtrInfoStructFree(CtrInfo_t *ctrInfo);
static VarInfo_t * VarInfoStructAlloc(void);
static void VarInfoStructFree(VarInfo_t *varInfo);
static void CtrItemStructFree(CtrItem_t *ctrItem);
static void VarItemStructFree(VarItem_t *varItem);
static int CtrInfoMaxIndexCompare(const void * p1, const void * p2);
static void PrintCtrInfoStruct(CtrInfo_t *ctrInfo);
static void PrintVarInfoStruct(VarInfo_t *varInfo);
static int CheckQuantificationSchedule(array_t *relationArray, array_t *arraySmoothVarBddArray);
static int CheckCtrInfoArray(array_t *ctrInfoArray, int numDomainVars, int numQuantifyVars, int numRangeVars);
static int CheckCtrInfo(CtrInfo_t *ctrInfo, int numDomainVars, int numQuantifyVars, int numRangeVars);
static int CheckVarInfoArray(array_t *varInfoArray, int numRelation);
static array_t * BddArrayDup(array_t *bddArray);
static array_t * BddArrayArrayDup(array_t *bddArray);
static mdd_t * BddLinearAndSmooth(mdd_manager *mddManager, bdd_t *fromSet, array_t *relationArray, array_t *arraySmoothVarBddArray, array_t *smoothVarCubeArray, int verbosity, ImgPartialImageOption_t *PIoption, boolean *partial, boolean lazySiftFlag);
static void HashIdToBddTable(st_table *table, array_t *idArray, array_t *bddArray);
static void PrintSmoothIntroducedCount(array_t *clusterArray, array_t **arraySmoothVarBddArrayPtr, array_t *psBddIdArray, array_t *nsBddIdArray);
static void PrintBddIdFromBddArray(array_t *bddArray);
static void PrintBddIdTable(st_table *idTable);
static void PartitionTraverseRecursively(mdd_manager *mddManager, vertex_t *vertex, int mddId, st_table *vertexTable, array_t *relationArray, st_table *domainQuantifyVarMddIdTable, array_t *quantifyVarMddIdArray);
static void PrintPartitionRecursively(vertex_t *vertex, st_table *vertexTable, int indent);
static bdd_t * RecomputeImageIfNecessary(ImgFunctionData_t *functionData, mdd_manager *mddManager, bdd_t *domainSubset, array_t *relationArray, array_t *arraySmoothVarBddArray, array_t *smoothVarCubeArray, int verbosity, ImgPartialImageOption_t *PIoption, array_t *toCareSetArray, boolean *partial, boolean lazySiftFlag);
static bdd_t * ComputeSubsetOfIntermediateProduct(bdd_t *product, ImgPartialImageOption_t *PIoption);
static bdd_t * ComputeClippedAndAbstract(bdd_t *product, bdd_t *relation, array_t *smoothVarBddArray, int nvars, int clippingDepth, boolean *partial, int verbosity);
static array_t * CopyArrayBddArray(array_t *arrayBddArray);
static void PrintPartitionedTransitionRelation(mdd_manager *mddManager, Iwls95Info_t *info, Img_DirectionType directionType);
static void ReorderPartitionedTransitionRelation(Iwls95Info_t *info, ImgFunctionData_t *functionData, Img_DirectionType directionType);
static void UpdateQuantificationSchedule(Iwls95Info_t *info, ImgFunctionData_t *functionData, Img_DirectionType directionType);
static array_t * MakeSmoothVarCubeArray(mdd_manager *mddManager, array_t *arraySmoothVarBddArray);
static mdd_t * TrmEliminateDependVars(Img_ImageInfo_t *imageInfo, array_t *relationArray, mdd_t *states, mdd_t **dependRelations);
static int TrmSignatureCompare(int *ptrX, int *ptrY);
static void SetupLazySifting(mdd_manager *mddManager, array_t *bddRelationArray, array_t *domainVarBddArray, array_t *quantifyVarBddArray, array_t *rangeVarBddArray, int verbosity);
static int MddSizeCompare(const void *ptr1, const void *ptr2);
/**AutomaticEnd***************************************************************/


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


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

  Synopsis    [Returns the result after existentially quantifying a
  set of variables after taking the product of the array of relations.]

  Description ["relationArray" is an array of mdd's which need to be
  multiplied and the variables in the "smoothVarMddIdArray" need to
  be quantified out from the product. "introducedVarMddIdArray" is the
  array of mddIds of the variables (other than the variables to be
  quantified out) in the support of the relations. This array is used
  to compute the product order such that the number of new variables
  introduced in the product is minimized. However passing an empty
  array or an array of mddIds of partial support will not result in any
  error (some optimality will be lost though). The computation consists of 2
  phases. In phase 1, an ordering of the relations and a schedule of
  quantifying variables is found (based on IWLS95) heuristic. In phase
  2, the relations are multiplied in order and the quantifying
  variables are quantified according to the schedule.]

  SideEffects [None]

******************************************************************************/
mdd_t*
Img_MultiwayLinearAndSmooth(mdd_manager *mddManager, array_t *relationArray,
                            array_t *smoothVarMddIdArray,
                            array_t *introducedVarMddIdArray,
                            Img_MethodType method,
                            Img_DirectionType direction)
{
  mdd_t *product;
  array_t *smoothVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                      smoothVarMddIdArray);
  array_t *introducedVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                   introducedVarMddIdArray);

  product = ImgMultiwayLinearAndSmooth(mddManager, relationArray,
                        smoothVarBddArray, introducedVarBddArray,
                        method, direction, NIL(ImgTrmOption_t));

  mdd_array_free(introducedVarBddArray);
  mdd_array_free(smoothVarBddArray);
  return product;
}


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

  Synopsis    [Prints info of the partitioned transition relation.]

  Description [Prints info of the partitioned transition relation.]

  SideEffects []

******************************************************************************/
void
Img_PrintPartitionedTransitionRelation(mdd_manager *mddManager,
                                        Img_ImageInfo_t *imageInfo,
                                        Img_DirectionType directionType)
{
  Iwls95Info_t *info = (Iwls95Info_t *) imageInfo->methodData;

  if (imageInfo->methodType != Img_Iwls95_c &&
      imageInfo->methodType != Img_Mlp_c &&
      imageInfo->methodType != Img_Linear_c)
    return;

  PrintPartitionedTransitionRelation(mddManager, info, directionType);
}


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

  Synopsis    [Reorder partitioned transition relation.]

  Description [Reorder partitioned transition relation.]

  SideEffects []

******************************************************************************/
void
Img_ReorderPartitionedTransitionRelation(Img_ImageInfo_t *imageInfo,
                                         Img_DirectionType directionType)
{
  if (imageInfo->methodType != Img_Iwls95_c &&
      imageInfo->methodType != Img_Mlp_c &&
      imageInfo->methodType != Img_Linear_c)
    return;

  ReorderPartitionedTransitionRelation((Iwls95Info_t *)imageInfo->methodData,
                                        &(imageInfo->functionData),
                                        directionType);
}


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

  Synopsis    [Updates quantification schedule.]

  Description [Updates quantification schedule.]

  SideEffects []

******************************************************************************/
void
Img_UpdateQuantificationSchedule(Img_ImageInfo_t *imageInfo,
                                 Img_DirectionType directionType)
{
  if (imageInfo->methodType != Img_Iwls95_c)
    return;

  UpdateQuantificationSchedule((Iwls95Info_t *)imageInfo->methodData,
                                &(imageInfo->functionData), directionType);
}


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

  Synopsis    [Clusters relation array and returns quantification schedule.]

  Description [Clusters relation array and returns quantification schedule.]

  SideEffects []

******************************************************************************/
void
Img_ClusterRelationArray(mdd_manager *mddManager,
                        Img_MethodType method,
                        Img_DirectionType direction,
                        array_t *relationArray,
                        array_t *domainVarMddIdArray,
                        array_t *rangeVarMddIdArray,
                        array_t *quantifyVarMddIdArray,
                        array_t **clusteredRelationArray,
                        array_t **arraySmoothVarBddArray,
                        array_t **smoothVarCubeArray,
                        boolean freeRelationArray)
{
  array_t *domainVarBddArray;
  array_t *quantifyVarBddArray;
  array_t *rangeVarBddArray;
  ImgTrmOption_t *option;

  option = TrmGetOptions();
  domainVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                domainVarMddIdArray);
  rangeVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                rangeVarMddIdArray);
  quantifyVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                  quantifyVarMddIdArray);

  *clusteredRelationArray = ImgClusterRelationArray(mddManager,
                                                    NIL(ImgFunctionData_t),
                                                    method,
                                                    direction,
                                                    option,
                                                    relationArray,
                                                    domainVarBddArray,
                                                    quantifyVarBddArray,
                                                    rangeVarBddArray,
                                                    arraySmoothVarBddArray,
                                                    smoothVarCubeArray,
                                                    freeRelationArray);
  
  ImgFreeTrmOptions(option);
  mdd_array_free(domainVarBddArray);
  mdd_array_free(rangeVarBddArray);
  mdd_array_free(quantifyVarBddArray);
}


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


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

  Synopsis    [Returns the result after existentially quantifying a
  set of variables after taking the product of the array of relations.]

  Description ["relationArray" is an array of mdd's which need to be
  multiplied and the variables in the "smoothVarBddArray" need to
  be quantified out from the product. "introducedVarBddArray" is the
  array of the bdd variables (other than the variables to be
  quantified out) in the support of the relations. This array is used
  to compute the product order such that the number of new variables
  introduced in the product is minimized. However passing an empty
  array or an array of mddIds of partial support will not result in any
  error (some optimality will be lost though). The computation consists of 2
  phases. In phase 1, an ordering of the relations and a schedule of
  quantifying variables is found (based on IWLS95) heuristic. In phase
  2, the relations are multiplied in order and the quantifying
  variables are quantified according to the schedule.]

  SideEffects [None]

******************************************************************************/
mdd_t*
ImgMultiwayLinearAndSmooth(mdd_manager *mddManager,
                           array_t *relationArray,
                           array_t *smoothVarBddArray,
                           array_t *introducedVarBddArray,
                           Img_MethodType method,
                           Img_DirectionType direction,
                           ImgTrmOption_t *option)
{
  mdd_t *product;
  array_t *domainVarBddArray = array_alloc(bdd_t*, 0);

  product = ImgMultiwayLinearAndSmoothWithFrom(mddManager, relationArray,
                                                NIL(mdd_t),
                                                smoothVarBddArray,
                                                domainVarBddArray,
                                                introducedVarBddArray,
                                                method, direction, option);

  array_free(domainVarBddArray);
  return product;
}


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

  Synopsis    [Returns the result after existentially quantifying a
  set of variables after taking the product of the array of relations.]

  Description ["relationArray" is an array of mdd's which need to be
  multiplied and the variables in the "smoothVarBddArray" need to
  be quantified out from the product. "introducedVarBddArray" is the
  array of the bdd variables (other than the variables to be
  quantified out) in the support of the relations. This array is used
  to compute the product order such that the number of new variables
  introduced in the product is minimized. However passing an empty
  array or an array of mddIds of partial support will not result in any
  error (some optimality will be lost though). The computation consists of 2
  phases. In phase 1, an ordering of the relations and a schedule of
  quantifying variables is found (based on IWLS95) heuristic. In phase
  2, the relations are multiplied in order and the quantifying
  variables are quantified according to the schedule. If "from" is not nil,
  we first order "relationArray" and make quantification schedule, then
  we multiply with "from". In this case, "domainVarBddArray" should contain
  the state variables to be quantified out. If "from" is nil, then
  "domainVarBddArray" should be an empty array.]

  SideEffects [None]

******************************************************************************/
mdd_t*
ImgMultiwayLinearAndSmoothWithFrom(mdd_manager *mddManager,
                                   array_t *relationArray,
                                   mdd_t *from,
                                   array_t *smoothVarBddArray,
                                   array_t *domainVarBddArray,
                                   array_t *introducedVarBddArray,
                                   Img_MethodType method,
                                   Img_DirectionType direction,
                                   ImgTrmOption_t *option)
{
  mdd_t *product;
  boolean partial = FALSE;
  array_t *orderedRelationArray = NIL(array_t);
  array_t *arraySmoothVarBddArray = NIL(array_t);
  int freeOptionFlag;
  int lazySiftFlag = bdd_is_lazy_sift(mddManager);

  if (!option) {
    option = TrmGetOptions();
    freeOptionFlag = 1;
  } else
    freeOptionFlag = 0;

  if (method == Img_Iwls95_c) {
    OrderRelationArray(mddManager, relationArray, domainVarBddArray,
                        smoothVarBddArray, introducedVarBddArray, option,
                        &orderedRelationArray, &arraySmoothVarBddArray);
  } else if (method == Img_Mlp_c) {
    if (direction == Img_Forward_c) {
      /*
      * domainVarBddArray : PS
      * introducedVarBddArray : NS
      */
      if (option->mlpMultiway && from == NIL(mdd_t)) {
        product = ImgMlpMultiwayAndSmooth(mddManager, NIL(ImgFunctionData_t),
                                          relationArray, domainVarBddArray,
                                          smoothVarBddArray,
                                          introducedVarBddArray, direction,
                                          option);
        array_free(domainVarBddArray);
        if (freeOptionFlag)
          ImgFreeTrmOptions(option);
        return product;
      } else {
        ImgMlpOrderRelationArray(mddManager, relationArray, domainVarBddArray,
                                 smoothVarBddArray, introducedVarBddArray,
                                 direction, &orderedRelationArray,
                                 &arraySmoothVarBddArray, option);
      }
    } else { /* direction == Img_Backward_c */
      /*
      * domainVarBddArray : NS
      * introducedVarBddArray : PS
      */
      if (option->mlpPreSwapStateVars) {
        if (option->mlpMultiway && from == NIL(mdd_t)) {
          product = ImgMlpMultiwayAndSmooth(mddManager, NIL(ImgFunctionData_t),
                                            relationArray,
                                            domainVarBddArray,
                                            smoothVarBddArray,
                                            introducedVarBddArray,
                                            Img_Forward_c, option);
          array_free(domainVarBddArray);
          if (freeOptionFlag)
            ImgFreeTrmOptions(option);
          return product;
        } else {
          ImgMlpOrderRelationArray(mddManager, relationArray,
                                   domainVarBddArray, smoothVarBddArray,
                                   introducedVarBddArray,
                                   Img_Forward_c, &orderedRelationArray,
                                   &arraySmoothVarBddArray, option);
        }
      } else {
        if (option->mlpMultiway && from == NIL(mdd_t)) {
          product = ImgMlpMultiwayAndSmooth(mddManager, NIL(ImgFunctionData_t),
                                            relationArray, domainVarBddArray,
                                            smoothVarBddArray,
                                            introducedVarBddArray,
                                            Img_Backward_c, option);
          array_free(domainVarBddArray);
          if (freeOptionFlag)
            ImgFreeTrmOptions(option);
          return product;
        } else {
          ImgMlpOrderRelationArray(mddManager, relationArray,
                                   domainVarBddArray, smoothVarBddArray,
                                   introducedVarBddArray,
                                   Img_Backward_c, &orderedRelationArray,
                                   &arraySmoothVarBddArray, option);
        }
      }
    }
  } else
    assert(0);

  product = BddLinearAndSmooth(mddManager, from,
                               orderedRelationArray,
                               arraySmoothVarBddArray,
                               NIL(array_t),
                               option->verbosity,
                               NULL, &partial, lazySiftFlag);

  if (freeOptionFlag) 
    ImgFreeTrmOptions(option);
  mdd_array_free(orderedRelationArray);
  mdd_array_array_free(arraySmoothVarBddArray);
  return product;
}


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

  Synopsis [Returns a BDD after taking the product of fromSet with the
  BDDs in the relationArray and quantifying out the variables using
  the schedule given in the arraySmoothVarBddArray.]

  Description [The product is taken from the left, i.e., fromSet is
  multiplied with relationArray[0]. The array of variables given by
  arraySmoothVarBddArray[i] are quantified when the relationArray[i]
  is multiplied in the product.]

  SideEffects [None]

******************************************************************************/
mdd_t*
ImgBddLinearAndSmooth(mdd_manager *mddManager,
                        mdd_t *from,
                        array_t *relationArray,
                        array_t *arraySmoothVarBddArray,
                        array_t *smoothVarCubeArray,
                        int verbosity)
{
  mdd_t *product;
  boolean partial = FALSE;
  int lazySiftFlag = bdd_is_lazy_sift(mddManager);

  product = BddLinearAndSmooth(mddManager, from,
                               relationArray,
                               arraySmoothVarBddArray,
                               smoothVarCubeArray,
                               verbosity,
                               NULL, &partial, lazySiftFlag);

  return product;
}


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

  Synopsis    [Clusters relation array and makes quantification schedule.]

  Description [Orders and clusters relation array and makes quantification
  schedule. (Order/Cluster/Order.)]

  SideEffects []

******************************************************************************/
array_t *
ImgClusterRelationArray(mdd_manager *mddManager,
                        ImgFunctionData_t *functionData,
                        Img_MethodType method,
                        Img_DirectionType direction,
                        ImgTrmOption_t *option,
                        array_t *relationArray,
                        array_t *domainVarBddArray,
                        array_t *quantifyVarBddArray,
                        array_t *rangeVarBddArray,
                        array_t **arraySmoothVarBddArray,
                        array_t **smoothVarCubeArray,
                        boolean freeRelationArray)
{
  array_t       *orderedRelationArray;
  array_t       *clusteredRelationArray;
  int           freeOptionFlag;
  long          initialTime, finalTime;

  if (!option) {
    option = TrmGetOptions();
    freeOptionFlag = 1;
  } else
    freeOptionFlag = 0;

  if (option->verbosity >= 2)
    initialTime = util_cpu_time();
  else /* to remove uninitialized variables warning */
    initialTime = 0;

  if (method == Img_Iwls95_c) {
    assert(direction != Img_Both_c);
    if (direction == Img_Forward_c) {
      OrderClusterOrder(mddManager, relationArray, domainVarBddArray,
                        rangeVarBddArray, quantifyVarBddArray,
                        &orderedRelationArray, arraySmoothVarBddArray,
                        option, freeRelationArray);
    } else {
      assert(direction == Img_Backward_c);
      OrderClusterOrder(mddManager, relationArray, rangeVarBddArray,
                        domainVarBddArray, quantifyVarBddArray,
                        &orderedRelationArray, arraySmoothVarBddArray,
                        option, freeRelationArray);

    }
    if (freeOptionFlag)
      ImgFreeTrmOptions(option);

    if (smoothVarCubeArray) {
      if (functionData->createVarCubesFlag) {
        *smoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                                   *arraySmoothVarBddArray);
      } else
        *smoothVarCubeArray = NIL(array_t);
    }
    if (option->verbosity >= 2) {
      finalTime = util_cpu_time();
      fprintf(vis_stdout, "time for clustering with IWLS95 = %10g\n",
              (double)(finalTime - initialTime) / 1000.0);
    }
    return(orderedRelationArray);
  } else if (method == Img_Mlp_c) {
    if ((direction == Img_Forward_c && option->mlpReorder == 2) ||
        (direction == Img_Backward_c && option->mlpPreReorder == 2)) {
      if (direction == Img_Forward_c)
        option->mlpReorder = 0;
      else
        option->mlpPreReorder = 0;
      ImgMlpClusterRelationArray(mddManager, functionData, relationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                direction, &clusteredRelationArray, NIL(array_t *), option);
      if (direction == Img_Forward_c)
        option->mlpReorder = 2;
      else
        option->mlpPreReorder = 2;
      ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                direction, &orderedRelationArray, arraySmoothVarBddArray,
                option);
      mdd_array_free(clusteredRelationArray);
    } else if ((direction == Img_Forward_c && option->mlpReorder == 3) ||
               (direction == Img_Backward_c && option->mlpPreReorder == 3)) {
      if (direction == Img_Forward_c)
        option->mlpReorder = 0;
      else
        option->mlpPreReorder = 0;
      ImgMlpClusterRelationArray(mddManager, functionData, relationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                direction, &clusteredRelationArray, NIL(array_t *), option);
      if (direction == Img_Forward_c)
        option->mlpReorder = 3;
      else
        option->mlpPreReorder = 3;
      OrderRelationArray(mddManager, clusteredRelationArray,
                domainVarBddArray, quantifyVarBddArray,
                rangeVarBddArray, option,
                &orderedRelationArray, arraySmoothVarBddArray);
      mdd_array_free(clusteredRelationArray);
    } else {
      ImgMlpClusterRelationArray(mddManager, functionData, relationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                direction, &orderedRelationArray, arraySmoothVarBddArray,
                option);
    }
    if (freeOptionFlag)
      ImgFreeTrmOptions(option);
    if (smoothVarCubeArray) {
      if (functionData->createVarCubesFlag) {
        *smoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                                   *arraySmoothVarBddArray);
      } else
        *smoothVarCubeArray = NIL(array_t);
    }
    if (option->verbosity >= 2) {
      finalTime = util_cpu_time();
      fprintf(vis_stdout, "time for clustering with MLP = %10g\n",
                (double)(finalTime - initialTime) / 1000.0);
    }
    return(orderedRelationArray);
  } else
    assert(0);
  return NIL(array_t); /* we should never get here */
}


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

  Synopsis    [Returns quantification schedule for a given relation array.]

  Description [Returns quantification schedule for a given relation array.]

  SideEffects []

******************************************************************************/
array_t *
ImgGetQuantificationSchedule(mdd_manager *mddManager,
                                ImgFunctionData_t *functionData,
                                Img_MethodType method,
                                Img_DirectionType direction,
                                ImgTrmOption_t *option,
                                array_t *relationArray,
                                array_t *smoothVarBddArray,
                                array_t *domainVarBddArray,
                                array_t *introducedVarBddArray,
                                boolean withClustering,
                                array_t **orderedRelationArrayPtr)
{
  array_t *orderedRelationArray = NIL(array_t);
  array_t *clusteredRelationArray;
  array_t *arraySmoothVarBddArray = NIL(array_t);
  int freeOptionFlag;

  if (!option) {
    option = TrmGetOptions();
    freeOptionFlag = 1;
  } else
    freeOptionFlag = 0;

  if (method == Img_Iwls95_c) {
    if (withClustering) {
      OrderRelationArray(mddManager, relationArray, domainVarBddArray,
                        smoothVarBddArray, introducedVarBddArray, option,
                        &orderedRelationArray, NIL(array_t *));

      clusteredRelationArray = CreateClusters(mddManager, orderedRelationArray,
                                                option);
      mdd_array_free(orderedRelationArray);

      OrderRelationArray(mddManager, clusteredRelationArray, domainVarBddArray,
                        smoothVarBddArray, introducedVarBddArray, option,
                        &orderedRelationArray, &arraySmoothVarBddArray);
      mdd_array_free(clusteredRelationArray);
    } else {
      OrderRelationArray(mddManager, relationArray, domainVarBddArray,
                        smoothVarBddArray, introducedVarBddArray, option,
                        &orderedRelationArray, &arraySmoothVarBddArray);
    }
  } else if (method == Img_Mlp_c) {
    if (withClustering) {
      ImgMlpClusterRelationArray(mddManager, functionData, relationArray,
                domainVarBddArray, smoothVarBddArray, introducedVarBddArray,
                direction, &orderedRelationArray, &arraySmoothVarBddArray,
                option);
    } else {
      ImgMlpOrderRelationArray(mddManager, relationArray, domainVarBddArray,
                        smoothVarBddArray, introducedVarBddArray, direction,
                        &orderedRelationArray, &arraySmoothVarBddArray,
                        option);
    }
  } else
    assert(0);

  if (freeOptionFlag)
    ImgFreeTrmOptions(option);
  if (orderedRelationArrayPtr)
    *orderedRelationArrayPtr = orderedRelationArray;
  else
    mdd_array_free(orderedRelationArray);
  return(arraySmoothVarBddArray);
}


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

  Synopsis    [Allocates an option structure for transition relation method.]

  Description [Allocates an option structure for transition relation method.]

  SideEffects []

******************************************************************************/
ImgTrmOption_t *
ImgGetTrmOptions(void)
{
  ImgTrmOption_t        *option;

  option = TrmGetOptions();
  return(option);
}


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

  Synopsis    [Frees the option structure for transition relation method.]

  Description [Frees the option structure for transition relation method.]

  SideEffects []

******************************************************************************/
void
ImgFreeTrmOptions(ImgTrmOption_t *option)
{
  if (option->readCluster)
    FREE(option->readCluster);
  if (option->writeCluster)
    FREE(option->writeCluster);
  if (option->writeDepMatrix)
    FREE(option->writeDepMatrix);
  FREE(option);
}


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

  Synopsis    [Turn off option to apply linear based image computation by comsidering range of state variables]

  Description [ Turn off option to apply linear based image computation by comsidering range of state variables.]

  SideEffects []

  SeeAlso [ImgImageInfoFreeIwls95]
******************************************************************************/

void 
ImgResetMethodDataLinearComputeRange(void *methodData)
{
  Iwls95Info_t *info = (Iwls95Info_t *) methodData;
  info->linearComputeRange = 0;
}

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

  Synopsis    [Turn on option to apply linear based image computation by comsidering range of state variables]

  Description [ Turn on option to apply linear based image computation by comsidering range of state variables.]

  SideEffects []

  SeeAlso []
******************************************************************************/
void 
ImgSetMethodDataLinearComputeRange(void *methodData)
{
  Iwls95Info_t *info = (Iwls95Info_t *) methodData;
  info->linearComputeRange = 1;
}

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

  Synopsis    [Get unreachable states]

  Description [Return unreachable states]

  SideEffects []

  SeeAlso []
******************************************************************************/
mdd_t *
Img_ImageGetUnreachableStates(Img_ImageInfo_t * imageInfo)
{
  array_t *relationArray;
  mdd_t *states;
  Iwls95Info_t *info = (Iwls95Info_t *) imageInfo->methodData;

  relationArray = info->fwdOptClusteredRelationArray;
  states = array_fetch(mdd_t *, relationArray, 0);
  states = ImgSubstitute(states, &(imageInfo->functionData), Img_R2D_c);
  return(states);
}


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

  Synopsis    [Initializes an image data structure for image
  computation using the Iwls95 technique.]

  Description [This process consists of following steps.
  1. The transition functions are built.
  2. An array of bit level relations are built using the function
  mdd_fn_array_to_bdd_rel_array.
  3. The relations are clustered using the threshold value.
  4. The clustered relations are ordered.]

  SideEffects []

  SeeAlso [ImgImageInfoFreeIwls95]
******************************************************************************/


void *
ImgImageInfoInitializeIwls95(void *methodData,
                             ImgFunctionData_t *functionData,
                             Img_DirectionType directionType,
                             Img_MethodType method)
{
  array_t *fwdClusteredRelationArray;
  array_t *bwdClusteredRelationArray;
  Iwls95Info_t *info = (Iwls95Info_t *) methodData;
  Ntk_Network_t *network;
  Ntk_Node_t *node;
  char *flagValue;
  int i, mddId ;

  if (info &&
      (((info->fwdClusteredRelationArray) &&
        (info->bwdClusteredRelationArray)) ||
       ((directionType == Img_Forward_c) &&
        (info->fwdClusteredRelationArray)) ||
       ((directionType == Img_Backward_c) &&
        (info->bwdClusteredRelationArray)))) {
    /* Nothing needs to be done. Return */
    return (void *) info;
  }
  else{
    array_t *rangeVarBddArray, *domainVarBddArray, *quantifyVarBddArray;
    array_t *newQuantifyVarMddIdArray;
    array_t *clusteredRelationArray;
    FILE *fin, *fout;
    array_t *bddRelationArray, *quantifyVarMddIdArray;
    mdd_manager *mddManager = Part_PartitionReadMddManager(functionData->mddNetwork);

    if (info == NIL(Iwls95Info_t)) {
      info = Iwls95InfoStructAlloc(method);
    }
    CreateBitRelationArray(info, functionData);

    /** if FAFWFlag then remove PI from quantifyVarMddIdArray **/
    bddRelationArray = info->bitRelationArray;
    quantifyVarMddIdArray = info->quantifyVarMddIdArray;
    if(functionData->FAFWFlag) {
      network = functionData->network;
      newQuantifyVarMddIdArray = array_alloc(int, 0);
      for(i=0; i<array_n(quantifyVarMddIdArray); i++) {
        mddId = array_fetch(int, quantifyVarMddIdArray, i);
        node = Ntk_NetworkFindNodeByMddId(network, mddId);
        if(!Ntk_NodeTestIsInput(node)) {
          array_insert_last(int, newQuantifyVarMddIdArray, mddId);
        }
      }
      array_free(quantifyVarMddIdArray);
      quantifyVarMddIdArray = newQuantifyVarMddIdArray;
      info->quantifyVarMddIdArray = newQuantifyVarMddIdArray;
    }

    info->bitRelationArray = NIL(array_t);
    info->quantifyVarMddIdArray = NIL(array_t);

    rangeVarBddArray = functionData->rangeBddVars;
    domainVarBddArray = functionData->domainBddVars;
    quantifyVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                    quantifyVarMddIdArray);

    array_free(quantifyVarMddIdArray);

    if (((directionType == Img_Forward_c) ||
         (directionType == Img_Both_c)) &&
        (info->fwdClusteredRelationArray == NIL(array_t))) {
      if (info->option->readCluster) {
        if (info->option->readReorderCluster) {
          fin = fopen(info->option->readCluster, "r");
          ImgMlpReadClusterFile(fin, mddManager,
                          functionData, bddRelationArray,
                          domainVarBddArray, rangeVarBddArray,
                          quantifyVarBddArray, Img_Forward_c,
                          &clusteredRelationArray, NIL(array_t *),
                          info->option);
          fclose(fin);
          if (method == Img_Iwls95_c) {
            OrderRelationArray(mddManager, clusteredRelationArray,
                          domainVarBddArray, quantifyVarBddArray,
                          rangeVarBddArray, info->option,
                          &info->fwdClusteredRelationArray,
                          &info->fwdArraySmoothVarBddArray);
          } else if (method == Img_Mlp_c) {
            ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                          domainVarBddArray, quantifyVarBddArray,
                          rangeVarBddArray, Img_Forward_c,
                          &info->fwdClusteredRelationArray,
                          &info->fwdArraySmoothVarBddArray,
                          info->option);
          } else
            assert(0);
          mdd_array_free(clusteredRelationArray);
        } else {
          fin = fopen(info->option->readCluster, "r");
          ImgMlpReadClusterFile(fin, mddManager,
                          functionData, bddRelationArray,
                          domainVarBddArray, rangeVarBddArray,
                          quantifyVarBddArray, Img_Forward_c,
                          &info->fwdClusteredRelationArray,
                          &info->fwdArraySmoothVarBddArray,
                          info->option);
          fclose(fin);
        }
      } else if (method == Img_Iwls95_c) {
        /*
         * Since clusters are formed by multiplying the latches starting
         * from one end we need to order the latches using some heuristics
         * first. Right now, we order the latches using the same heuristic
         * as the one used for ordering the clusters for early quantifiction.
         * However, since we are not interested in the quantification
         * schedule at this stage, we will pass a NIL(array_t*) as the last
         * argument.
         */
        OrderClusterOrder(mddManager, bddRelationArray, domainVarBddArray,
                          rangeVarBddArray, quantifyVarBddArray,
                          &info->fwdClusteredRelationArray,
                          &info->fwdArraySmoothVarBddArray,
                          info->option, 0);
      } else if (method == Img_Mlp_c) {
        if (info->option->mlpReorder == 2) {
          info->option->mlpReorder = 0;
          ImgMlpClusterRelationArray(mddManager, functionData,
                        bddRelationArray, domainVarBddArray,
                        quantifyVarBddArray, rangeVarBddArray,
                        Img_Forward_c, &clusteredRelationArray,
                        NIL(array_t *), info->option);
          info->option->mlpReorder = 2;
          ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, Img_Forward_c,
                        &info->fwdClusteredRelationArray,
                        &info->fwdArraySmoothVarBddArray, info->option);
          mdd_array_free(clusteredRelationArray);
        } else if (info->option->mlpReorder == 3) {
          info->option->mlpReorder = 0;
          ImgMlpClusterRelationArray(mddManager, functionData,
                        bddRelationArray, domainVarBddArray,
                        quantifyVarBddArray, rangeVarBddArray,
                        Img_Forward_c, &clusteredRelationArray,
                        NIL(array_t *), info->option);
          info->option->mlpReorder = 3;
          OrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, info->option,
                        &info->fwdClusteredRelationArray,
                        &info->fwdArraySmoothVarBddArray);
          mdd_array_free(clusteredRelationArray);
        } else {
          ImgMlpClusterRelationArray(mddManager, functionData,
                        bddRelationArray, domainVarBddArray,
                        quantifyVarBddArray, rangeVarBddArray,
                        Img_Forward_c, &info->fwdClusteredRelationArray,
                        &info->fwdArraySmoothVarBddArray, info->option);
        }
      } else if (method == Img_Linear_c || method == Img_Linear_Range_c) {
          if(method == Img_Linear_Range_c) {
            info->option->linearComputeRange = 1;
            method = Img_Linear_c;
          }
          ImgLinearClusterRelationArray(mddManager, functionData,
                        bddRelationArray, 
                        Img_Forward_c, 
                        &info->fwdClusteredRelationArray,
                        &info->fwdArraySmoothVarBddArray, 
                        &info->fwdOptClusteredRelationArray,
                        &info->fwdOptArraySmoothVarBddArray, 
                        info->option);
      } else
        assert(0);

      if (info->option->verbosity > 2) 
        PrintPartitionedTransitionRelation(mddManager, info, Img_Forward_c);
      if(method == Img_Linear_c)
        ImgPrintPartitionedTransitionRelation(mddManager,
                info->fwdOptClusteredRelationArray , 
                info->fwdOptArraySmoothVarBddArray);

      if (info->option->printDepMatrix || info->option->writeDepMatrix) {
        if (info->option->writeDepMatrix) {
          fout = fopen(info->option->writeDepMatrix, "w");
          if (!fout) {
            fprintf(vis_stderr, "** img error: can't open file [%s]\n",
                    info->option->writeDepMatrix);
          }
        } else
          fout = NIL(FILE);
        ImgMlpPrintDependenceMatrix(mddManager,
                info->fwdClusteredRelationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                Img_Forward_c, info->option->printDepMatrix, fout,
                info->option);
        if (fout)
          fclose(fout);
      }
      if (info->option->writeCluster) {
        fout = fopen(info->option->writeCluster, "w");
        ImgMlpWriteClusterFile(fout, mddManager, 
                info->fwdClusteredRelationArray,
                domainVarBddArray, rangeVarBddArray);
        fclose(fout);
      }

      if (functionData->createVarCubesFlag && info->fwdArraySmoothVarBddArray) 
        info->fwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->fwdArraySmoothVarBddArray);
      else
        info->fwdSmoothVarCubeArray = NIL(array_t);

      if (functionData->createVarCubesFlag && info->fwdOptArraySmoothVarBddArray) 
        info->fwdOptSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->fwdOptArraySmoothVarBddArray);
      else
        info->fwdOptSmoothVarCubeArray = NIL(array_t);

#ifndef NDEBUG
      if(info->fwdClusteredRelationArray)
        assert(CheckQuantificationSchedule(info->fwdClusteredRelationArray,
                                         info->fwdArraySmoothVarBddArray));
      if(info->fwdOptClusteredRelationArray)
        assert(CheckQuantificationSchedule(info->fwdOptClusteredRelationArray,
                                         info->fwdOptArraySmoothVarBddArray));
#endif
    }

    if (((directionType == Img_Backward_c) ||
         (directionType == Img_Both_c)) &&
        (info->bwdClusteredRelationArray == NIL(array_t))) {
      if (info->option->readCluster) {
        if (info->option->readReorderCluster) {
          fin = fopen(info->option->readCluster, "r");
          ImgMlpReadClusterFile(fin, mddManager,
                          functionData, bddRelationArray,
                          domainVarBddArray, rangeVarBddArray,
                          quantifyVarBddArray, Img_Backward_c,
                          &clusteredRelationArray, NIL(array_t *),
                          info->option);
          fclose(fin);
          if (method == Img_Iwls95_c) {
            OrderRelationArray(mddManager, clusteredRelationArray,
                          rangeVarBddArray, quantifyVarBddArray,
                          domainVarBddArray, info->option,
                          &info->bwdClusteredRelationArray,
                          &info->bwdArraySmoothVarBddArray);
          } else if (method == Img_Mlp_c) {
            ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                          domainVarBddArray, quantifyVarBddArray,
                          rangeVarBddArray, Img_Backward_c,
                          &info->bwdClusteredRelationArray,
                          &info->bwdArraySmoothVarBddArray,
                          info->option);
          } else
            assert(0);
          mdd_array_free(clusteredRelationArray);
        } else {
          fin = fopen(info->option->readCluster, "r");
          ImgMlpReadClusterFile(fin, mddManager,
                          functionData, bddRelationArray,
                          domainVarBddArray, rangeVarBddArray,
                          quantifyVarBddArray, Img_Backward_c,
                          &info->bwdClusteredRelationArray,
                          &info->bwdArraySmoothVarBddArray,
                          info->option);
          fclose(fin);
        }
      } else if (method == Img_Iwls95_c) {
        /*
         * Since clusters are formed by multiplying the latches starting
         * from one end we need to order the latches using some heuristics
         * first. Right now, we order the latches using the same heuristic
         * as the one used for ordering the clusters for early quantifiction.
         * However, since we are not interested in the quantification
         * schedule at this stage, we will pass a NIL(array_t*) as the last
         * argument.
         */
        OrderClusterOrder(mddManager, bddRelationArray, rangeVarBddArray,
                          domainVarBddArray, quantifyVarBddArray,
                          &info->bwdClusteredRelationArray,
                          &info->bwdArraySmoothVarBddArray,
                          info->option, 0);
      } else if (method == Img_Mlp_c) {
        if (info->option->mlpPreReorder == 2) {
          info->option->mlpPreReorder = 0;
          ImgMlpClusterRelationArray(mddManager, functionData,
                        bddRelationArray, domainVarBddArray,
                        quantifyVarBddArray, rangeVarBddArray,
                        Img_Backward_c, &clusteredRelationArray,
                        NIL(array_t *), info->option);
          info->option->mlpPreReorder = 2;
          ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, Img_Backward_c,
                        &info->bwdClusteredRelationArray,
                        &info->bwdArraySmoothVarBddArray, info->option);
          mdd_array_free(clusteredRelationArray);
        } else if (info->option->mlpPreReorder == 3) {
          info->option->mlpPreReorder = 0;
          ImgMlpClusterRelationArray(mddManager, functionData,
                        bddRelationArray, domainVarBddArray,
                        quantifyVarBddArray, rangeVarBddArray,
                        Img_Backward_c, &clusteredRelationArray,
                        NIL(array_t *), info->option);
          info->option->mlpPreReorder = 3;
          OrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, info->option,
                        &info->bwdClusteredRelationArray,
                        &info->bwdArraySmoothVarBddArray);
          mdd_array_free(clusteredRelationArray);
        } else {
          ImgMlpClusterRelationArray(mddManager, functionData,
                        bddRelationArray, domainVarBddArray,
                        quantifyVarBddArray, rangeVarBddArray,
                        Img_Backward_c, &info->bwdClusteredRelationArray,
                        &info->bwdArraySmoothVarBddArray, info->option);
        }
      } else
        assert(0);
      
      info->careSetArray = NIL(array_t);
      info->bwdClusteredCofactoredRelationArray = NIL(array_t);
      
      if (info->option->verbosity > 2) 
        PrintPartitionedTransitionRelation(mddManager, info, Img_Backward_c);

      if (info->option->printDepMatrix || info->option->writeDepMatrix) {
        if (info->option->writeDepMatrix) {
          fout = fopen(info->option->writeDepMatrix, "w");
          if (!fout) {
            fprintf(vis_stderr, "** img error: can't open file [%s]\n",
                    info->option->writeDepMatrix);
          }
        } else
          fout = NIL(FILE);
        ImgMlpPrintDependenceMatrix(mddManager,
                info->bwdClusteredRelationArray,
                domainVarBddArray, quantifyVarBddArray, rangeVarBddArray,
                Img_Backward_c, info->option->printDepMatrix, fout,
                info->option);
      
        if (fout)
          fclose(fout);
      }
      if (info->option->writeCluster) {
        fout = fopen(info->option->writeCluster, "w");
        ImgMlpWriteClusterFile(fout, mddManager, 
                info->bwdClusteredRelationArray,
                domainVarBddArray, rangeVarBddArray);
        fclose(fout);
      }

      if (functionData->createVarCubesFlag)
        info->bwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->bwdArraySmoothVarBddArray);
      else
        info->bwdSmoothVarCubeArray = NIL(array_t);

#ifndef NDEBUG
      assert(CheckQuantificationSchedule(info->bwdClusteredRelationArray,
                                         info->bwdArraySmoothVarBddArray));
#endif
    }

    
    info->lazySiftFlag = bdd_is_lazy_sift(mddManager);
    if (info->lazySiftFlag) {
      SetupLazySifting(mddManager, bddRelationArray, domainVarBddArray,
                       quantifyVarBddArray, rangeVarBddArray,
                       info->option->verbosity);
    }
    
    /* Free the bddRelationArray */
    mdd_array_free(bddRelationArray);
    
    if ((info->option)->verbosity > 0) {
      fwdClusteredRelationArray = info->fwdClusteredRelationArray;
      bwdClusteredRelationArray = info->bwdClusteredRelationArray;
      if (method == Img_Iwls95_c)
        fprintf(vis_stdout,"Computing Image Using IWLS95 technique.\n");
      else if (method == Img_Mlp_c)
        fprintf(vis_stdout,"Computing Image Using MLP technique.\n");
      else if (method == Img_Linear_c) {
        fprintf(vis_stdout,"Computing Image Using LINEAR technique.\n");
        fwdClusteredRelationArray = info->fwdOptClusteredRelationArray;
        /*bwdClusteredRelationArray = info->bwdOptClusteredRelationArray;*/
      }
      else
        assert(0);
      fprintf(vis_stdout,"Total # of domain binary variables = %3d\n",
              array_n(domainVarBddArray));
      fprintf(vis_stdout,"Total # of range binary variables = %3d\n",
              array_n(rangeVarBddArray));
      fprintf(vis_stdout,"Total # of quantify binary variables = %3d\n",
              array_n(quantifyVarBddArray));
      if ((directionType == Img_Forward_c) || (directionType ==
                                               Img_Both_c)) {
        fprintf(vis_stdout, "Shared size of transition relation for forward image computation is %10ld BDD nodes (%-4d components)\n",
                bdd_size_multiple(fwdClusteredRelationArray),
                array_n(fwdClusteredRelationArray));
      }
      if ((directionType == Img_Backward_c) || (directionType == Img_Both_c)) {
        fprintf(vis_stdout, "Shared size of transition relation for backward image computation is %10ld BDD nodes (%-4d components)\n",
                bdd_size_multiple(bwdClusteredRelationArray),
                array_n(bwdClusteredRelationArray));
      }
    }
    
    /* Update nvars field in partial image options structure with
       total number of variables. */
    if (info->PIoption != NULL) {
      (info->PIoption)->nvars = array_n(domainVarBddArray) +
        array_n(quantifyVarBddArray) +
        array_n(rangeVarBddArray);
    }
    /* Free the Bdd arrays */
    mdd_array_free(quantifyVarBddArray);
    
    flagValue = Cmd_FlagReadByName("image_eliminate_depend_vars");
    if (flagValue == NIL(char))
      info->eliminateDepend = 0; /* the default value */
    else
      info->eliminateDepend = atoi(flagValue);
    if (info->eliminateDepend && bdd_get_package_name() != CUDD) {
      fprintf(vis_stderr,
"** trm error : image_eliminate_depend_vars is available for only CUDD.\n");
      info->eliminateDepend = 0;
    }
    info->nPrevEliminatedFwd = -1;
    
    return (void *)info;
  }
}


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

  Synopsis    [Computes the forward image of a set of states.]

  Description [Computes the forward image of a set of states. If
  partial image requested and the first image computed doesn't
  intersect toCareSet, then recompute image. ASSUME that toCareSet is
  in terms of functionData->domainVars.]

  SideEffects [ IMPORTANT NOTE: image may be in terms of domain
  variables even when the reverse is expected IF partial images are
  allowed. The clean way would be to convert toCareset in terms of
  functionData->rangeVars but it may end up being quite inefficient.]

  SeeAlso     [ImgImageInfoComputeFwdWithDomainVarsIwls95]

******************************************************************************/
mdd_t *
ImgImageInfoComputeFwdIwls95(ImgFunctionData_t *functionData,
                             Img_ImageInfo_t *imageInfo,
                             mdd_t *fromLowerBound,
                             mdd_t *fromUpperBound,
                             array_t *toCareSetArray)
{
  mdd_t *image, *domainSubset, *newFrom;
  mdd_manager *mddManager;
  boolean partial;
  Iwls95Info_t *info = (Iwls95Info_t *)imageInfo->methodData;
  int eliminateDepend;
  Img_OptimizeType optDir;
  array_t *fwdOriClusteredRelationArray;
  array_t *fwdOriArraySmoothVarBddArray;
  array_t *fwdOriSmoothVarCubeArray;
  array_t *fwdClusteredRelationArray;
  Relation_t *head;
  int bOptimize;
  boolean farSideImageIsEnabled;
  char    *userSpecifiedMethod;

  mddManager = Part_PartitionReadMddManager(functionData->mddNetwork);

  if(imageInfo->methodType == Img_Linear_c) {
    if(info->fwdClusteredRelationArray == 0 || 
       imageInfo->useOptimizedRelationFlag) {
      fwdOriClusteredRelationArray = info->fwdOptClusteredRelationArray;
      fwdOriArraySmoothVarBddArray = info->fwdOptArraySmoothVarBddArray;
      fwdOriSmoothVarCubeArray = info->fwdOptSmoothVarCubeArray;
    }
    else {
      fwdOriClusteredRelationArray = info->fwdClusteredRelationArray;
      fwdOriArraySmoothVarBddArray = info->fwdArraySmoothVarBddArray;
      fwdOriSmoothVarCubeArray = info->fwdSmoothVarCubeArray;
    }
  }
  else {
    fwdOriClusteredRelationArray = info->fwdClusteredRelationArray;
    fwdOriArraySmoothVarBddArray = info->fwdArraySmoothVarBddArray;
    fwdOriSmoothVarCubeArray = info->fwdSmoothVarCubeArray;
  }

  if (fwdOriClusteredRelationArray == NIL(array_t)) {
    fprintf(vis_stderr, "** img error: The data structure has not been ");
    fprintf(vis_stderr, "initialized for forward image computation\n");
    return NIL(mdd_t);
  }

  if (info->eliminateDepend == 1 ||
      (info->eliminateDepend == 2 && info->nPrevEliminatedFwd > 0)) {
    eliminateDepend = 1;
  } else
    eliminateDepend = 0;

  domainSubset = bdd_between(fromLowerBound, fromUpperBound);

  if (eliminateDepend) {
    fwdClusteredRelationArray = mdd_array_duplicate(fwdOriClusteredRelationArray);
    newFrom = TrmEliminateDependVars(imageInfo, fwdClusteredRelationArray,
                                     domainSubset, NIL(mdd_t *));
    mdd_free(domainSubset);
    domainSubset = newFrom;
  } else
    fwdClusteredRelationArray = fwdOriClusteredRelationArray;

  /* Record if partial images are allowed. */
  partial = info->allowPartialImage;
  info->wasPartialImage = FALSE;

  /* bias field is used to bias the intermediate results towards care states */
  if ((info->PIoption)->partialImageApproxMethod == BDD_APPROX_BIASED_RUA) {
    /* assume that this array has only one mdd */
    mdd_t *toCareSet = array_fetch(mdd_t *, toCareSetArray, 0);
    (info->PIoption)->bias = mdd_dup(toCareSet);
  }

  /* when the "compositional far side image approach" is selected, and
   * there are more than one transition relation clusters,
   * over-approximate images (w.r.t individual TR clusters) will be
   * used to minimize the transition relation.
   */
  farSideImageIsEnabled = FALSE;
  userSpecifiedMethod = Cmd_FlagReadByName("image_farside_method");
  if (userSpecifiedMethod != NIL(char)) {
    if (!strcmp(userSpecifiedMethod, "1")) 
      farSideImageIsEnabled = TRUE;
  }

  if (!farSideImageIsEnabled || array_n(fwdClusteredRelationArray) < 2) {

    image = BddLinearAndSmooth(mddManager, domainSubset,
                               fwdClusteredRelationArray,
                               fwdOriArraySmoothVarBddArray,
                               fwdOriSmoothVarCubeArray,
                               (info->option)->verbosity,
                               info->PIoption, &partial, info->lazySiftFlag);

  }else {

    array_t  *fwdFarArraySmoothVarBddArray = fwdOriArraySmoothVarBddArray;
    array_t  *fwdFarSmoothVarCubeArray      = fwdOriSmoothVarCubeArray;
    array_t  *fwdFarClusteredRelationArray = mdd_array_duplicate(fwdClusteredRelationArray);
    array_t  *upperBoundImages = array_alloc(mdd_t *, 0);
    st_table *domainSupportVarTable, *trClusterSupportVarTable;
    mdd_t    *trCluster, *newDomainset, *newTRcluster, *upperBoundImage;
    mdd_t    *holdMdd, *tmpMdd;
    array_t  *tempArray;
    int      beforeSize, afterSize, middleSize;
    int      size1, size2, size3;
    boolean  flag1, flag2;
    int      i,j, k;
    long     bddId;

    /* compute a series of over-approximate images based on individual
     * TR clusters
     */
    domainSupportVarTable = st_init_table(st_ptrcmp, st_ptrhash);
    tempArray = mdd_get_bdd_support_ids(mddManager, domainSubset);
    arrayForEachItem(int, tempArray, j, bddId) {
      st_insert(domainSupportVarTable, (char *)bddId, (char *)0);
    }
    array_free(tempArray);

    beforeSize = afterSize = middleSize =0;
    arrayForEachItem(mdd_t *, fwdFarClusteredRelationArray, i, trCluster) {
      array_t *preQuantifiedVarsInDomainset = array_alloc(mdd_t *, 0);
      array_t *preQuantifiedVarsInTRcluster = array_alloc(mdd_t *, 0);
      array_t *quantifiedVarsInProduct = array_alloc(mdd_t *, 0);

      /* some variables can be quantified before taking the conjunction */
      trClusterSupportVarTable = st_init_table(st_ptrcmp, st_ptrhash);
      tempArray = mdd_get_bdd_support_ids(mddManager, trCluster);
      arrayForEachItem(int, tempArray, j, bddId) {
        st_insert(trClusterSupportVarTable, (char *)bddId, (char *)0);
      }
      array_free(tempArray);
      arrayForEachItem(array_t *, fwdFarArraySmoothVarBddArray, j, tempArray) {
        arrayForEachItem(mdd_t *, tempArray, k, holdMdd) {
          bddId = (int)bdd_top_var_id(holdMdd);
          flag1 = st_is_member(domainSupportVarTable, (char *)bddId);
          flag2 = st_is_member(trClusterSupportVarTable, (char *)bddId);
          if (flag1 && flag2)
            array_insert_last(mdd_t *, quantifiedVarsInProduct, holdMdd);
          else if (flag1)
            array_insert_last(mdd_t *, preQuantifiedVarsInDomainset, holdMdd);
          else if (flag2)
            array_insert_last(mdd_t *, preQuantifiedVarsInTRcluster, holdMdd);
          else 
            ;
        }
      }
      st_free_table(trClusterSupportVarTable);

      if (array_n(preQuantifiedVarsInDomainset) == 0)
        newDomainset = bdd_dup(domainSubset);
      else
        newDomainset = bdd_smooth(domainSubset, preQuantifiedVarsInDomainset);
      array_free(preQuantifiedVarsInDomainset);

      if (array_n(preQuantifiedVarsInTRcluster) == 0)
        newTRcluster = bdd_dup(trCluster);
      else
        newTRcluster = bdd_smooth(trCluster, preQuantifiedVarsInTRcluster);
      array_free(preQuantifiedVarsInTRcluster);

      /* compute the over-approximated image */
      if (array_n(quantifiedVarsInProduct) == 0)
        upperBoundImage = bdd_and(newTRcluster, newDomainset, 1, 1);
      else
        upperBoundImage = bdd_and_smooth(newTRcluster, newDomainset, quantifiedVarsInProduct);
      array_free(quantifiedVarsInProduct);
      mdd_free(newTRcluster);
      mdd_free(newDomainset);
      
      /* minimize the TR cluster with the new upper bound image */
      newTRcluster = Img_MinimizeImage(trCluster, upperBoundImage, 
                                       Img_DefaultMinimizeMethod_c, TRUE);

      size1 = mdd_size(trCluster);
      beforeSize += size1;
      size2 = mdd_size(upperBoundImage);
      middleSize += size2;
      size3 = mdd_size(newTRcluster);
      afterSize += size3;
      if (info->option->verbosity >= 3) {
        fprintf(vis_stdout, "farSideImg: minimize TR cluster[%d]  %d |%d => %d\n", i, 
                size1, size2, size3);
      }
      
      if (mdd_equal(newTRcluster, trCluster)) {
        mdd_free(newTRcluster);
        mdd_free(upperBoundImage);
      }else {
        mdd_free(trCluster);
        array_insert(mdd_t *, fwdFarClusteredRelationArray, i, newTRcluster);
        array_insert_last(mdd_t *, upperBoundImages, upperBoundImage);
      }
    }
    st_free_table(domainSupportVarTable);

    if (info->option->verbosity >= 2) {
      fprintf(vis_stdout, "farSideImg: minimize TR total  %d |%d => %d\n",
              beforeSize, middleSize, afterSize);
    }
    
    /* compute the image with the simplified transition relation */
    image = BddLinearAndSmooth(mddManager, domainSubset,
                               fwdFarClusteredRelationArray,
                               fwdFarArraySmoothVarBddArray,
                               fwdFarSmoothVarCubeArray,
                               (info->option)->verbosity,
                               info->PIoption, &partial, info->lazySiftFlag);
    mdd_array_free(fwdFarClusteredRelationArray);

    /* for performance concerns, the following operations are carried
     *  out in the domain space
     */
    arrayForEachItem(mdd_t *, upperBoundImages, i, holdMdd) { 
      tmpMdd = ImgSubstitute(holdMdd, functionData, Img_R2D_c); 
      mdd_free(holdMdd); 
      array_insert(mdd_t *, upperBoundImages, i, tmpMdd); 
    } 
    holdMdd = image; 
    image = ImgSubstitute(image, functionData,Img_R2D_c); 
    mdd_free(holdMdd); 
    
    /* remove the "bad states" introduced by the TR cluster minimization */
    beforeSize = mdd_size(image);
    array_sort(upperBoundImages, MddSizeCompare);    
    arrayForEachItem(mdd_t *, upperBoundImages, i, tmpMdd) {
      size1 = mdd_size(image);
      size2 = mdd_size(tmpMdd);
      holdMdd = image;
      image = mdd_and(image, tmpMdd, 1, 1);
      mdd_free(tmpMdd);
      mdd_free(holdMdd);
      size3 = mdd_size(image);

      if (info->option->verbosity >= 3) {
        fprintf(vis_stdout, "farSideImg: clip image  %d | %d => %d\n",
                size1, size2, size3);
      }
    }
    array_free(upperBoundImages);
    
    if (info->option->verbosity >= 2) {
      afterSize = mdd_size(image);
      fprintf(vis_stdout, "farSideImg: clip image total  %d => %d \n",
              beforeSize, afterSize);
    }

  }


  /* check that no partial image was computed if not allowed */
  assert(!((!info->allowPartialImage) && partial));

  /* if partial image computed, and if no new states, recompute
   * image with relaxed parameters.
   */
  if (partial) {
    if (bdd_leq_array(image, toCareSetArray, 1, 0)) {
      mdd_free(image);
      partial = info->allowPartialImage;
      image = RecomputeImageIfNecessary(functionData,
                                        mddManager,
                                        domainSubset,
                                        fwdClusteredRelationArray,
                                        fwdOriArraySmoothVarBddArray,
                                        fwdOriSmoothVarCubeArray,
                                        (info->option)->verbosity,
                                        info->PIoption,
                                        toCareSetArray,
                                        &partial, info->lazySiftFlag);
    }
  }

  if(imageInfo->methodType == Img_Linear_c) {
    if(imageInfo->useOptimizedRelationFlag) {
      if(info->option->linearOptimize == Opt_NS ||
         info->option->linearOptimize == Opt_Both)
        optDir = Opt_Both;
      else
        optDir = Opt_CS;
  
      head = ImgLinearRelationInit(mddManager, fwdOriClusteredRelationArray, 
            functionData->domainBddVars, functionData->rangeBddVars, 
            functionData->quantifyBddVars, info->option);
      bOptimize = ImgLinearOptimizeAll(head, optDir, 0);
      if(bOptimize) {
        ImgLinearRefineRelation(head);
        if(info->fwdClusteredRelationArray == 0 && optDir == Opt_Both) {
          info->fwdClusteredRelationArray = fwdOriClusteredRelationArray;
          info->fwdArraySmoothVarBddArray = 
                   BddArrayArrayDup(fwdOriArraySmoothVarBddArray);
          info->fwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->fwdArraySmoothVarBddArray);
        }
        else {
          mdd_array_free(fwdOriClusteredRelationArray);
        }
        fwdOriClusteredRelationArray = ImgLinearExtractRelationArrayT(head);
        info->fwdOptClusteredRelationArray = fwdOriClusteredRelationArray;
      }
      ImgLinearRelationQuit(head);
    }
  }

  /* free the bias function created. */
  if ((info->PIoption)->partialImageApproxMethod == BDD_APPROX_BIASED_RUA) {
    mdd_free((info->PIoption)->bias);
    (info->PIoption)->bias = NIL(mdd_t);
  }
  /* indicate whether this image computation was partial or exact. */
  info->wasPartialImage = partial;
  info->allowPartialImage = FALSE;

  if (eliminateDepend)
    mdd_array_free(fwdClusteredRelationArray);
  mdd_free(domainSubset);

  return image;
}

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

  Synopsis    [Computes the forward image of a set of states in terms
  of domain variables.]

  Description [First the forward image computation function is called
  which returns an image on range vars. Later, variable substitution
  is used to obtain image on domain vars.]

  SideEffects []

  SeeAlso     [ImgImageInfoComputeFwdIwls95]

******************************************************************************/
mdd_t *
ImgImageInfoComputeFwdWithDomainVarsIwls95(ImgFunctionData_t *functionData,
                                           Img_ImageInfo_t *imageInfo,
                                           mdd_t *fromLowerBound,
                                           mdd_t *fromUpperBound,
                                           array_t *toCareSetArray)
{
  int i;
  array_t *toCareSetArrayRV = NIL(array_t);
  mdd_t *toCareSet, *toCareSetRV;
  mdd_t *imageRV, *imageDV;
  int useToCareSetFlag;
  Iwls95Info_t *info = (Iwls95Info_t *)imageInfo->methodData;

  if (toCareSetArray && info->allowPartialImage)
    useToCareSetFlag = 1;
  else
    useToCareSetFlag = 0;

  if (useToCareSetFlag) {
    toCareSetArrayRV = array_alloc(mdd_t *, 0);
    arrayForEachItem(bdd_t *, toCareSetArray, i, toCareSet) {
      toCareSetRV = ImgSubstitute(toCareSet, functionData, Img_D2R_c);
      array_insert(bdd_t *, toCareSetArrayRV, i, toCareSetRV);
    }
  }

  imageRV = ImgImageInfoComputeFwdIwls95(functionData,
                                         imageInfo,
                                         fromLowerBound,
                                         fromUpperBound,
                                         toCareSetArrayRV);

  imageDV = ImgSubstitute(imageRV, functionData, Img_R2D_c);
  mdd_free(imageRV);
  if (useToCareSetFlag)
    mdd_array_free(toCareSetArrayRV);
  return imageDV;
}


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

  Synopsis    [Computes the backward image of domainSubset.]

  Description [First bdd_cofactor is used to compute the simplified set of
  states to compute the image of. Next, this simplified set is multiplied with
  the relation arrays given in the bwdClusteredRelationArray in order and the
  variables are quantified according to the schedule in
  bwdArraySmoothVarBddArray.]

  SideEffects []

  SeeAlso     [ImgImageInfoComputeBwdWithDomainVarsIwls95]

******************************************************************************/
mdd_t *
ImgImageInfoComputeBwdIwls95(ImgFunctionData_t *functionData,
                             Img_ImageInfo_t *imageInfo,
                             mdd_t *fromLowerBound,
                             mdd_t *fromUpperBound,
                             array_t *toCareSetArray)
{
  mdd_t *image, *domainSubset, *simplifiedImage;
  mdd_manager *mddManager;
  Iwls95Info_t *info = (Iwls95Info_t *)imageInfo->methodData;
  boolean partial;
  
  if (info->bwdClusteredRelationArray == NIL(array_t)) {
    fprintf(vis_stderr, "** img error: The data structure has not been ");
    fprintf(vis_stderr, "initialized for backward image computation\n");
    return NIL(mdd_t);
  }
  mddManager = Part_PartitionReadMddManager(functionData->mddNetwork);
  domainSubset = bdd_between(fromLowerBound, fromUpperBound);

  assert(mddManager != NIL(mdd_manager));

  /* Check if we can use the stored simplified relations */
  if (info->careSetArray == NIL(array_t) ||
      !mdd_array_equal(info->careSetArray, toCareSetArray)){
    FreeClusteredCofactoredRelationArray(info);
    info->careSetArray = mdd_array_duplicate(toCareSetArray);
    info->bwdClusteredCofactoredRelationArray =
      ImgMinimizeImageArrayWithCareSetArray(info->bwdClusteredRelationArray,
                                            toCareSetArray,
                                            Img_DefaultMinimizeMethod_c, TRUE,
                                            (info->option->verbosity > 0),
                                            Img_Backward_c);
  }

  /* if the partial images are allowed, compute the care set (or its superset)
   * so that partial products can be minimized with respect to it.
   */
  partial = info->allowPartialImage;
  info->wasPartialImage = FALSE;

  /* bias is used to bias the intermediate results towards care states */
  if ((info->PIoption)->partialImageApproxMethod == BDD_APPROX_BIASED_RUA) {
    /* assume that this array has only one mdd */
    mdd_t *toCareSet = array_fetch(mdd_t *, toCareSetArray, 0);
    (info->PIoption)->bias = mdd_dup(toCareSet);
  }


  image = BddLinearAndSmooth(mddManager, domainSubset,
                             info->bwdClusteredCofactoredRelationArray,
                             info->bwdArraySmoothVarBddArray,
                             info->bwdSmoothVarCubeArray,
                             (info->option)->verbosity,
                             info->PIoption, &partial, info->lazySiftFlag);


  /* check that no partial image was computed if not allowed */
  assert(!((!info->allowPartialImage) && partial));

  /* if partial image computed, and if no new states, recompute
   * image with relaxed parameters.
   */
  if (partial && bdd_leq_array(image, toCareSetArray, 1, 0)) {
    bdd_free(image);
    partial = info->allowPartialImage;
    image = RecomputeImageIfNecessary(functionData,
                                      mddManager,
                                      domainSubset,
                                      info->bwdClusteredRelationArray,
                                      info->bwdArraySmoothVarBddArray,
                                      info->bwdSmoothVarCubeArray,
                                      (info->option)->verbosity,
                                      info->PIoption,
                                      toCareSetArray,
                                      &partial, info->lazySiftFlag);
  }

  /* free the bias function created. */
  if ((info->PIoption)->partialImageApproxMethod == BDD_APPROX_BIASED_RUA) {
    mdd_free((info->PIoption)->bias);
    (info->PIoption)->bias = NIL(mdd_t);
  }
  /* indicate whether this image computation was partial or exact. */
  info->wasPartialImage = partial;
  info->allowPartialImage = FALSE;

  mdd_free(domainSubset);
  simplifiedImage = bdd_minimize_array(image, toCareSetArray);
  bdd_free(image);
  return simplifiedImage;
}

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

  Synopsis    [Computes the backward image of a set of states.]

  Description [Identical to ImgImageInfoComputeBwdIwls95 except in the
  fromLowerBound and fromUpperBound, domainVars are substituted by
  rangeVars.]

  SideEffects []

  SeeAlso     [ImgImageInfoComputeBwdIwls95]

******************************************************************************/
mdd_t *
ImgImageInfoComputeBwdWithDomainVarsIwls95(ImgFunctionData_t *functionData,
                                           Img_ImageInfo_t *imageInfo,
                                           mdd_t *fromLowerBound,
                                           mdd_t *fromUpperBound,
                                           array_t *toCareSetArray)
{
  mdd_t *fromLowerBoundRV;
  mdd_t *fromUpperBoundRV;
  mdd_t *image;

  fromLowerBoundRV = ImgSubstitute(fromLowerBound, functionData, Img_D2R_c);
  if (mdd_equal(fromLowerBound, fromUpperBound))
    fromUpperBoundRV = fromLowerBoundRV;
  else
    fromUpperBoundRV = ImgSubstitute(fromUpperBound, functionData, Img_D2R_c);

  image =  ImgImageInfoComputeBwdIwls95(functionData,
                                        imageInfo,
                                        fromLowerBoundRV,
                                        fromUpperBoundRV,
                                        toCareSetArray);

  mdd_free(fromLowerBoundRV);
  if (!mdd_equal(fromLowerBound, fromUpperBound))
    mdd_free(fromUpperBoundRV);
  return image;
}

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

  Synopsis    [Frees the memory associated with imageInfo.]

  Description [Frees the memory associated with imageInfo.]

  SideEffects []

  SeeAlso [ImgImageInfoInitializeIwls95]

******************************************************************************/
void
ImgImageInfoFreeIwls95(void *methodData)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;

  if (info == NIL(Iwls95Info_t)) {
    fprintf(vis_stderr, "** img error: Trying to free the NULL image info\n");
    return;
  }
  if (info->bitRelationArray != NIL(array_t)) {
    mdd_array_free(info->bitRelationArray);
  }
  if (info->quantifyVarMddIdArray != NIL(array_t)) {
    array_free(info->quantifyVarMddIdArray);
  }
  if (info->fwdOptClusteredRelationArray != NIL(array_t)) {
    mdd_array_free(info->fwdOptClusteredRelationArray);
    mdd_array_array_free(info->fwdOptArraySmoothVarBddArray);
  }
  if (info->fwdOptSmoothVarCubeArray)
    mdd_array_free(info->fwdOptSmoothVarCubeArray);
  if (info->fwdClusteredRelationArray != NIL(array_t)) {
    mdd_array_free(info->fwdClusteredRelationArray);
    mdd_array_array_free(info->fwdArraySmoothVarBddArray);
  }
  if (info->fwdSmoothVarCubeArray)
    mdd_array_free(info->fwdSmoothVarCubeArray);
  if (info->bwdClusteredRelationArray != NIL(array_t)) {
    mdd_array_free(info->bwdClusteredRelationArray);
    mdd_array_array_free(info->bwdArraySmoothVarBddArray);
  }
  if (info->bwdSmoothVarCubeArray)
    mdd_array_free(info->bwdSmoothVarCubeArray);
  FreeClusteredCofactoredRelationArray(info);
  
  if (info->option != NULL) ImgFreeTrmOptions(info->option);
  if (info->PIoption != NULL) FREE(info->PIoption);

  if (info->savedArraySmoothVarBddArray != NIL(array_t))
    mdd_array_array_free(info->savedArraySmoothVarBddArray);
  if (info->savedSmoothVarCubeArray != NIL(array_t))
    mdd_array_free(info->savedSmoothVarCubeArray);

  FREE(info);
}

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

  Synopsis    [Prints information about the IWLS95 method.]

  Description [This function is used to obtain the information about
  the parameters used to initialize the imageInfo. If the file pointer
  argument is not NULL, the options are printed out. The shared size
  of the transition relation is printed out on vis_stdout.]

  SideEffects []

******************************************************************************/
void
ImgImageInfoPrintMethodParamsIwls95(void *methodData, FILE *fp)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;
  if (fp == NULL) return;
  PrintOption(info->method, info->option, fp);
  if (info->fwdClusteredRelationArray != NIL(array_t)) {
    (void) fprintf(fp, "Shared size of transition relation for forward image computation is %10ld BDD nodes (%-4d components)\n",
                   bdd_size_multiple(info->fwdClusteredRelationArray),
                   array_n(info->fwdClusteredRelationArray));
  }
  if (info->bwdClusteredRelationArray != NIL(array_t)) {
    (void) fprintf(fp, "Shared size of transition relation for backward image computation is %10ld BDD nodes (%-4d components)\n",
                   bdd_size_multiple(info->bwdClusteredRelationArray),
                   array_n(info->bwdClusteredRelationArray));
  }
}

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

  Synopsis    [Returns the st_table containing the bdd id of the support
  variables of the function.]

  Description [Returns the st_table containing the bdd id of the support
  variables of the function.]

  SideEffects []

******************************************************************************/
st_table *
ImgBddGetSupportIdTable(bdd_t *function)
{
  st_table *supportTable;
  var_set_t *supportVarSet;
  int i;

  supportTable = st_init_table(st_numcmp, st_numhash);
  supportVarSet = bdd_get_support(function);
  for(i=0; i<supportVarSet->n_elts; i++) {
    if (var_set_get_elt(supportVarSet, i) == 1) {
      st_insert(supportTable, (char *)(long) i, NIL(char));
    }
  }
  var_set_free(supportVarSet);
  return supportTable;
}

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

  Synopsis    [Checks if the last image/preimage was partial.]

  Description [Checks if the last image/preimage was partial. This
  flag is valid only for the last image/preimage computation. Reset
  before every image/preimage computation. ]

  SideEffects [Reset before every image/preimage computation.]

  SeeAlso     [ImgImageAllowPartialImageIwls95]

******************************************************************************/
boolean
ImgImageWasPartialIwls95(void *methodData)
{
  Iwls95Info_t *info = (Iwls95Info_t *) methodData;
  return (info->wasPartialImage);
}

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

  Synopsis    [Sets the flag to allow partial images.]

  Description [Sets the flag to allow partial images.  Default is to
  not allow partial image computations. This flag is reset at the end
  of every image/preimage computation i.e., has to be specified before
  every image/preimage computation if partial image is desired. ]

  SideEffects [Flag is reset at the end of every image/preimage
  computation.]

  SeeAlso     [ImgImageWasPartialIwls95]

******************************************************************************/
void
ImgImageAllowPartialImageIwls95(void *methodData, boolean value)
{
  Iwls95Info_t *info = (Iwls95Info_t *) methodData;
  info->allowPartialImage = value;
}

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

  Synopsis    [Restores original transition relation from saved.]

  Description [Restores original transition relation from saved.]

  SideEffects []

  SeeAlso[ImgDuplicateTransitionRelationIwls95]
******************************************************************************/
void
ImgRestoreTransitionRelationIwls95(Img_ImageInfo_t *imageInfo,
        Img_DirectionType directionType)
{
  array_t *relationArray;
  Iwls95Info_t *iwlsInfo = (Iwls95Info_t *) imageInfo->methodData;
  ImgFunctionData_t *functionData = &(imageInfo->functionData);
  mdd_manager *mgr = NIL(mdd_manager);
  
  mgr =  Part_PartitionReadMddManager(imageInfo->functionData.mddNetwork);
  
  if (directionType == Img_Both_c) {
    fprintf(vis_stderr,
            "** img error : can't restore fwd and bwd at the same time.\n");
    return;
  }

  relationArray = ImgGetTransitionRelationIwls95(iwlsInfo, directionType);
  mdd_array_free(relationArray);
  ImgUpdateTransitionRelationIwls95(iwlsInfo,
                                    (array_t *)imageInfo->savedRelation,
                                    directionType);
  imageInfo->savedRelation = NIL(void);

  assert(iwlsInfo->savedArraySmoothVarBddArray);
  if (directionType == Img_Forward_c) {
    mdd_array_array_free(iwlsInfo->fwdArraySmoothVarBddArray);
    iwlsInfo->fwdArraySmoothVarBddArray = iwlsInfo->savedArraySmoothVarBddArray;
    if (functionData->createVarCubesFlag) {
      mdd_array_free(iwlsInfo->fwdSmoothVarCubeArray);
      iwlsInfo->fwdSmoothVarCubeArray = iwlsInfo->savedSmoothVarCubeArray;
    }
  } else {
    mdd_array_array_free(iwlsInfo->bwdArraySmoothVarBddArray);
    iwlsInfo->bwdArraySmoothVarBddArray = iwlsInfo->savedArraySmoothVarBddArray;
    if (functionData->createVarCubesFlag) {
      mdd_array_free(iwlsInfo->bwdSmoothVarCubeArray);
      iwlsInfo->bwdSmoothVarCubeArray = iwlsInfo->savedSmoothVarCubeArray;
    }
  }
  iwlsInfo->savedArraySmoothVarBddArray = NIL(array_t);
  iwlsInfo->savedSmoothVarCubeArray = NIL(array_t);

  if(directionType == Img_Backward_c || directionType == Img_Both_c)
    ResetClusteredCofactoredRelationArray(mgr, iwlsInfo);
}

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

  Synopsis    [Duplicates transition relation.]

  Description [Duplicates transition relation. Assumes that the
  Transition Relation is an array.]

  SideEffects []

  SeeAlso[ImgRestoreTransitionRelationIwls95]
******************************************************************************/
void
ImgDuplicateTransitionRelationIwls95(Img_ImageInfo_t *imageInfo,
                                     Img_DirectionType directionType)
{
  array_t *relationArray, *copiedRelationArray;
  Iwls95Info_t *iwlsInfo = (Iwls95Info_t *) imageInfo->methodData;
  ImgFunctionData_t *functionData = &(imageInfo->functionData);

  if (directionType == Img_Both_c) {
    fprintf(vis_stderr,
            "** img error : can't duplicate fwd and bwd at the same time.\n");
    return;
  }

  relationArray = ImgGetTransitionRelationIwls95(iwlsInfo, directionType);
  copiedRelationArray = BddArrayDup(relationArray);
  assert(!imageInfo->savedRelation);
  imageInfo->savedRelation = (void *)relationArray;
  ImgUpdateTransitionRelationIwls95(iwlsInfo,
                                    copiedRelationArray, directionType);

  assert(!iwlsInfo->savedArraySmoothVarBddArray);
  if (directionType == Img_Forward_c) {
    iwlsInfo->savedArraySmoothVarBddArray = iwlsInfo->fwdArraySmoothVarBddArray;
    iwlsInfo->fwdArraySmoothVarBddArray =
      CopyArrayBddArray(iwlsInfo->savedArraySmoothVarBddArray);
    if (functionData->createVarCubesFlag) {
      iwlsInfo->savedSmoothVarCubeArray = iwlsInfo->fwdSmoothVarCubeArray;
      iwlsInfo->fwdSmoothVarCubeArray =
        BddArrayDup(iwlsInfo->savedSmoothVarCubeArray);
    }
  } else {
    iwlsInfo->savedArraySmoothVarBddArray = iwlsInfo->bwdArraySmoothVarBddArray;
    iwlsInfo->bwdArraySmoothVarBddArray =
      CopyArrayBddArray(iwlsInfo->savedArraySmoothVarBddArray);
    if (functionData->createVarCubesFlag) {
      iwlsInfo->savedSmoothVarCubeArray = iwlsInfo->bwdSmoothVarCubeArray;
      iwlsInfo->bwdSmoothVarCubeArray =
        BddArrayDup(iwlsInfo->savedSmoothVarCubeArray);
    }
  }
}


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

  Synopsis    [Minimizes transition relation with given set of constraints.]

  Description [Minimizes transition relation with given set of
  constraints. This invalidates the bwdclusteredcofactoredrelationarray.

  <p>RB: Can Kavita or In-Ho please explain the difference between these four
  related functions?? (other three in SeeAlso)

  The constraint should be in terms of present-state and input variables only.

  This function takes care that the quantification schedule is correct
  afterwards.  If reorderClusters is true, the clusters are reordered after
  minimization. 
  ]

  SideEffects []

  SeeAlso [ImgMinimizeTransitionRelationIwls95,
  ImgApproximateTransitionRelationIwls95, ImgAbstractTransitionRelationIwls95,
  Img_MinimizeImage] 
******************************************************************************/
void
ImgMinimizeTransitionRelationIwls95(
  void *methodData,
  ImgFunctionData_t  *functionData,
  array_t *constrainArray,
  Img_MinimizeType minimizeMethod,
  Img_DirectionType directionType,
  boolean reorderClusters,
  int printStatus)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;

  if (minimizeMethod == Img_DefaultMinimizeMethod_c)
    minimizeMethod = Img_ReadMinimizeMethod();  
  
  if (directionType == Img_Forward_c || directionType == Img_Both_c) {
    ImgMinimizeImageArrayWithCareSetArrayInSitu(
      info->fwdClusteredRelationArray, constrainArray, minimizeMethod, TRUE,
      (printStatus == 2), Img_Forward_c); 

    if(reorderClusters)
      ReorderPartitionedTransitionRelation(info, functionData,
                                           Img_Forward_c); 
    else{
      /* any minimization method other than restrict may introduce new
         variables, which invalidates the quantification schedule */
      assert(minimizeMethod != Img_DefaultMinimizeMethod_c);
      if(minimizeMethod != Img_Restrict_c)
        UpdateQuantificationSchedule(info, functionData, Img_Forward_c);
    }
     assert(CheckQuantificationSchedule(info->fwdClusteredRelationArray,
                                     info->fwdArraySmoothVarBddArray));
  }

  if (directionType == Img_Backward_c || directionType == Img_Both_c) {
    ImgMinimizeImageArrayWithCareSetArrayInSitu(
      info->bwdClusteredRelationArray, constrainArray, minimizeMethod, TRUE,
      (printStatus == 2), Img_Backward_c); 

    if(reorderClusters)
      ReorderPartitionedTransitionRelation(info, functionData,
                                           Img_Backward_c); 
    /* Minimization can only introduce PS variables, which need not be
       quantified out, hence no need to recompute quantification schedule
       in backward case. */
    assert(CheckQuantificationSchedule(info->bwdClusteredRelationArray,
                                       info->bwdArraySmoothVarBddArray))
    }

  if(directionType == Img_Backward_c || directionType == Img_Both_c)
    FreeClusteredCofactoredRelationArray(info);
  
  return;
}


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

  Synopsis    [Abstracts out given variables from transition relation.]

  Description [Abstracts out given variables from transition relation.
  abstractVars should be an array of bdd variables.  This invalidates the
  bwdclusteredcofactoredrelationarray. ]

  SideEffects []
  
  SeeAlso [ImgMinimizeTransitionRelationIwls95,
  ImgApproximateTransitionRelationIwls95, ImgAbstractTransitionRelationIwls95]

******************************************************************************/
void
ImgAbstractTransitionRelationIwls95(Img_ImageInfo_t *imageInfo,
        array_t *abstractVars, mdd_t *abstractCube,
        Img_DirectionType directionType, int printStatus)
{
  Iwls95Info_t *info = (Iwls95Info_t *)imageInfo->methodData;
  int i, fwd_size, bwd_size;
  bdd_t *relation, *abstractedRelation;
  
  if (!abstractVars || array_n(abstractVars) == 0)
    return;

  fwd_size = bwd_size = 0;
  if (printStatus) {
    if (directionType == Img_Forward_c || directionType == Img_Both_c)
      fwd_size = bdd_size_multiple(info->fwdClusteredRelationArray);
    if (directionType == Img_Backward_c || directionType == Img_Both_c)
      bwd_size = bdd_size_multiple(info->bwdClusteredRelationArray);
  } 
  if (directionType == Img_Forward_c || directionType == Img_Both_c) {
    arrayForEachItem(bdd_t*, info->fwdClusteredRelationArray, i, relation) {
      if (abstractCube)
        abstractedRelation = bdd_smooth_with_cube(relation, abstractCube);
      else
        abstractedRelation = bdd_smooth(relation, abstractVars);
      if (printStatus == 2) {
        (void) fprintf(vis_stdout,
          "IMG Info: abstracted fwd relation %d => %d in component %d\n",
                         bdd_size(relation), bdd_size(abstractedRelation), i);
      }
      mdd_free(relation);
      array_insert(bdd_t*, info->fwdClusteredRelationArray, i,
                   abstractedRelation);
    }
  }
  if (directionType == Img_Backward_c || directionType == Img_Both_c) {
    arrayForEachItem(bdd_t*, info->bwdClusteredRelationArray, i, relation) {
      if (abstractCube)
        abstractedRelation = bdd_smooth_with_cube(relation, abstractCube);
      else
        abstractedRelation = bdd_smooth(relation, abstractVars);
      if (printStatus == 2) {
        (void) fprintf(vis_stdout,
          "IMG Info: abstracted bwd relation %d => %d in component %d\n",
                         bdd_size(relation), bdd_size(abstractedRelation), i);
      }
      mdd_free(relation);
      array_insert(bdd_t*, info->bwdClusteredRelationArray, i,
                   abstractedRelation);
    }
  }
  if (printStatus) {
    if (directionType == Img_Forward_c || directionType == Img_Both_c) {
      (void) fprintf(vis_stdout,
     "IMG Info: abstracted fwd relation [%d => %ld] with %d components\n",
                     fwd_size,
                     bdd_size_multiple(info->fwdClusteredRelationArray),
                     array_n(info->fwdClusteredRelationArray));
    }
    if (directionType == Img_Backward_c || directionType == Img_Both_c) {
      (void) fprintf(vis_stdout,
     "IMG Info: abstracted bwd relation [%d => %ld] with %d components\n",
                     bwd_size,
                     bdd_size_multiple(info->bwdClusteredRelationArray),
                     array_n(info->bwdClusteredRelationArray));
    }
  }

  if(directionType == Img_Backward_c || directionType == Img_Both_c) 
    FreeClusteredCofactoredRelationArray(info);
}


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

  Synopsis    [Approximate transition relation.]

  Description [Approximate transition relation.  This invalidates the
  bwdclusteredcofactoredrelationarray. ]

  SideEffects [ImgMinimizeTransitionRelationIwls95,
  ImgApproximateTransitionRelationIwls95, ImgAbstractTransitionRelationIwls95,
  Img_ApproximateImage] 

******************************************************************************/
int
ImgApproximateTransitionRelationIwls95(mdd_manager *mgr, void *methodData,
                    bdd_approx_dir_t approxDir, bdd_approx_type_t approxMethod,
                    int approxThreshold, double approxQuality,
                    double approxQualityBias,
                    Img_DirectionType directionType, mdd_t *bias)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;
  int i;
  bdd_t *relation, *approxRelation;
  int   unchanged = 0;
  
  if (directionType == Img_Forward_c || directionType == Img_Both_c) {
    arrayForEachItem(bdd_t*, info->fwdClusteredRelationArray, i, relation) {
      approxRelation = Img_ApproximateImage(mgr, relation,
                                            approxDir, approxMethod,
                                            approxThreshold, approxQuality,
                                            approxQualityBias, bias);
      if (bdd_is_tautology(approxRelation, 1)) {
        fprintf(vis_stdout,
                "** img warning : bdd_one from subsetting in fwd[%d].\n", i);
        mdd_free(approxRelation);
        unchanged++;
      } else if (bdd_is_tautology(approxRelation, 0)) {
        fprintf(vis_stdout,
                "** img warning : bdd_zero from subsetting in fwd[%d].\n", i);
        mdd_free(approxRelation);
        unchanged++;
      } else if (mdd_equal(approxRelation, relation)) {
        mdd_free(approxRelation);
        unchanged++;
      } else {
        mdd_free(relation);
        array_insert(bdd_t*, info->fwdClusteredRelationArray, i,
                     approxRelation);
      }
    }
  }
  if (directionType == Img_Backward_c || directionType == Img_Both_c) {
    arrayForEachItem(bdd_t*, info->bwdClusteredRelationArray, i, relation) {
      approxRelation = Img_ApproximateImage(mgr, relation,
                                            approxDir, approxMethod,
                                            approxThreshold, approxQuality,
                                            approxQualityBias, bias);
      if (bdd_is_tautology(approxRelation, 1)) {
        fprintf(vis_stdout,
                "** img warning : bdd_one from subsetting in bwd[%d].\n", i);
        mdd_free(approxRelation);
        unchanged++;
      } else if (bdd_is_tautology(approxRelation, 0)) {
        fprintf(vis_stdout,
                "** img warning : bdd_zero from subsetting in bwd[%d].\n", i);
        mdd_free(approxRelation);
        unchanged++;
      } else if (mdd_equal(approxRelation, relation)) {
        mdd_free(approxRelation);
        unchanged++;
      } else {
        mdd_free(relation);
        array_insert(bdd_t*, info->bwdClusteredRelationArray, i,
                     approxRelation);
      }
    }
  }
  
  if(directionType == Img_Backward_c || directionType == Img_Both_c)
    FreeClusteredCofactoredRelationArray(info);

  return(unchanged);
}


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

  Synopsis    [Returns current transition relation.]

  Description [Returns current transition relation.]

  SideEffects []

******************************************************************************/
array_t *
ImgGetTransitionRelationIwls95(void *methodData,
                               Img_DirectionType directionType)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;

  if (directionType == Img_Forward_c)
    return(info->fwdClusteredRelationArray);
  else if (directionType == Img_Backward_c)
    return(info->bwdClusteredRelationArray);
  return(NIL(array_t));
}


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

  Synopsis    [Overwrites transition relation with given.]

  Description [Overwrites transition relation with given.]

  SideEffects []

******************************************************************************/
void
ImgUpdateTransitionRelationIwls95(void *methodData, array_t *relationArray,
                                  Img_DirectionType directionType)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;

  if (directionType == Img_Forward_c)
    info->fwdClusteredRelationArray = relationArray;
  else if (directionType == Img_Backward_c)
    info->bwdClusteredRelationArray = relationArray;
}


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

  Synopsis    [Replace ith partitioned transition relation.]

  Description [Replace ith partitioned transition relation.]

  SideEffects []

******************************************************************************/
void
ImgReplaceIthPartitionedTransitionRelationIwls95(void *methodData,
        int i, mdd_t *relation, Img_DirectionType directionType)
{
  array_t       *relationArray;
  mdd_t         *old;

  relationArray = ImgGetTransitionRelationIwls95(methodData, directionType);
  if (!relationArray)
    return;
  old = array_fetch(mdd_t *, relationArray, i);
  mdd_free(old);
  array_insert(mdd_t *, relationArray, i, relation);
}


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

  Synopsis [Frees current transition relation and associated
  quantification schedule for the given direction.]

  Description [Frees current transition relation and associated
  quantification schedule for the given direction.]

  SideEffects []

******************************************************************************/
void
ImgImageFreeClusteredTransitionRelationIwls95(void *methodData,
                               Img_DirectionType directionType)
{
  Iwls95Info_t *info = (Iwls95Info_t *)methodData;

  if (directionType == Img_Forward_c || directionType == Img_Both_c) {
    if (info->fwdClusteredRelationArray != NIL(array_t)) {
      mdd_array_free(info->fwdClusteredRelationArray);
      info->fwdClusteredRelationArray = NIL(array_t);
    }
    if (info->fwdArraySmoothVarBddArray != NIL(array_t)) {
      mdd_array_array_free(info->fwdArraySmoothVarBddArray);
      info->fwdArraySmoothVarBddArray = NIL(array_t);
    }
    if (info->fwdSmoothVarCubeArray != NIL(array_t)) {
      mdd_array_free(info->fwdSmoothVarCubeArray);
      info->fwdSmoothVarCubeArray = NIL(array_t);
    }
  }
  if (directionType == Img_Backward_c || directionType == Img_Both_c) {
    if (info->bwdClusteredRelationArray != NIL(array_t)) {
      mdd_array_free(info->bwdClusteredRelationArray);
      info->bwdClusteredRelationArray = NIL(array_t);
    }
    if (info->bwdArraySmoothVarBddArray != NIL(array_t)) {
      mdd_array_array_free(info->bwdArraySmoothVarBddArray);
      info->bwdArraySmoothVarBddArray = NIL(array_t);
    }
    if (info->bwdSmoothVarCubeArray != NIL(array_t)) {
      mdd_array_free(info->bwdSmoothVarCubeArray);
      info->bwdSmoothVarCubeArray = NIL(array_t);
    }
  }
}



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

  Synopsis [Constrains/Adds dont-cares to the bit relation and creates
  a new transition relation.]

  Description [Constrains/adds dont-cares to the bit relation and creates a new
  transition relation. First, the existing transition relation is freed.  The
  bit relation is created if it isn't stored already. The bit relation is then
  copied into the clustered relation field of the given direction.  Then, it is
  constrained by the given constraint (underapprox) or the complement of the
  constraint is added (overapprox) to the bit relation.  Then the modified bit
  relation is clustered.  The underApprox flag indicates whether to create an
  over or an underapproximation of the transition relation (TRUE for
  underapproximation, FALSE for overapproximation).  Although the procedure
  Img_MinimizeImage is used, the minimizeMethod flag for underapproximations
  should be generated by Img_GuidedSearchReadUnderApproxMinimizeMethod.  The
  clean up flag indicates freeing of the bit relations.  If the forceReorder
  flag is set, variable reordering is fired after freeing old relations and
  before creating new relations.]

  SideEffects [Frees current transition relation and quantification
  schedule. Frees bit relations if indicated. ]

  SeeAlso [Img_MultiwayLinearAndSmooth ImgImageInfoInitializeMono
  Img_MinimizeImage]
  
******************************************************************************/
void
ImgImageConstrainAndClusterTransitionRelationIwls95OrMlp(
  Img_ImageInfo_t *imageInfo,
  Img_DirectionType direction,
  mdd_t *constrain,
  Img_MinimizeType minimizeMethod,
  boolean underApprox,
  boolean cleanUp,
  boolean forceReorder,
    int printStatus)
{
  ImgFunctionData_t *functionData = &(imageInfo->functionData);
  Iwls95Info_t *info = (Iwls95Info_t *)imageInfo->methodData;
  array_t *constrainArray;
  mdd_manager *mddManager = Part_PartitionReadMddManager(functionData->mddNetwork);
  array_t *relationArray    = NIL(array_t);
  array_t *relationArrayFwd = NIL(array_t);
  array_t *relationArrayBwd = NIL(array_t);

  assert(imageInfo->methodType == Img_Iwls95_c ||
         imageInfo->methodType == Img_Mlp_c);

  /* free existing relation */
  ImgImageFreeClusteredTransitionRelationIwls95(imageInfo->methodData,
                                                direction);
  if (forceReorder) bdd_reorder(mddManager);
  /* create a bit relation and quantifiable variables (PIs and
     intermediate variables) if necessary */
  CreateBitRelationArray(info, functionData);

  /* make a work copy of the bit array */
  relationArray = BddArrayDup(info->bitRelationArray);

  if (cleanUp) FreeBitRelationArray(info);
  
  /* apply the constraint to the bit relation */
  if (constrain) {
    constrainArray = array_alloc(mdd_t *, 1);
    array_insert(mdd_t *, constrainArray, 0, constrain);
    ImgMinimizeImageArrayWithCareSetArrayInSitu(relationArray,
                                                constrainArray,
                                                minimizeMethod, underApprox,
                                                printStatus, direction);
    array_free(constrainArray);
  }

  if (direction == Img_Forward_c) {
    relationArrayFwd = relationArray;
  } else if (direction == Img_Backward_c) {
    relationArrayBwd = relationArray;
  } else if (direction == Img_Both_c) {
    relationArrayFwd = relationArray;
    relationArrayBwd = BddArrayDup(relationArray);
  } else {
    assert(0);
  }
  relationArray = NIL(array_t);
  
  /* return order and recluster the relation */
  if (direction == Img_Forward_c || direction == Img_Both_c) {
    Img_ClusterRelationArray(mddManager,
                             imageInfo->methodType,
                             direction,
                             relationArrayFwd,
                             functionData->domainVars,
                             functionData->rangeVars,
                             info->quantifyVarMddIdArray,
                             &info->fwdClusteredRelationArray,
                             &info->fwdArraySmoothVarBddArray,
                             NULL,
                             1);
    if (functionData->createVarCubesFlag) {
      info->fwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->fwdArraySmoothVarBddArray);
    }
  }
  
  if (direction == Img_Backward_c || direction == Img_Both_c) {
    Img_ClusterRelationArray(mddManager,
                             imageInfo->methodType,
                             direction,
                             relationArrayBwd,
                             functionData->domainVars,
                             functionData->rangeVars,
                             info->quantifyVarMddIdArray,
                             &info->bwdClusteredRelationArray,
                             &info->bwdArraySmoothVarBddArray,
                             NULL,
                             1);
    if (functionData->createVarCubesFlag) {
      info->bwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->bwdArraySmoothVarBddArray);
    }
  }

  if (cleanUp) {
    array_free(info->quantifyVarMddIdArray);
    info->quantifyVarMddIdArray = NIL(array_t);
  }

  /* Print information */
  if (info->option->verbosity > 0){
    fprintf(vis_stdout,"Computing Image Using %s technique.\n",
            imageInfo->methodType == Img_Mlp_c ? "MLP" : "IWLS95");
    fprintf(vis_stdout,"Total # of domain binary variables = %3d\n",
            array_n(functionData->domainBddVars));
    fprintf(vis_stdout,"Total # of range binary variables = %3d\n",
            array_n(functionData->rangeBddVars));
    if (info->quantifyVarMddIdArray) {
      array_t *bddArray;

      bddArray = mdd_id_array_to_bdd_array(mddManager,
                                                  info->quantifyVarMddIdArray);
      fprintf(vis_stdout,"Total # of quantify binary variables = %3d\n",
              array_n(bddArray));
      mdd_array_free(bddArray);
    } else {
      fprintf(vis_stdout,"Total # of quantify binary variables = %3d\n",
              array_n(functionData->quantifyBddVars));
    }
    if (((direction == Img_Forward_c) || (direction == Img_Both_c)) &&
        (info->fwdClusteredRelationArray != NIL(array_t))) {
      (void) fprintf(vis_stdout,
                     "Shared size of transition relation for forward image ");
      (void) fprintf(vis_stdout,
                     "computation is %10ld BDD nodes (%-4d components)\n",
                     bdd_size_multiple(info->fwdClusteredRelationArray),
                     array_n(info->fwdClusteredRelationArray));
    }
    if (((direction == Img_Backward_c) || (direction == Img_Both_c)) &&
      (info->bwdClusteredRelationArray != NIL(array_t))) {
      (void) fprintf(vis_stdout,
                     "Shared size of transition relation for backward image ");
      (void) fprintf(vis_stdout,
                     "computation is %10ld BDD nodes (%-4d components)\n",
                     bdd_size_multiple(info->bwdClusteredRelationArray),
                     array_n(info->bwdClusteredRelationArray));
    }
  }/* if verbosity > 0 */

  if(direction == Img_Backward_c || direction == Img_Both_c)
    ResetClusteredCofactoredRelationArray(mddManager, info);
}



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

  Synopsis    [Prints integers from an array.]

  Description [Prints integers from an array.]

  SideEffects []

******************************************************************************/
void
ImgPrintIntegerArray(array_t *idArray)
{
  int i;
  fprintf(vis_stdout,
    "**************** printing bdd ids from the bdd array ***********\n");
  fprintf(vis_stdout,"%d::\t", array_n(idArray));
  for (i=0;i<array_n(idArray); i++) {
    fprintf(vis_stdout," %d ", array_fetch(int, idArray, i));
  }
}


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

  Synopsis    [Prints the given parition.]

  Description [Prints the given parition.]

  SideEffects []

******************************************************************************/
void
ImgPrintPartition(graph_t *partition)
{
  vertex_t *vertex;
  lsList vertexList = g_get_vertices(partition);
  lsGen gen;
  st_table *vertexTable = st_init_table(st_ptrcmp, st_ptrhash);
  fprintf(vis_stdout,"\n");
  lsForEachItem(vertexList, gen, vertex) {
    if (lsLength(g_get_out_edges(vertex)) == 0)
      PrintPartitionRecursively(vertex,vertexTable,0);
  }
  st_free_table(vertexTable);
}


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

  Synopsis    [Prints info of the partitioned transition relation.]

  Description [Prints info of the partitioned transition relation.]

  SideEffects []

******************************************************************************/
void
ImgPrintPartitionedTransitionRelation(mdd_manager *mddManager,
                                array_t *relationArray,
                                array_t *arraySmoothVarBddArray)
{
  int   i, j, n;
  mdd_t *relation;
  mdd_t *bdd;
  array_t *smoothVarBddArray;
  array_t *bvar_list = mdd_ret_bvar_list(mddManager);
  var_set_t *vset;
  bvar_type bv;

  n = array_n(relationArray);

  fprintf(vis_stdout, "# of relations = %d\n", n);
  if (arraySmoothVarBddArray) {
    smoothVarBddArray = array_fetch(array_t *, arraySmoothVarBddArray, 0);
    fprintf(vis_stdout, "Early smooth(%d) = ", array_n(smoothVarBddArray));
    for (j = 0; j < array_n(smoothVarBddArray); j++) {
      bdd = array_fetch(mdd_t *, smoothVarBddArray, j);
      mdd_print_support_to_file(vis_stdout, " %s", bdd);
    }
    fprintf(vis_stdout, "\n");
  }
  for (i = 0; i < n; i++) {
    relation = array_fetch(mdd_t *, relationArray, i);
    fprintf(vis_stdout, "T[%d] : bdd_size = %d\n", i, bdd_size(relation));
    fprintf(vis_stdout, "  support(%d) = ",
            mdd_get_number_of_bdd_support(mddManager, relation));
    vset = bdd_get_support(relation);
    for (j = 0; j < array_n(bvar_list); j++) {
      if (var_set_get_elt(vset, j) == 1) {
        bv = array_fetch(bvar_type, bvar_list, j);
        mdd_print_support_to_file(vis_stdout, " %s", bv.node);
      }
    }
    var_set_free(vset);
    fprintf(vis_stdout, "\n");
    if (arraySmoothVarBddArray) {
      smoothVarBddArray = array_fetch(array_t *, arraySmoothVarBddArray, i + 1);
      fprintf(vis_stdout, "  smooth(%d) = ", array_n(smoothVarBddArray));
      for (j = 0; j < array_n(smoothVarBddArray); j++) {
        bdd = array_fetch(mdd_t *, smoothVarBddArray, j);
        mdd_print_support_to_file(vis_stdout, " %s", bdd);
      }
      fprintf(vis_stdout, "\n");
    }
  }
}


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

  Synopsis    [Eliminates dependent variables from transition relation.]

  Description [Eliminates dependent variables from a transition
  relation.  Returns a simplified copy of the given states if successful;
  NULL otherwise.]

  SideEffects [relationArray is also modified.]

  SeeAlso     []

******************************************************************************/
mdd_t *
ImgTrmEliminateDependVars(ImgFunctionData_t *functionData,
                          array_t *relationArray,
                          mdd_t *states, mdd_t **dependRelations,
                          int *nDependVars)
{
  int           i, j;
  int           howMany = 0;    /* number of latches that can be eliminated */
  mdd_t         *var, *newStates, *abs, *positive, *phi, *tmp, *relation;
  int           nSupports;    /* vars in the support of the state set */
  int           *candidates;    /* vars to be considered for elimination */
  double        minStates;
  var_set_t     *supportVarSet;
  graph_t       *mddNetwork;
  mdd_manager   *mddManager;

  mddNetwork = functionData->mddNetwork;
  mddManager = Part_PartitionReadMddManager(mddNetwork);

  if (dependRelations)
    *dependRelations = mdd_one(mddManager);
  newStates = mdd_dup(states);

  nSupports = 0;
  supportVarSet = bdd_get_support(newStates);
  for (i = 0; i < supportVarSet->n_elts; i++) {
    if (var_set_get_elt(supportVarSet, i) == 1)
      nSupports++;
  }
  candidates = ALLOC(int, nSupports);
  if (candidates == NULL) {
    var_set_free(supportVarSet);
    return(NULL);
  }
  nSupports = 0;
  for (i = 0; i < supportVarSet->n_elts; i++) {
    if (var_set_get_elt(supportVarSet, i) == 1) {
      candidates[nSupports] = i;
      nSupports++;
    }
  }
  var_set_free(supportVarSet);

  /* The signatures of the variables in a function are the number
  ** of minterms of the positive cofactors with respect to the
  ** variables themselves. */
  signatures = bdd_cof_minterm(newStates);
  if (signatures == NULL) {
    FREE(candidates);
    return(NULL);
  }
  /* We now extract a positive quantity which is higher for those
  ** variables that are closer to being essential. */
  minStates = signatures[nSupports];
  for (i = 0; i < nSupports; i++) {
    double z = signatures[i] / minStates - 1.0;
    signatures[i] = (z < 0.0) ? -z : z;    /* make positive */
  }
  qsort((void *)candidates, nSupports, sizeof(int),
      (int (*)(const void *, const void *))TrmSignatureCompare);
  FREE(signatures);

  /* Now process the candidates in the given order. */
  for (i = 0; i < nSupports; i++) {
    var = bdd_var_with_index(mddManager, candidates[i]);
    if (bdd_var_is_dependent(newStates, var)) {
      abs = bdd_smooth_with_cube(newStates, var);
      if (abs == NULL)
        return(NULL);
      positive = bdd_cofactor(newStates, var);
      if (positive == NULL)
        return(NULL);
      phi = Img_MinimizeImage(positive, abs, Img_DefaultMinimizeMethod_c, 1);
      if (phi == NULL)
        return(NULL);
      mdd_free(positive);
      if (bdd_size(phi) < IMG_MAX_DEP_SIZE) {
        howMany++;
        for (j = 0; j < array_n(relationArray); j++) {
          relation = array_fetch(mdd_t *, relationArray, j);
          if (dependRelations)
            tmp = bdd_smooth_with_cube(relation, var);
          else
            tmp = bdd_compose(relation, var, phi);
          mdd_free(relation);
          array_insert(mdd_t *, relationArray, j, tmp);
        }
        mdd_free(newStates);
        newStates = abs;
        if (dependRelations) {
          relation = mdd_xnor(var, phi);
          tmp = ImgSubstitute(relation, functionData, Img_R2D_c);
          mdd_free(relation);
          relation = mdd_and(*dependRelations, tmp, 1, 1);
          mdd_free(*dependRelations);
          *dependRelations = relation;
        }
      } else {
        mdd_free(abs);
      }
      mdd_free(phi);
    }
  }
  FREE(candidates);

  *nDependVars = howMany;
  return(newStates);
} /* end of ImgTrmEliminateDependVars */


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

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

  Synopsis    [Reset the clusteredcofactoredrelationarray to the
  clusteredrelationarray.]

  Description [Copy the clustered array into the clusteredcofactored array, and
  set the caresetarray to (1).  Assumes either both the caresetarray and the
  cofactoredarray are NIL, or that they both contain sensible data.]

  SideEffects []

******************************************************************************/
static void
ResetClusteredCofactoredRelationArray(
   mdd_manager *mddManager,
   Iwls95Info_t *info)
{
  int i;
  mdd_t *cluster;

  assert(info->bwdClusteredRelationArray != NIL(array_t));
  
  FreeClusteredCofactoredRelationArray(info);

  info->careSetArray = array_alloc(mdd_t *, 1);
  array_insert_last(mdd_t *, info->careSetArray, mdd_one(mddManager));
  info->bwdClusteredCofactoredRelationArray = array_alloc(mdd_t *, 0);

  arrayForEachItem(mdd_t *, info->bwdClusteredRelationArray, i, cluster){
    array_insert_last(mdd_t *, info->bwdClusteredCofactoredRelationArray,
                      mdd_dup(cluster)); 
  }
}


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

  Synopsis    [Frees the clusteredcofactoredrelationarray]


  Description [Frees the clusteredCofactoredRelationArray, and the
  careSetArray, and sets them to NIL to signify that they are invalid.]

  SideEffects []

******************************************************************************/
static void
FreeClusteredCofactoredRelationArray(
   Iwls95Info_t *info)
{
  /* one should only be NIL if the other is, too */
  assert((info->bwdClusteredCofactoredRelationArray != NIL(array_t)) ==
         (info->careSetArray != NIL(array_t)));
  
  if(info->bwdClusteredCofactoredRelationArray == NIL(array_t))
    return;
  
  mdd_array_free(info->bwdClusteredCofactoredRelationArray);
  info->bwdClusteredCofactoredRelationArray = NIL(array_t);
  mdd_array_free(info->careSetArray);
  info->careSetArray = NIL(array_t);
}


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

  Synopsis    [Gets the necessary options for computing the image and returns
               in the option structure.]

  Description [Gets the necessary options for computing the image and returns
               in the option structure.]

  SideEffects []

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

static ImgTrmOption_t *
TrmGetOptions(void)
{
  char *flagValue;
  ImgTrmOption_t *option;

  option = ALLOC(ImgTrmOption_t, 1);

  flagValue = Cmd_FlagReadByName("image_cluster_size");
  if (flagValue == NIL(char)) {
    option->clusterSize = 5000; /* the default value */
  }
  else {
    option->clusterSize = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_verbosity");
  if (flagValue == NIL(char)) {
    option->verbosity = 0; /* the default value */
  }
  else {
    option->verbosity = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_W1");
  if (flagValue == NIL(char)) {
    option->W1 = 6; /* the default value */
  }
  else {
    option->W1 = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_W2");
  if (flagValue == NIL(char)) {
    option->W2 = 1; /* the default value */
  }
  else {
    option->W2 = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_W3");
  if (flagValue == NIL(char)) {
    option->W3 = 1; /* the default value */
  }
  else {
    option->W3 = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_W4");
  if (flagValue == NIL(char)) {
    option->W4 = 2; /* the default value */
  }
  else {
    option->W4 = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_print_depend_matrix");
  if (flagValue == NIL(char)) {
    option->printDepMatrix = 0; /* the default value */
  }
  else {
    option->printDepMatrix = atoi(flagValue);
  }


  flagValue = Cmd_FlagReadByName("mlp_method");
  if (flagValue == NIL(char)) {
    option->mlpMethod = 0; /* the default value */
  }
  else {
    option->mlpMethod = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_cluster");
  if (flagValue == NIL(char)) {
    option->mlpCluster = 1; /* the default value */
  }
  else {
    option->mlpCluster = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_cluster_dynamic");
  if (flagValue == NIL(char)) {
    option->mlpClusterDynamic = 1; /* the default value */
  }
  else {
    option->mlpClusterDynamic = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_affinity_threshold");
  if (flagValue == NIL(char)) {
    option->mlpAffinityThreshold = 0.0; /* the default value */
  }
  else {
    sscanf(flagValue, "%f", &option->mlpAffinityThreshold);
  }

  flagValue = Cmd_FlagReadByName("mlp_cluster_quantify_vars");
  if (flagValue == NIL(char)) {
    option->mlpClusterQuantifyVars = 0; /* the default value */
  }
  else {
    option->mlpClusterQuantifyVars = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_cluster_sorted_list");
  if (flagValue == NIL(char)) {
    option->mlpClusterSortedList = 1; /* the default value */
  }
  else {
    option->mlpClusterSortedList = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_initial_cluster");
  if (flagValue == NIL(char)) {
    option->mlpInitialCluster = 0; /* the default value */
  }
  else {
    option->mlpInitialCluster = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_cs_first");
  if (flagValue == NIL(char)) {
    option->mlpCsFirst = 0; /* the default value */
  }
  else {
    option->mlpCsFirst = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_skip_rs");
  if (flagValue == NIL(char)) {
    option->mlpSkipRs = 0; /* the default value */
  }
  else {
    option->mlpSkipRs = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_skip_cs");
  if (flagValue == NIL(char)) {
    option->mlpSkipCs = 1; /* the default value */
  }
  else {
    option->mlpSkipCs = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_row_based");
  if (flagValue == NIL(char)) {
    option->mlpRowBased = 0; /* the default value */
  }
  else {
    option->mlpRowBased = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_reorder");
  if (flagValue == NIL(char)) {
    option->mlpReorder = 0; /* the default value */
  }
  else {
    option->mlpReorder = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_pre_reorder");
  if (flagValue == NIL(char)) {
    option->mlpPreReorder = 0; /* the default value */
  }
  else {
    option->mlpPreReorder = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_postprocess");
  if (flagValue == NIL(char)) {
    option->mlpPostProcess = 0; /* the default value */
  }
  else {
    option->mlpPostProcess = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_decompose_scc");
  if (flagValue == NIL(char)) {
    option->mlpDecomposeScc = 1; /* the default value */
  }
  else {
    option->mlpDecomposeScc = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_cluster_scc");
  if (flagValue == NIL(char)) {
    option->mlpClusterScc = 1; /* the default value */
  }
  else {
    option->mlpClusterScc = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_sort_scc");
  if (flagValue == NIL(char)) {
    option->mlpSortScc = 1; /* the default value */
  }
  else {
    option->mlpSortScc = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_sort_scc_mode");
  if (flagValue == NIL(char)) {
    option->mlpSortSccMode = 2; /* the default value */
  }
  else {
    option->mlpSortSccMode = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_cluster_merge");
  if (flagValue == NIL(char)) {
    option->mlpClusterMerge = 0; /* the default value */
  }
  else {
    option->mlpClusterMerge = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_multiway");
  if (flagValue == NIL(char)) {
    option->mlpMultiway = 0; /* the default value */
  }
  else {
    option->mlpMultiway = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_lambda_mode");
  if (flagValue == NIL(char)) {
    option->mlpLambdaMode = 0; /* the default value */
  }
  else {
    option->mlpLambdaMode = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_sorted_rows");
  if (flagValue == NIL(char)) {
    option->mlpSortedRows = 1; /* the default value */
  }
  else {
    option->mlpSortedRows = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_sorted_cols");
  if (flagValue == NIL(char)) {
    option->mlpSortedCols = 1; /* the default value */
  }
  else {
    option->mlpSortedCols = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_pre_swap_state_vars");
  if (flagValue == NIL(char)) {
    option->mlpPreSwapStateVars = 0; /* the default value */
  }
  else {
    option->mlpPreSwapStateVars = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_print_matrix");
  if (flagValue == NIL(char)) {
    option->mlpPrintMatrix = 0; /* the default value */
  }
  else {
    option->mlpPrintMatrix = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_write_order");
  if (flagValue == NIL(char)) {
    option->mlpWriteOrder = NIL(char); /* the default value */
  }
  else {
    option->mlpWriteOrder = util_strsav(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_verbosity");
  if (flagValue == NIL(char)) {
    option->mlpVerbosity = 0; /* the default value */
  }
  else {
    option->mlpVerbosity = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("mlp_debug");
  if (flagValue == NIL(char)) {
    option->mlpDebug = 0; /* the default value */
  }
  else {
    option->mlpDebug = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_read_reorder_cluster");
  if (flagValue == NIL(char)) {
    option->readReorderCluster = 0; /* the default value */
  }
  else {
    option->readReorderCluster = atoi(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_read_cluster");
  if (flagValue == NIL(char)) {
    option->readCluster = NIL(char); /* the default value */
  }
  else {
    option->readCluster = util_strsav(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_write_cluster");
  if (flagValue == NIL(char)) {
    option->writeCluster = NIL(char); /* the default value */
  }
  else {
    option->writeCluster = util_strsav(flagValue);
  }

  flagValue = Cmd_FlagReadByName("image_write_depend_matrix");
  if (flagValue == NIL(char)) {
    option->writeDepMatrix = NIL(char); /* the default value */
  }
  else {
    option->writeDepMatrix = util_strsav(flagValue);
  }

  flagValue = Cmd_FlagReadByName("linear_grain_type");
  if (flagValue == NIL(char)) {
    option->linearFineGrainFlag = 0; /* the default value */
  }
  else {
    option->linearFineGrainFlag = 1;
  }

  flagValue = Cmd_FlagReadByName("linear_optimize");
  if (flagValue == NIL(char)) {
    option->linearOptimize = Opt_None; /* the default value */
  }
  else {
    option->linearOptimize = Opt_NS;
  }

  option->linearOrderVariable = 0;
  flagValue = Cmd_FlagReadByName("linear_order_variable");
  if (flagValue == NIL(char)) {
    option->linearOrderVariable = 0; /* the default value */
  }
  else {
    option->linearOrderVariable = 1;
  }

  option->linearGroupStateVariable = 0; 
  flagValue = Cmd_FlagReadByName("linear_group_state_variable");
  if (flagValue == NIL(char)) {
    option->linearGroupStateVariable = 0; /* the default value */
  }
  else {
    option->linearGroupStateVariable = 1;
  }
  
  option->linearIncludeCS = 1;
  option->linearIncludeNS = 1;
  option->linearQuantifyCS = 0; 
  option->linearComputeRange = 0;
  flagValue = Cmd_FlagReadByName("linear_exclude_cs");
  if (flagValue) option->linearIncludeCS = 0;
  flagValue = Cmd_FlagReadByName("linear_exclude_ns");
  if (flagValue) option->linearIncludeNS = 0;

  return option;
}

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

  Synopsis [Creates the bit relations and the mdd id array for the
  quantifiable variables and stores it in the ImgIwls95
  structure. ]

  Description [Creates the bit relations and the mdd id array for the
  quantifiable variables and stores it in the ImgIwls95
  structure. This is identical to the initialize procedure. Read the
  roots (next state mdds), compute the bit relations for each root. If
  there are some intermediate variables, add them to quantifiable
  vars. ]
  
  SideEffects [bitRelationArray and quantifyVarMddIdArray fields are
  modified of the Iwls95Info_t struct.]
  
  SeeAlso [ImgImageInfoInitializeIwls95]
  
******************************************************************************/

static void
CreateBitRelationArray(Iwls95Info_t *info, ImgFunctionData_t *functionData)
{
  array_t *bddRelationArray = NIL(array_t);
  array_t *domainVarMddIdArray = functionData->domainVars;
  array_t *rangeVarMddIdArray = functionData->rangeVars;
  array_t *quantifyVarMddIdArray = NIL(array_t);
  graph_t *mddNetwork = functionData->mddNetwork;
  array_t *roots = functionData->roots;
  int i;
  int n1, n2, nIntermediateVars = 0;
  mdd_manager *mddManager = NIL(mdd_manager);
  st_table *vertexTable   = NIL(st_table);
  st_table *domainQuantifyVarMddIdTable = NIL(st_table);
  int mddId;
  char *nodeName;

  if (info->bitRelationArray != NIL(array_t)) {
    assert(info->quantifyVarMddIdArray != NIL(array_t));
    return;
  }
  
  bddRelationArray = array_alloc(mdd_t*, 0);
  quantifyVarMddIdArray = array_dup(functionData->quantifyVars);
  mddManager = Part_PartitionReadMddManager(mddNetwork);
  vertexTable = st_init_table(st_ptrcmp, st_ptrhash);
  domainQuantifyVarMddIdTable = st_init_table(st_numcmp, st_numhash);
  
  arrayForEachItem(int, domainVarMddIdArray, i, mddId) {
    st_insert(domainQuantifyVarMddIdTable, (char *)(long)mddId, NIL(char));
  }
  arrayForEachItem(int, quantifyVarMddIdArray, i, mddId) {
    st_insert(domainQuantifyVarMddIdTable, (char *)(long)mddId, NIL(char));
  }
  
  arrayForEachItem(char *, roots, i, nodeName) {
    vertex_t *vertex = Part_PartitionFindVertexByName(mddNetwork, nodeName);
    Mvf_Function_t *mvf = Part_VertexReadFunction(vertex);
    int mddId = array_fetch(int, rangeVarMddIdArray, i);
    
    /*mdd_t *relation = Mvf_FunctionBuildRelationWithVariable(mvf, mddId);*/
    array_t *varBddRelationArray = mdd_fn_array_to_bdd_rel_array(mddManager,
                                                                 mddId, mvf);
    array_append(bddRelationArray, varBddRelationArray);
    array_free(varBddRelationArray);
    /*array_insert_last(mdd_t *, bddRelationArray, relation);*/
    if (info->option->verbosity)
      n1 = array_n(bddRelationArray);
    else /* to remove uninitialized variable warning */
      n1 = 0;
    PartitionTraverseRecursively(mddManager, vertex, -1, vertexTable,
                                 bddRelationArray,
                                 domainQuantifyVarMddIdTable,
                                 quantifyVarMddIdArray);
    if (info->option->verbosity) {
      n2 = array_n(bddRelationArray);
      nIntermediateVars += n2 - n1;
    }
  }
  if (info->option->verbosity) {
    (void)fprintf(vis_stdout, "** img info: %d intermediate variables\n",
                  nIntermediateVars);
  }
  st_free_table(vertexTable);
  st_free_table(domainQuantifyVarMddIdTable);
  info->bitRelationArray = bddRelationArray;
  info->quantifyVarMddIdArray = quantifyVarMddIdArray;
  return;

}

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

  Synopsis    [Frees bit relation array and the quantification variables.]


  Description [Frees bit relation array and quantification
  variables. Important: to set these field to NIL in case of reuse.]

  SideEffects [Also frees the quantifVarMddIdArray field.]
******************************************************************************/
static void
FreeBitRelationArray (Iwls95Info_t *info)
{
  mdd_array_free(info->bitRelationArray);
  info->bitRelationArray = NIL(array_t);
  return;
}



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

  Synopsis [Takes an array of relations and does order, cluster,
  order on it according to the IWLS95 method. This forms a
  heuristically optimized clustered transition relation, plus a
  quantification schedule that goes with it.]

  Description [Takes an array of relations and does order,
  cluster, order on it. This forms a heuristically optimized clustered
  transition relation, plus a quantification schedule that goes with
  it. Requires the relation array, quantify variables, a
  "from" variable array (contained in the quantify variables), a "to"
  variable array and direction of computation (image or
  preimage). Eventually orderedClusteredRelationArray holds an
  ordered-clustered-ordered version of the relationArray and
  smoothVarBddArrayPtr holds the quantification schedule that goes
  with orderClusteredRelationArray.  If the freeRelationArray flag is
  set, the relationArray is freed.  This code is identical to
  ImgInfoInitialize.]

  SideEffects [The existing contents of
  orderedClusteredRelationArrayPtr or smoothVarBddArrayPtr are lost
  (not freed). relationArray is freed if the freeRelationArray flag
  is set.]

  SeeAlso     [ImgImageInfoInitializeIwls95]

******************************************************************************/
static void
OrderClusterOrder(mdd_manager *mddManager,                            
                  array_t     *relationArray, 
                  array_t     *fromVarBddArray, 
                  array_t     *toVarBddArray, 
                  array_t     *quantifyVarBddArray,
                  array_t    **orderedClusteredRelationArrayPtr,
                  array_t    **smoothVarBddArrayPtr,
                  ImgTrmOption_t   *option,
                  boolean freeRelationArray)
{
  array_t *clusteredRelationArray; /*clustered version of orderedRelationArray*/

  /*
   * Since clusters are formed by multiplying the latches starting
   * from one end we need to order the latches using some heuristics
   * first. Right now, we order the latches using the same heuristic
   * as the one used for ordering the clusters for early quantifiction.
   * However, since we are not interested in the quantification
   * schedule at this stage, we will pass a NIL(array_t*) as the last
   * argument. 
   */      
  *orderedClusteredRelationArrayPtr = NIL(array_t);
  OrderRelationArray(mddManager, relationArray, fromVarBddArray,
                     quantifyVarBddArray, toVarBddArray, option, 
                     orderedClusteredRelationArrayPtr, NIL(array_t *));

  /* free relationArray if told to */
  if (freeRelationArray) mdd_array_free(relationArray);
  
  
  /* Create the clusters */
  clusteredRelationArray = CreateClusters(mddManager,
                                          *orderedClusteredRelationArrayPtr,
                                          option);

  /* Free the orderedRelationArray */
  mdd_array_free(*orderedClusteredRelationArrayPtr);
  
  /* Order the clusters for image and pre-image computation. */
  OrderRelationArray(mddManager, clusteredRelationArray, fromVarBddArray, 
                     quantifyVarBddArray, toVarBddArray, option,
                     orderedClusteredRelationArrayPtr, smoothVarBddArrayPtr);
  /* Free the clusteredRelationArray. */
  mdd_array_free(clusteredRelationArray);
}


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

  Synopsis    [Forms the clusters of relations based on BDD size heuristic.]

  Description [The clusters are formed by taking the product in order. Once the
               BDD size of the cluster reaches a threshold, a new cluster is
               started.]

  SideEffects []
******************************************************************************/

static array_t *
CreateClusters(bdd_manager *bddManager, array_t *relationArray,
               ImgTrmOption_t *option)
{
  array_t *clusterArray;
  int i;
  bdd_t *cluster, *relation, *tempCluster;
  int flag;
  int size;

  clusterArray = array_alloc(bdd_t *, 0);

  for (i=0; i<array_n(relationArray);) {
    cluster = bdd_one(bddManager);
    flag = 0;
    do{
      relation = array_fetch(bdd_t *,relationArray, i);
      i++;
      tempCluster = bdd_and_with_limit(cluster, relation, 1, 1,
                                       option->clusterSize);
      assert(flag || tempCluster != NIL(mdd_t));
      if (tempCluster != NIL(mdd_t)) {
        size = bdd_size(tempCluster);
        if (size <= option->clusterSize || flag == 0) {
          bdd_free(cluster);
          cluster = tempCluster;
          if (flag == 0 && size >= option->clusterSize)
            break;
          flag = 1;
        }
        else{
          bdd_free(tempCluster);
          i--;
          break;
        }
      }
      else {
        i--;
        break;
      }
    } while (i < array_n(relationArray));
    array_insert_last(bdd_t *, clusterArray, cluster);
  }
  return clusterArray;
}

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

  Synopsis    [Returns an array of ordered relations and an array
               of BDD cubes (array of BDDs).]

  Description [This function returns an array of ordered relations and an array
               of BDD cubes (array of BDDs).
               This consists of following steps:
               a. Initialize the array of ctrInfoStructs and varInfoStructs.
               b. Fill in the list of varItemStruct's of ctrInfo's and
                  ctrItemStruct's of varInfo's.
               c. Simplify the relations by quantifying out the quantify
                  variables local to a particular relation.
               d. Order the relations according to the cost function described
                  in "CalculateRelationBenefit".]
  SideEffects []

******************************************************************************/
static void
OrderRelationArray(mdd_manager *mddManager,
                   array_t *relationArray,
                   array_t *domainVarBddArray,
                   array_t *quantifyVarBddArray,
                   array_t *rangeVarBddArray,
                   ImgTrmOption_t *option,
                   array_t **orderedRelationArrayPtr,
                   array_t **arraySmoothVarBddArrayPtr)
{
  array_t *quantifyVarBddIdArray, *domainVarBddIdArray, *rangeVarBddIdArray;
  array_t *domainAndQuantifyVarBddArray;
  array_t *ctrInfoArray, *varInfoArray, *simplifiedRelationArray;
  array_t *sortedMaxIndexArray;
  array_t *arrayDomainQuantifyVarsWithZeroNumCtr = NIL(array_t);
  int numRelation, numDomainVars, numQuantifyVars, numRangeVars;
  int numSmoothVars, numVars;
  int bddId;
  int i;
  int *sortedMaxIndexVector;
  int numSmoothVarsRemaining, numIntroducedVarsRemaining;
  st_table  *bddIdToVarInfoTable, *bddIdToBddTable;
  bdd_t *relation, *bdd;
  CtrInfo_t *ctrInfo;
  VarInfo_t *varInfo;
  lsList remainingCtrInfoList;

  numRelation = array_n(relationArray);
  numDomainVars = array_n(domainVarBddArray);
  numRangeVars = array_n(rangeVarBddArray);
  numQuantifyVars = array_n(quantifyVarBddArray);
  numSmoothVars = numDomainVars + numQuantifyVars;
  numVars = numSmoothVars + numRangeVars;

  domainVarBddIdArray = bdd_get_varids(domainVarBddArray);
  rangeVarBddIdArray = bdd_get_varids(rangeVarBddArray);
  quantifyVarBddIdArray = bdd_get_varids(quantifyVarBddArray);
  domainAndQuantifyVarBddArray = array_join(domainVarBddArray,
                                            quantifyVarBddArray);

  bddIdToBddTable = st_init_table(st_numcmp, st_numhash);
  HashIdToBddTable(bddIdToBddTable, domainVarBddIdArray, domainVarBddArray);
  HashIdToBddTable(bddIdToBddTable, quantifyVarBddIdArray, quantifyVarBddArray);
  HashIdToBddTable(bddIdToBddTable, rangeVarBddIdArray, rangeVarBddArray);

  bddIdToVarInfoTable = st_init_table(st_numcmp, st_numhash);

  /*
   * Create the array of ctrInfo's for each component
   */
  ctrInfoArray = array_alloc(CtrInfo_t*, 0);
  varInfoArray = array_alloc(VarInfo_t*, 0);

  /*
   * Create the array of varInfo for each variable
   */
  for (i=0; i<numVars; i++) {
    varInfo = VarInfoStructAlloc();
    if (i<numDomainVars) {
      bddId = array_fetch(int, domainVarBddIdArray, i);
      varInfo->varType = domainVar_c;
    }
    else if (i < (numDomainVars+numQuantifyVars)) {
      bddId = array_fetch(int, quantifyVarBddIdArray, i-numDomainVars);
      varInfo->varType = quantifyVar_c;
    }
    else{
      bddId = array_fetch(int, rangeVarBddIdArray,
                          (i-(numDomainVars+numQuantifyVars)));
      varInfo->varType = rangeVar_c;
    }
    array_insert_last(VarInfo_t*, varInfoArray, varInfo);
    varInfo->bddId = bddId;
    st_insert(bddIdToVarInfoTable, (char *)(long) bddId, (char *)varInfo);
  }

  /*
   * Fill in the data structure of the ctrInfo and varInfo
   */
  for (i=0; i<numRelation; i++) {
    st_table *supportTable;
    st_generator *stGen;
    VarItem_t *varItem;
    CtrItem_t *ctrItem;
    lsHandle varItemHandle, ctrItemHandle;
    long varId;

    ctrInfo = CtrInfoStructAlloc();
    array_insert_last(CtrInfo_t*, ctrInfoArray, ctrInfo);
    ctrInfo->ctrId = i;
    relation = array_fetch(bdd_t*, relationArray, i);
    supportTable = ImgBddGetSupportIdTable(relation);
    st_foreach_item(supportTable, stGen, &varId, NULL) {
      varInfo = NIL(VarInfo_t);
      if (st_lookup(bddIdToVarInfoTable, (char *) varId, &varInfo) == 0) {
        /* This variable is of no interest, because it is not present
         * in the bddIdToVarInfoTable.
         */
        continue;
      }
      varItem = ALLOC(VarItem_t,1);
      varItem->varInfo = varInfo;
      (void) lsNewBegin(ctrInfo->varItemList, (lsGeneric)varItem,
                        (lsHandle *)&varItemHandle);
      ctrItem = ALLOC(CtrItem_t,1);
      ctrItem->ctrInfo = ctrInfo;
      (void) lsNewBegin(varInfo->ctrItemList, (lsGeneric)ctrItem,
                        (lsHandle *)&ctrItemHandle);
      varItem->ctrItemHandle = ctrItemHandle;
      ctrItem->varItemHandle = varItemHandle;
      varInfo->numCtr++;
    }
    st_free_table(supportTable);
  } /* Initialization of ctrInfoArray and varInfoArray complete */

  /*
   * Smooth out the quantify variables which are local to a particular cluster.
   */
  simplifiedRelationArray = RelationArraySmoothLocalVars(relationArray,
                                                         ctrInfoArray,
                                                         varInfoArray,
                                                         bddIdToBddTable);

  assert(CheckCtrInfoArray(ctrInfoArray, numDomainVars, numQuantifyVars,
                           numRangeVars));
  assert(CheckVarInfoArray(varInfoArray, numRelation));

  /*
   * In the ordering scheme of clusters, the number of variables yet to be
   * quantified out and number of variables which are yet to be introduced is
   * taken into account.  The number of smooth variables which are yet to be
   * quantified out could have changed. Also some of the range variables may
   * not be in the support of any of the clusters. This can happen while doing
   * the ordering for clusters for backward image (when a present state
   * variable does not appear in the support of any of the next state
   * functions. Recalculate these numbers.  */

  numSmoothVarsRemaining = numSmoothVars;
  numIntroducedVarsRemaining = numRangeVars;

  /*
   * Find out which variables do not appear in the support of any cluster.
   * Update the numSmoothVarsRemaining and numIntroducedVarsRemaining
   * accordingly.
   * Also, these domainVars and quantifyVars can be quantified out as soon as
   * possible.
   */
  if (arraySmoothVarBddArrayPtr != NIL(array_t *))
    arrayDomainQuantifyVarsWithZeroNumCtr = array_alloc(bdd_t*, 0);
  for (i=0; i<numVars; i++) {
    varInfo = array_fetch(VarInfo_t*, varInfoArray, i);
    if (varInfo->numCtr == 0) {
      if ((varInfo->varType == domainVar_c) ||
          (varInfo->varType == quantifyVar_c)) {
        numSmoothVarsRemaining--;
        if (arraySmoothVarBddArrayPtr != NIL(array_t *)) {
          st_lookup(bddIdToBddTable, (char *)(long)varInfo->bddId, &bdd);
          array_insert_last(bdd_t*, arrayDomainQuantifyVarsWithZeroNumCtr,
                            bdd_dup(bdd));
        }
      }
      else numIntroducedVarsRemaining--;
    }
  }

  /*
   * The ordering scheme also depends on the value of the maximum index of a
   * smooth variable which is yet to be quantified. For efficiency reasons, we
   * maintain a vector indicating the maximum index of the clusters. This
   * vector is sorted in the decreasing order of index. The rank of a cluster
   * is stored in its "rankMaxSmoothVarIndex" field.
   */
  sortedMaxIndexArray = array_dup(ctrInfoArray);
  array_sort(sortedMaxIndexArray, CtrInfoMaxIndexCompare);
  sortedMaxIndexVector = ALLOC(int, numRelation);
  for (i=0; i<numRelation; i++) {
    ctrInfo = array_fetch(CtrInfo_t*, sortedMaxIndexArray, i);
    ctrInfo->rankMaxSmoothVarIndex = i;
    sortedMaxIndexVector[i] = ctrInfo->maxSmoothVarIndex;
  }
  array_free(sortedMaxIndexArray);

  /*
   * Create a list of clusters which are yet to be ordered. Right now
   * put all the clusters. Later on the list will contain only those
   * clusters which have not yet been ordered.
   */
  remainingCtrInfoList = lsCreate();
  for(i=0; i<numRelation; i++) {
    lsHandle ctrHandle;
    ctrInfo = array_fetch(CtrInfo_t*, ctrInfoArray, i);
    lsNewBegin(remainingCtrInfoList, (lsGeneric)ctrInfo, &ctrHandle);
    ctrInfo->ctrInfoListHandle = ctrHandle;
  }

  /* Call the auxiliary routine to do the ordering. */
  OrderRelationArrayAux(simplifiedRelationArray, remainingCtrInfoList,
                        ctrInfoArray, varInfoArray, sortedMaxIndexVector,
                        numSmoothVarsRemaining, numIntroducedVarsRemaining,
                        bddIdToBddTable, option, domainAndQuantifyVarBddArray,
                        orderedRelationArrayPtr, arraySmoothVarBddArrayPtr,
                        arrayDomainQuantifyVarsWithZeroNumCtr);

  lsDestroy(remainingCtrInfoList,0);

  /* Free the info arrays */

  for (i=0; i<numRelation; i++) {
    ctrInfo = array_fetch(CtrInfo_t*, ctrInfoArray, i);
    CtrInfoStructFree(ctrInfo);
  }
  array_free(ctrInfoArray);

  for (i=0; i<numVars; i++) {
    varInfo = array_fetch(VarInfo_t*, varInfoArray, i);
    assert(varInfo->numCtr == 0);
    VarInfoStructFree(varInfo);
  }
  array_free(varInfoArray);

  if (option->verbosity >= 3) {
    int numRelation = array_n(*orderedRelationArrayPtr);

    if (arraySmoothVarBddArrayPtr != NIL(array_t*)) {
      PrintSmoothIntroducedCount(*orderedRelationArrayPtr,
                                 arraySmoothVarBddArrayPtr,
                                 domainVarBddIdArray,
                                 rangeVarBddIdArray);
    }
    if (option->verbosity >= 4) {
      for (i= 0; i <numRelation; i++) {
        st_table *supportTable;
        (void) fprintf(vis_stdout, "Cluster # %d\n",i);
        supportTable = ImgBddGetSupportIdTable(
                       array_fetch(bdd_t*, *orderedRelationArrayPtr, i));
        PrintBddIdTable(supportTable);
        if (arraySmoothVarBddArrayPtr != NIL(array_t*)) {
          (void) fprintf(vis_stdout, "Exist cube # %d\n",i);
          PrintBddIdFromBddArray(array_fetch(array_t*,
                                             *arraySmoothVarBddArrayPtr, i));
        }
      }
    }
  }
  FREE(sortedMaxIndexVector);
  array_free(simplifiedRelationArray);

  /* Free BDD Id arrays */
  array_free(domainVarBddIdArray);
  array_free(quantifyVarBddIdArray);
  array_free(rangeVarBddIdArray);
  array_free(domainAndQuantifyVarBddArray);
  /* Free the st_tables */
  st_free_table(bddIdToBddTable);
  st_free_table(bddIdToVarInfoTable);
}

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

  Synopsis    [This function takes an array of relations and quantifies out the
               quantify variables which are local to a particular relation.]

  Description [This function takes an array of relations and quantifies out the
               quantify variables which are local to a particular relation.]

  SideEffects [This function fills in the "numSmoothVars",
               "numLocalSmoothVars", "numIntroducedVars", "maxSmoothVarIndex"
               fields of ctrInfoStruct corresponding to each relation. It also
               alters the "numCtr" and "ctrItemList" field of those quantify
               variables which are quantified out in this function.]

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

static array_t *
RelationArraySmoothLocalVars(array_t *relationArray, array_t *ctrInfoArray,
                             array_t *varInfoArray,
                             st_table *bddIdToBddTable)
{
  array_t *arraySmoothVarBddArray, *smoothVarBddArray;
  array_t *simplifiedRelationArray;
  bdd_t *simplifiedRelation, *relation, *varBdd;
  int maxSmoothVarIndexForAllCtr, i, numRelation;

  numRelation = array_n(relationArray);

  /*
   * Initialize the array of cubes (of quantified variables which occur in
   * only one relation), to be smoothed out.
   */
  arraySmoothVarBddArray = array_alloc(array_t*, 0);
  for (i=0; i<numRelation; i++) {
    smoothVarBddArray = array_alloc(bdd_t*, 0);
    array_insert_last(array_t*, arraySmoothVarBddArray, smoothVarBddArray);
  }

  maxSmoothVarIndexForAllCtr = -1;
  for (i=0; i<numRelation; i++) {
    VarItem_t *varItem;
    lsHandle varItemHandle;
    CtrInfo_t *ctrInfo = array_fetch(CtrInfo_t*, ctrInfoArray, i);
    int maxSmoothVarIndex = -1;
    lsGen varItemListGen = lsStart(ctrInfo->varItemList);

    smoothVarBddArray = array_fetch(array_t *, arraySmoothVarBddArray, i);
    while (lsNext(varItemListGen, &varItem, &varItemHandle) ==
           LS_OK) {
      int  varBddId;
      VarInfo_t *varInfo = varItem->varInfo;
      if (varInfo->varType == rangeVar_c) {
        ctrInfo->numIntroducedVars++;
      }
      else if (varInfo->numCtr == 1) {
        assert(lsLength(varInfo->ctrItemList) == 1);
        if (varInfo->varType == quantifyVar_c) {
          /*
           * The variable can be smoothed out.
           * Put it in an array of variables to be
           * smoothed out from the relation.
           */
          CtrItem_t *ctrItem;
          varBddId = varInfo->bddId;
          st_lookup(bddIdToBddTable, (char *)(long)varBddId,&varBdd);
          array_insert_last(bdd_t*, smoothVarBddArray, bdd_dup(varBdd));
          varInfo->numCtr = 0;
          /* Remove the cluster from the ctrItemList of varInfo */
          lsRemoveItem(varItem->ctrItemHandle, &ctrItem);
          CtrItemStructFree(ctrItem);
          /* Remove the variable from the varItemList of ctrInfo.*/
          lsRemoveItem(varItemHandle, &varItem);
          VarItemStructFree(varItem);
          continue;
        }
        else {
          /* Increase the numLocalSmoothVars count of the corresponding ctr. */
          ctrInfo->numLocalSmoothVars++;
          ctrInfo->numSmoothVars++;
          if (maxSmoothVarIndex < varInfo->bddId) {
            maxSmoothVarIndex = varInfo->bddId;
          }
        }
      }
      else{ /* varInfo->numCtr > 1 */
        assert(varInfo->numCtr > 1);
        ctrInfo->numSmoothVars++;
        if (maxSmoothVarIndex < varInfo->bddId) {
          maxSmoothVarIndex = varInfo->bddId;
        }
      }
    }
    lsFinish(varItemListGen);
    if (maxSmoothVarIndex >= 0) {
      ctrInfo->maxSmoothVarIndex = maxSmoothVarIndex;
    }
    else{
      ctrInfo->maxSmoothVarIndex = 0;
    }
    if (maxSmoothVarIndexForAllCtr < maxSmoothVarIndex) {
      maxSmoothVarIndexForAllCtr = maxSmoothVarIndex;
    }
  }

  /*
   * Initialization Finished. We need to smooth out those quantify
   * variables which appear in only one ctr.
   */
  simplifiedRelationArray = array_alloc(bdd_t*, 0);
  for (i=0; i<numRelation; i++) {
    relation = array_fetch(bdd_t*, relationArray, i);
    smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, i);
    if (array_n(smoothVarBddArray) != 0)
      simplifiedRelation = bdd_smooth(relation, smoothVarBddArray);
    else
      simplifiedRelation = bdd_dup(relation);
    array_insert_last(bdd_t*, simplifiedRelationArray, simplifiedRelation);
  }
  mdd_array_array_free(arraySmoothVarBddArray);
  return simplifiedRelationArray;
}

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

  Synopsis    [An auxiliary routine for orderRelationArray to order the
  clusters.]

  Description [This routine is called by orderRelationArray after the
  initialization is complete and the cluster relations are simplified wrt to
  the quantify variable local to the cluster. The algorithm is:
      While (there exists a cluster to be ordered) {
         Calculate benefit for each unordered cluster.
         Choose the cluster with the maximum benefit.
         Update the cost function parameters.
      }
  ]

  SideEffects [ctrInfo and varInfo structures are modified.]

******************************************************************************/
static void
OrderRelationArrayAux(array_t *relationArray,
                      lsList remainingCtrInfoList,
                      array_t *ctrInfoArray,
                      array_t *varInfoArray,
                      int *sortedMaxIndexVector,
                      int numSmoothVarsRemaining,
                      int numIntroducedVarsRemaining,
                      st_table *bddIdToBddTable,
                      ImgTrmOption_t *option,
                      array_t *domainAndQuantifyVarBddArray,
                      array_t **orderedRelationArrayPtr,
                      array_t **arraySmoothVarBddArrayPtr,
                      array_t *arrayDomainQuantifyVarsWithZeroNumCtr)
{
  int numRelation, ctrCount, currentRankMaxSmoothVarIndex;
  int maxIndex = 0;
  array_t *orderedRelationArray, *arraySmoothVarBddArray;
  array_t *smoothVarBddArray;
  int *bestCtrIdVector;
  int maxNumLocalSmoothVars, maxNumSmoothVars, maxNumIntroducedVars;
  lsGen lsgen;

  ctrCount = 0;
  currentRankMaxSmoothVarIndex = 0;

  numRelation = array_n(relationArray);

  arraySmoothVarBddArray = NIL(array_t);
  orderedRelationArray = array_alloc(bdd_t*, 0);
  *orderedRelationArrayPtr = orderedRelationArray;
  if (arraySmoothVarBddArrayPtr != NIL(array_t *)) {
    arraySmoothVarBddArray = array_alloc(array_t*, 0);
    *arraySmoothVarBddArrayPtr = arraySmoothVarBddArray;
    array_insert_last(array_t*, arraySmoothVarBddArray,
                        arrayDomainQuantifyVarsWithZeroNumCtr);
  }
  bestCtrIdVector = ALLOC(int, numRelation);

  while (ctrCount < numRelation) {
    float bestBenefit, benefit;
    int bestCtrId;
    bdd_t *bestRelation;
    lsGen ctrInfoListGen;
    CtrInfo_t *ctrInfo, *ctrInfoAux;

    bestBenefit = -1.0e20;      /* a safely small number */
    bestCtrId = -1;
    /* Find the maximum index of the remaining clusters */
    while (currentRankMaxSmoothVarIndex < numRelation &&
      (maxIndex = sortedMaxIndexVector[currentRankMaxSmoothVarIndex++]) == -1);

    /* Calculate the maximum values of local smooth variables etc. */
    maxNumLocalSmoothVars = 0;
    maxNumSmoothVars = 0;
    maxNumIntroducedVars = 0;
    lsForEachItem(remainingCtrInfoList, lsgen, ctrInfo) {
      if (ctrInfo->numLocalSmoothVars > maxNumLocalSmoothVars) {
        maxNumLocalSmoothVars = ctrInfo->numLocalSmoothVars;
      }
      if (ctrInfo->numSmoothVars > maxNumSmoothVars) {
        maxNumSmoothVars = ctrInfo->numSmoothVars;
      }
      if (ctrInfo->numIntroducedVars > maxNumIntroducedVars) {
        maxNumIntroducedVars = ctrInfo->numIntroducedVars;
      }
    }

    if (option->verbosity >= 4) {
      fprintf(vis_stdout, "maxNumLocalSmoothVars = %3d maxNumSmoothVars = %3d ",
              maxNumLocalSmoothVars, maxNumSmoothVars);
      fprintf(vis_stdout, "maxNumIntroducedVars = %3d\n", maxNumIntroducedVars);
    }
    /* Calculate the cost function of each of the cluster */
    ctrInfoListGen = lsStart(remainingCtrInfoList);
    while (lsNext(ctrInfoListGen, &ctrInfo, NIL(lsHandle)) ==
           LS_OK) {
      benefit = CalculateBenefit(ctrInfo, maxNumLocalSmoothVars,
                                 maxNumSmoothVars, maxIndex,
                                 maxNumIntroducedVars, option);

      if (option->verbosity >= 4) {
        fprintf(vis_stdout,
          "id = %d: numLocalSmoothVars = %d numSmoothVars = %d benefit = %f\n",
                ctrInfo->ctrId, ctrInfo->numLocalSmoothVars,
                ctrInfo->numSmoothVars, benefit);
      }
      if (benefit > bestBenefit) {
        bestBenefit = benefit;
        bestCtrId = ctrInfo->ctrId;
      }
    }
    lsFinish(ctrInfoListGen);
    /*
     * Insert the relation in the ordered array of relations and put in the
     * appropriate cubes.
     */

    bestCtrIdVector[ctrCount] = bestCtrId;
    ctrInfo = array_fetch(CtrInfo_t*, ctrInfoArray, bestCtrId);
    lsRemoveItem(ctrInfo->ctrInfoListHandle, &ctrInfoAux);
    assert(ctrInfo == ctrInfoAux);
    bestRelation = array_fetch(bdd_t*, relationArray, bestCtrId);
    array_insert_last(bdd_t*, orderedRelationArray, bestRelation);
    if (option->verbosity >= 4) {
      fprintf(vis_stdout,
       "Best id = %d benefit = %f numLocalSmoothVars = %d numSmoothVars = %d\n",
              bestCtrId, bestBenefit, ctrInfo->numLocalSmoothVars,
              ctrInfo->numSmoothVars);
    }
    /*
     * Update the remaining ctrInfo's and the varInfo's affected by choosing
     * this cluster. Also get the array of smooth variables which can be
     * quantified once this cluster is multiplied in the product.
     */
    smoothVarBddArray = UpdateInfoArrays(ctrInfo, bddIdToBddTable,
                                         &numSmoothVarsRemaining,
                                         &numIntroducedVarsRemaining);
    assert(CheckCtrInfo(ctrInfo, numSmoothVarsRemaining, 0,
                        numIntroducedVarsRemaining));

    if (arraySmoothVarBddArrayPtr != NIL(array_t *)) {
      if (ctrCount == numRelation-1) {
        /*
         * In the last element of arraySmoothVarBddArray, put all the domain
         * and quantify variables (in case some of them were not in the support
         * of any cluster).
         */
        int             i, j;
        st_table        *varsTable = st_init_table((ST_PFICPCP)bdd_ptrcmp,
                                                   (ST_PFICPI)bdd_ptrhash);
        array_t         *tmpVarBddArray;
        mdd_t           *tmpVar;

        arrayForEachItem(array_t *, arraySmoothVarBddArray, i, tmpVarBddArray) {
          arrayForEachItem(mdd_t *, tmpVarBddArray, j, tmpVar) {
            st_insert(varsTable, (char *)tmpVar, NIL(char));
          }
        }
        arrayForEachItem(mdd_t *, smoothVarBddArray, i, tmpVar) {
          st_insert(varsTable, (char *)tmpVar, NIL(char));
        }
        arrayForEachItem(mdd_t *, domainAndQuantifyVarBddArray, i, tmpVar) {
          if (st_lookup(varsTable, (char *)tmpVar, NIL(char *)) == 0)
            array_insert_last(mdd_t *, smoothVarBddArray, mdd_dup(tmpVar));
        }
        st_free_table(varsTable);
      }
      array_insert_last(array_t*, arraySmoothVarBddArray, smoothVarBddArray);
    }
    else{
      mdd_array_free(smoothVarBddArray);
    }
    sortedMaxIndexVector[ctrInfo->rankMaxSmoothVarIndex] = -1;
    ctrCount++;
  }
  assert(numIntroducedVarsRemaining == 0);
  assert(numSmoothVarsRemaining == 0);
  if (option->verbosity >= 3) {
    int i;
    (void) fprintf(vis_stdout,"Cluster Sequence = ");
    for (i=0; i<numRelation; i++) {
      (void) fprintf(vis_stdout, "%d ", bestCtrIdVector[i]);
    }
    (void) fprintf(vis_stdout,"\n");
  }

  FREE(bestCtrIdVector);
}


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

  Synopsis    [This function computes the set of variables which can be
  smoothed out once a particular cluster is chosen to be multiplied in the
  product.]

  Description [The support variable list (varItemList) of the cluster (ctrInfo)
  is traversed and depending upon the type of the variable and the number of
  unordered clustered relations which depend on that variable following actions
  are taken:

  a. If the variable is of range type:
     This implies that this variable is appearing for the first time in the
     product. Since it is already introduced in the product, the cost function
     of any other unordered relation which depends on this variable will get
     modified. The numIntroducedVariables field of each of the cluster which
     depends on this variable is decreased by 1. Also this varItemStruct
     corresponding to this variable is removed from the varItemList field of
     the cluster.

  b. If the variable is of domain or quantify type:
    b1. If number of clusters which depend on this variable is 1 (numCtr == 1):
           In this case, this variable can be quantified out once the chosen
           cluster is multiplied in the product. Hence the variable is put in
           the smoothVarBddArray.
    b2. If  (numCtr == 2)
           In this case, there is one more unordered cluster which depends on
           this variable. But once the current cluster is multiplied in the
           product, the "numLocalSmoothVars" field of ctrInfo corresponding to
           the other dependent cluster needs to be increased by 1.
     b3. If (numCtr > 2)
           In this case, we just need to decrease the numCtr of the variable by
           1.

  In any case, for each varInfo in the support variable list (varItemList) of
  the cluster (ctrInfo) the following invariant is maintained:
  varInfo->numCtr == lsLength(varInfo->ctrItemList)
  ]

  SideEffects [The fields of ctrInfo and varInfo are changed as mentioned
  above.]

******************************************************************************/
static array_t *
UpdateInfoArrays(CtrInfo_t *ctrInfo, st_table *bddIdToBddTable,
                 int *numSmoothVarsRemainingPtr,
                 int *numIntroducedVarsRemainingPtr)
{
  array_t *smoothVarBddArray;
  lsGen varItemListGen;
  int numSmoothVarsRemaining = *numSmoothVarsRemainingPtr;
  int numIntroducedVarsRemaining = *numIntroducedVarsRemainingPtr;
  VarItem_t *varItem;
  lsHandle varItemHandle;

  smoothVarBddArray = array_alloc(bdd_t*, 0);
  varItemListGen = lsStart(ctrInfo->varItemList);
  while (lsNext(varItemListGen, &varItem, &varItemHandle) ==
         LS_OK) {
    VarInfo_t *varInfo = varItem->varInfo;
    CtrItem_t *ctrItem;
    assert(varInfo->numCtr == lsLength(varInfo->ctrItemList));
    lsRemoveItem(varItem->ctrItemHandle,  &ctrItem);
    CtrItemStructFree(ctrItem);
    varInfo->numCtr--;
    lsRemoveItem(varItemHandle, &varItem);
    VarItemStructFree(varItem);

    /*
     * If this variable is to be smoothed (domain or quantify type) and
     * the numCtr is 1, then it should be added to the smoothVarBddArray,
     * otherwise, the numCtr should be decreased by 1.
     * If the variable is of the range type then the number of introduced
     * vars remaining and the number of introduced vars should be
     * appropriately modified.
     */
    if ((varInfo->varType == domainVar_c) ||
        (varInfo->varType == quantifyVar_c)) {
      if (varInfo->numCtr == 0) {
        bdd_t *varBdd;
        st_lookup(bddIdToBddTable, (char *)(long)(varInfo->bddId), &varBdd);
        array_insert_last(bdd_t*, smoothVarBddArray, bdd_dup(varBdd));
        numSmoothVarsRemaining--;
        ctrInfo->numLocalSmoothVars--;
        ctrInfo->numSmoothVars--;
      }
      else{
        if (varInfo->numCtr == 1) {
          /*
           * We need to update the numLocalSmoothVars of the ctr
           * which depends on it.
           */
          lsFirstItem(varInfo->ctrItemList, &ctrItem, LS_NH);
          ctrItem->ctrInfo->numLocalSmoothVars++;
        }
        /*
         * else varInfo->numCtr > 1 : Nothing to be done because neither it
         * can be quantified out, nor it is effecting any cost function.
         */
        ctrInfo->numSmoothVars--;
      }
    }
    else{/* The variable is of range type, so need to appropriately modify the
          * numIntroducedVars of those clusters which contain this range
          * variable in their support.
          */
      lsHandle ctrItemHandle;
      lsGen ctrItemListGen = lsStart(varInfo->ctrItemList);
      ctrInfo->numIntroducedVars--;
      while (lsNext(ctrItemListGen, &ctrItem, &ctrItemHandle) == LS_OK) {
        ctrItem->ctrInfo->numIntroducedVars--;
        lsRemoveItem(ctrItem->varItemHandle,&varItem);
        lsRemoveItem(ctrItemHandle,&ctrItem);
        VarItemStructFree(varItem);
        CtrItemStructFree(ctrItem);
      }
      lsFinish(ctrItemListGen);
      numIntroducedVarsRemaining--;
      varInfo->numCtr = 0;
    }
    assert(varInfo->numCtr == lsLength(varInfo->ctrItemList));
  }
  lsFinish(varItemListGen);
  *numIntroducedVarsRemainingPtr = numIntroducedVarsRemaining;
  *numSmoothVarsRemainingPtr = numSmoothVarsRemaining;
  return smoothVarBddArray;
}


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

  Synopsis    [Gets the necessary options for computing the image and returns
               in the option structure.]

  Description [
       The strategy
       Choose the cost function:
       The objective function attached with each Ti is
       Ci =  W1 C1i + W2 C2i + W3 C3i - W4C4i
       where:
       W1 = weight attached with variable getting smoothed
       W2 = weight attached with the support count of the Ti
       W3 = weight attached with the max id of the Ti
       W4 = weight attached with variable getting introduced
       C1i = Ui/Vi
       C2i = Vi/Wi
       C3i = Mi/Max
       C4i = Xi/Yi
       Ui = number of variables getting smoothed
       Vi = number of ps support variables of Ti
       Wi = total number of ps variables remaining which have not been
            smoothed out
       Mi = value of Max id of Ti
       Max = value of Max id across all the Ti's remaining to be multiplied
       Xi = number of ns variables introduced
       Yi = total number of ns variables which have not been introduced so
            far.
       Get the weights from the global option]

  SideEffects []

******************************************************************************/
static float
CalculateBenefit(CtrInfo_t *ctrInfo, int maxNumLocalSmoothVars, int
                 maxNumSmoothVars, int maxIndex, int
                 maxNumIntroducedVars, ImgTrmOption_t *option)
{
  int W1, W2, W3, W4;
  float benefit;

  W1 = option->W1;
  W2 = option->W2;
  W3 = option->W3;
  W4 = option->W4;

  benefit = 0.0;
  benefit += (maxNumLocalSmoothVars ?
              ((float)W1 * ((float)ctrInfo->numLocalSmoothVars /
                (float)maxNumLocalSmoothVars)) : 0.0);

  benefit += (maxNumSmoothVars ?
              ((float)W2 * ((float)ctrInfo->numSmoothVars /
                (float)maxNumSmoothVars)) : 0.0);

  benefit += (maxIndex ?
              ((float)W3 * ((float)ctrInfo->maxSmoothVarIndex /
                (float)maxIndex)) : 0.0);

  benefit -= (maxNumIntroducedVars ?
              ((float)W4 * ((float)ctrInfo->numIntroducedVars /
                (float)maxNumIntroducedVars)) : 0.0);

  return benefit;
}

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

  Synopsis    [Prints the option values used in IWLS95 techinique for
  image computation.]

  Description [Prints the option values used in IWLS95 techinique for
  image computation.]

  SideEffects []

******************************************************************************/
static void
PrintOption(Img_MethodType method, ImgTrmOption_t *option, FILE *fp)
{
  if (method == Img_Iwls95_c)
    (void)fprintf(fp, "Printing Information about Image method: IWLS95\n");
  else if (method == Img_Mlp_c)
    (void)fprintf(fp, "Printing Information about Image method: MLP\n");
  else if (method == Img_Linear_c)
    (void)fprintf(fp, "Printing Information about Image method: LINEAR\n");
  else
    assert(0);
  (void)fprintf(fp,
                "\tThreshold Value of Bdd Size For Creating Clusters = %d\n",
                option->clusterSize);
  (void)fprintf(fp,
  "\t\t(Use \"set image_cluster_size value\" to set this to desired value)\n");
  (void)fprintf(fp, "\tVerbosity = %d\n", option->verbosity);
  (void)fprintf(fp,
    "\t\t(Use \"set image_verbosity value\" to set this to desired value)\n");
  if (method == Img_Iwls95_c) {
    (void)fprintf(fp, "\tW1 =%3d W2 =%2d W3 =%2d W4 =%2d\n",
                  option->W1, option->W2, option->W3, option->W4);
    (void)fprintf(fp,
        "\t\t(Use \"set image_W? value\" to set these to desired values)\n");
  }
}

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

  Synopsis    [Allocates and initializes the info structure for IWLS95
  technique.]

  Description [Allocates and initializes the info structure for IWLS95
  technique.]

  SideEffects []

******************************************************************************/
static Iwls95Info_t *
Iwls95InfoStructAlloc(Img_MethodType method)
{
  Iwls95Info_t *info;
  info = ALLOC(Iwls95Info_t, 1);
  memset(info, 0, sizeof(Iwls95Info_t));
  info->method = method;
  info->option = TrmGetOptions();
  info->PIoption = ImgGetPartialImageOptions();
  return info;
}

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

  Synopsis    [Allocates and initializes memory for ctrInfoStruct.]

  Description [Allocates and initializes memory for ctrInfoStruct.]

  SideEffects []

******************************************************************************/
static CtrInfo_t *
CtrInfoStructAlloc(void)
{
  CtrInfo_t *ctrInfo;
  ctrInfo = ALLOC(CtrInfo_t, 1);
  ctrInfo->ctrId = -1;
  ctrInfo->numLocalSmoothVars = 0;
  ctrInfo->numSmoothVars = 0;
  ctrInfo->maxSmoothVarIndex = -1;
  ctrInfo->numIntroducedVars = 0;
  ctrInfo->varItemList = lsCreate();
  ctrInfo->ctrInfoListHandle = NULL;
  return ctrInfo;
}


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

  Synopsis    [Frees the memory associated with ctrInfoStruct.]

  Description [Frees the memory associated with ctrInfoStruct.]

  SideEffects []

******************************************************************************/
static void
CtrInfoStructFree(CtrInfo_t *ctrInfo)
{
  lsDestroy(ctrInfo->varItemList, 0);
  FREE(ctrInfo);
}

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

  Synopsis    [Allocates and initializes memory for varInfoStruct.]

  Description [Allocates and initializes memory for varInfoStruct.]

  SideEffects []

******************************************************************************/
static VarInfo_t *
VarInfoStructAlloc(void)
{
  VarInfo_t *varInfo;
  varInfo = ALLOC(VarInfo_t, 1);
  varInfo->bddId = -1;
  varInfo->numCtr = 0;
  varInfo->varType = -1;
  varInfo->ctrItemList = lsCreate();
  return varInfo;
}

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

  Synopsis    [Frees the memory associated with varInfoStruct.]

  Description [Frees the memory associated with varInfoStruct.]

  SideEffects []

******************************************************************************/
static void
VarInfoStructFree(VarInfo_t *varInfo)
{
  lsDestroy(varInfo->ctrItemList,0);
  FREE(varInfo);
}

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

  Synopsis    [Frees the memory associated with CtrItemStruct]

  Description [Frees the memory associated with CtrItemStruct]

  SideEffects []

******************************************************************************/
static void
CtrItemStructFree(CtrItem_t *ctrItem)
{
  FREE(ctrItem);
}

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

  Synopsis    [Frees the memory associated with VarItemStruct]

  Description [Frees the memory associated with VarItemStruct]

  SideEffects []

******************************************************************************/
static void
VarItemStructFree(VarItem_t *varItem)
{
  FREE(varItem);
}

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

  Synopsis    [Compare function used to sort the ctrInfoStruct based on the
  maxSmoothVarIndex field.]

  Description [This function is used to sort the array of clusters based on the
  maximum index of the support variable.]

  SideEffects []

******************************************************************************/
static int
CtrInfoMaxIndexCompare(const void * p1, const void * p2)
{
  CtrInfo_t **infoCtr1 = (CtrInfo_t **) p1;
  CtrInfo_t **infoCtr2 = (CtrInfo_t **) p2;
  return ((*infoCtr1)->maxSmoothVarIndex > (*infoCtr2)->maxSmoothVarIndex);
}

/* **************************************************************************
 * Debugging Related Routines.
 ***************************************************************************/


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

  Synopsis    [Prints the CtrInfo_t data structure.]

  Description [Prints the CtrInfo_t data structure.]

  SideEffects []

******************************************************************************/
static void
PrintCtrInfoStruct(CtrInfo_t *ctrInfo)
{
  lsGen lsgen;
  VarItem_t *varItem;

  (void) fprintf(vis_stdout,"Ctr ID = %d\tNumLocal = %d\tNumSmooth=%d\tNumIntro=%d\tmaxSmooth=%d\trank=%d\n",
                 ctrInfo->ctrId, ctrInfo->numLocalSmoothVars,
                 ctrInfo->numSmoothVars, ctrInfo->numIntroducedVars,
                 ctrInfo->maxSmoothVarIndex, ctrInfo->rankMaxSmoothVarIndex);
  lsgen = lsStart(ctrInfo->varItemList);
  while (lsNext(lsgen, &varItem, NIL(lsHandle)) == LS_OK) {
    (void) fprintf(vis_stdout,"%d\t", varItem->varInfo->bddId);
  }
  lsFinish(lsgen);
  fprintf(vis_stdout,"\n");
}
/**Function********************************************************************

  Synopsis    [Prints the VarInfo_t structure.]

  Description [Prints the VarInfo_t structure.]

  SideEffects []

******************************************************************************/
static void
PrintVarInfoStruct(VarInfo_t *varInfo)
{
  lsGen lsgen;
  CtrItem_t *ctrItem;

  (void) fprintf(vis_stdout,"Var ID = %d\tNumCtr = %d\tVarType=%d\n",
                 varInfo->bddId, varInfo->numCtr, varInfo->varType);
  lsgen = lsStart(varInfo->ctrItemList);
  while (lsNext(lsgen, &ctrItem, NIL(lsHandle)) == LS_OK) {
    (void) fprintf(vis_stdout,"%d\t", ctrItem->ctrInfo->ctrId);
  }
  lsFinish(lsgen);
  fprintf(vis_stdout,"\n");
}

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

  Synopsis    [Given an array of BDD's representing the relations and another
  array of BDD cubes (in the form of array of bdd_t of variable) representing
  the quantification schedule, this function checks if the schedule is correct.
  ]

  Description [Assume Ci represents the ith cube in the array and the Ti
  represents the ith relation. The correctness of the schedule is defined as
  follows:
      a. For all Tj: j > i, SUP(Tj) and SUP(Ci) do not intersect, i.e., the
      variables which are quantified in Ci should not appear in the Tj for j>i.
      b. For any i, j, SUP(Ci) and SUP(Cj) do not intersect. However this is
      not true for the last cube (Cn-1). This is because the last cube contains
      all the smooth variables.
  ]
  SideEffects []

******************************************************************************/
static int
CheckQuantificationSchedule(array_t *relationArray,
                            array_t *arraySmoothVarBddArray)
{
  int i, j;
  long varId;
  st_table *smoothVarTable, *supportTable;
  bdd_t *relation;
  array_t *smoothVarBddArray, *smoothVarBddIdArray;
  st_generator *stgen;

  assert(array_n(relationArray) + 1 == array_n(arraySmoothVarBddArray));

  smoothVarTable = st_init_table(st_numcmp, st_numhash);

  smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, 0);
  smoothVarBddIdArray = bdd_get_varids(smoothVarBddArray);
  for (j=0; j<array_n(smoothVarBddIdArray); j++) {
    varId = (long) array_fetch(int, smoothVarBddIdArray, j);
    assert(st_insert(smoothVarTable, (char *) varId, NIL(char))==0);
  }
  array_free(smoothVarBddIdArray);

  for (i=0; i<array_n(relationArray); i++) {
    relation = array_fetch(bdd_t*, relationArray, i);
    supportTable = ImgBddGetSupportIdTable(relation);
    /*
     * Check that none of the variables in the support have already been
     * quantified.
     */
    st_foreach_item(supportTable, stgen, &varId, NULL) {
      assert(st_is_member(smoothVarTable, (char *) varId) == 0);

    }
    st_free_table(supportTable);
    if (i == (array_n(relationArray)-1)) {
      /* Since last element of arraySmoothVarBddArray has all the
       * smooth variables, it will not satisfy the condition.
       */
      continue;
    }

    /*
     * Put each of the variables quantified at this step in the
     * smoothVarTable. And check that none of these smooth variables
     * have been introduced before.
     */
    smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, i + 1);
    smoothVarBddIdArray = bdd_get_varids(smoothVarBddArray);
    for (j=0; j<array_n(smoothVarBddIdArray); j++) {
      varId = (long) array_fetch(int, smoothVarBddIdArray, j);
      assert(st_insert(smoothVarTable, (char *) varId, NIL(char))==0);
    }
    array_free(smoothVarBddIdArray);
  }
  st_free_table(smoothVarTable);
  return 1;
}

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

  Synopsis    [Checks the validity of an array of CtrInfoStructs.]

  Description [Checks the validity of an array of CtrInfoStructs.]

  SideEffects []

******************************************************************************/
static int
CheckCtrInfoArray(array_t *ctrInfoArray, int numDomainVars, int
                  numQuantifyVars, int numRangeVars)
{
  int i;
  for (i=0; i<array_n(ctrInfoArray); i++) {
    CheckCtrInfo(array_fetch(CtrInfo_t *, ctrInfoArray, i),
                 numDomainVars, numQuantifyVars, numRangeVars);
  }
  return 1;
}

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

  Synopsis    [Checks the validity of a CtrInfoStruct.]

  Description [Checks the validity of a CtrInfoStruct.]

  SideEffects []

******************************************************************************/
static int
CheckCtrInfo(CtrInfo_t *ctrInfo, int numDomainVars,
             int numQuantifyVars, int numRangeVars)
{
  assert(ctrInfo->numLocalSmoothVars <= (numDomainVars +
                                         numQuantifyVars));
  assert(ctrInfo->numSmoothVars <= (numDomainVars +
                                    numQuantifyVars));
  assert(ctrInfo->numIntroducedVars <= numRangeVars);
  assert(lsLength(ctrInfo->varItemList) ==
         (ctrInfo->numSmoothVars+ctrInfo->numIntroducedVars));
  return 1;
}

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

  Synopsis    [Checks the validity of an array of VarInfoStruct.]

  Description [Checks the validity of an array of VarInfoStruct.]

  SideEffects []

******************************************************************************/
static int
CheckVarInfoArray(array_t *varInfoArray, int numRelation)
{
  int i;
  for (i=0; i<array_n(varInfoArray); i++) {
    VarInfo_t *varInfo;
    varInfo = array_fetch(VarInfo_t *, varInfoArray, i);
    assert(varInfo->numCtr <= numRelation);
    assert(varInfo->varType >= 0);
    assert(varInfo->varType <= 2);
    assert(lsLength(varInfo->ctrItemList) == varInfo->numCtr);
  }
  return 1;
}

/* ****************************************************************
 * Utility Routines.
 ****************************************************************/

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

  Synopsis    [Duplicates an array of BDDs.]

  Description [Duplicates an array of BDDs.]

  SideEffects []

******************************************************************************/
static array_t *
BddArrayDup(array_t *bddArray)
{
  int i;
  array_t *resultArray;
  bdd_t *bdd;
  resultArray = array_alloc(bdd_t*, 0);
  arrayForEachItem(bdd_t *, bddArray, i, bdd) {
    array_insert_last(bdd_t*, resultArray, bdd_dup(bdd));
  }
  return resultArray;
}

static array_t *
BddArrayArrayDup(array_t *bddArray)
{
  int i;
  array_t *resultArray;
  array_t *arr, *tarr;

  resultArray = array_alloc(bdd_t*, 0);
  arrayForEachItem(array_t *, bddArray, i, arr) {
    tarr = BddArrayDup(arr);
    array_insert_last(array_t*, resultArray, tarr);
  }
  return resultArray;
}


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

  Synopsis [Returns a BDD after taking the product of fromSet with the
  BDDs in the relationArray and quantifying out the variables using
  the schedule given in the arraySmoothVarBddArray.]

  Description [The product is taken from the left, i.e., fromSet is
  multiplied with relationArray[0]. The array of variables given by
  arraySmoothVarBddArray[i] are quantified when the relationArray[i]
  is multiplied in the product. We notice that no simplification is
  used in the computation. It was found to be the fastest way to
  compute the image. When partial image is allowed, a partial image is
  computed based on the method specified. If the method is "approx",
  approximation is applied to some intermediate product. If the method
  is "clipping", clipping is introduced during and-abstract. A flag
  partial is set to indicate to the calling procedure that a partial
  image was computed as against an exact image. Partial flag comes in
  indicating whether partial images are allowed, and leave indicating
  whether the image was partial or not. ]

  SideEffects [partial is set when some approximation is applied]

******************************************************************************/
static mdd_t *
BddLinearAndSmooth(mdd_manager *mddManager,
                   bdd_t *fromSet,
                   array_t *relationArray,
                   array_t *arraySmoothVarBddArray,
                   array_t *smoothVarCubeArray,
                   int verbosity,
                   ImgPartialImageOption_t *PIoption,
                   boolean *partial, boolean lazySiftFlag)
{
  int i;
  int intermediateSize = 0;
  int maxSize = 0;
  mdd_t *newProduct;
  boolean partialImageAllowed;
  bdd_t *product = fromSet;
  int clippingDepth = -1;
  int partialImageThreshold = 0;
  int partialImageMethod = Img_PIApprox_c;
  array_t *smoothVarBddArray = NIL(array_t);
  mdd_t *smoothVarCube = NIL(mdd_t);
  int nvars = -2;

  assert(CheckQuantificationSchedule(relationArray, arraySmoothVarBddArray));

   /* Copy (*partial) which indicates whether partial image is allowed. */
  /* (*partial) is set to TRUE later if partial image was performed. */
  partialImageAllowed = *partial;
  *partial = FALSE;

  /* partial image options */
  if (PIoption != NULL) {
    nvars = PIoption->nvars;
    clippingDepth = (int) (PIoption->clippingDepthFrac*nvars);
    partialImageThreshold = PIoption->partialImageThreshold;
    partialImageMethod = PIoption->partialImageMethod;
  }

  if (fromSet) {
    if (smoothVarCubeArray) {
      smoothVarCube = array_fetch(mdd_t*, smoothVarCubeArray, 0);
      if (bdd_is_tautology(smoothVarCube, 1))
        product = bdd_dup(fromSet);
      else
        product = bdd_smooth_with_cube(fromSet, smoothVarCube);
    } else {
      smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, 0);
      if (array_n(smoothVarBddArray) == 0)
        product = bdd_dup(fromSet);
      else
        product = bdd_smooth(fromSet, smoothVarBddArray);
    }
  } else
    product = bdd_one(mddManager);

  for (i=0; i<array_n(relationArray); i++) {
    bdd_t *relation, *tmpProduct;
    boolean allowClipping = FALSE;

    /* fetch relation to intersect with the next product. */
    relation = array_fetch(bdd_t*, relationArray, i);
    if (smoothVarCubeArray)
      smoothVarCube = array_fetch(mdd_t*, smoothVarCubeArray, i + 1);
    else
      smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, i + 1);

    /* Clipping allowed only if clipping option is specified, partial
     * image is allowed, clipping depth is meaningful and sizes are not
     * under control.
     */
    if ((PIoption != NULL)  &&
        (partialImageAllowed) &&
        (partialImageMethod == Img_PIClipping_c) &&
        (nvars >= clippingDepth)) {
      /* if conjuncts too small, dont clip. */
      if ((bdd_size(product) > partialImageThreshold) ||
          (bdd_size(relation) > partialImageThreshold)) {
        allowClipping = TRUE;
      }
    }

    if (lazySiftFlag) {
      int j;
      int nVars = bdd_num_vars(mddManager);
      var_set_t *sup_ids = bdd_get_support(relation);

      for (j = 0; j < nVars; j++) {
        if (var_set_get_elt(sup_ids, j) == 1) {
          if (bdd_is_ns_var(mddManager, j))
            bdd_reset_var_to_be_grouped(mddManager, j);
          else if (!bdd_is_ps_var(mddManager, j))
            bdd_reset_var_to_be_grouped(mddManager, j);
        }
      }
      var_set_free(sup_ids);
    }

    /* if clipping not allowed, do regular and-abstract */
    if (!allowClipping) {
      if (smoothVarCubeArray) {
        if (bdd_is_tautology(smoothVarCube, 1))
          tmpProduct = bdd_and(product, relation, 1, 1);
        else {
          tmpProduct = bdd_and_smooth_with_cube(product, relation,
                                                smoothVarCube);
        }
      } else {
        if (array_n(smoothVarBddArray) == 0)
          tmpProduct = bdd_and(product, relation, 1, 1);
        else
          tmpProduct = bdd_and_smooth(product, relation, smoothVarBddArray);
      }
    } else {
      /* the conditions to clip are that partial image is allowed,
       * partial image method is set to clipping, clipping depth
       * is meaningful and either the product or the relation is
       * larger than the partial image threshold.
       */
      if (smoothVarCubeArray)
        smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, i+1);
      tmpProduct = ComputeClippedAndAbstract(product, relation,
                                             smoothVarBddArray, nvars,
                                             clippingDepth, partial, verbosity);
    }

    if (lazySiftFlag) {
      int k, index;
      bdd_t *svar;

      if (smoothVarCubeArray)
        smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, i+1);
      for (k = 0; k < array_n(smoothVarBddArray); k++) {
        svar = array_fetch(bdd_t *, smoothVarBddArray, k);
        index = bdd_top_var_id(svar);
        bdd_set_var_to_be_grouped(mddManager, index);
      }
    }

    bdd_free(product);
    product = tmpProduct;

    /* subset if necessary */
    if ((i != array_n(relationArray)) &&
        (partialImageAllowed) &&
        (partialImageMethod == Img_PIApprox_c)) {
      intermediateSize = bdd_size(product);
      /* approximate if the intermediate product is too large. */
      if (intermediateSize > partialImageThreshold) {
        if (verbosity >= 2) {
          (void) fprintf(vis_stdout, "Intermediate Block %d Size = %10d\n", i,
                         intermediateSize);
        }
        newProduct = ComputeSubsetOfIntermediateProduct(product, PIoption);
        if (verbosity >= 2) {
          if (!bdd_equal(newProduct, product)) {
            (void)fprintf(vis_stdout,
                          "Intermediate Block Subset Size = %10d\n",
                          bdd_size(newProduct));
          } else {
            (void)fprintf(vis_stdout, "Intermediate Block unchanged\n");
          }
        }
        /* record if partial image occurred */
        if (!bdd_equal(newProduct, product)) {
          *partial = TRUE;
          intermediateSize = bdd_size(product);
        }
        bdd_free(product);
        product = newProduct;
      }
    }

    /* keep track of intermediate product size */
    if (verbosity >= 2) {
      if ((partialImageMethod  != Img_PIApprox_c) ||
          (!partialImageAllowed) ||
          ( i == array_n(relationArray))) {
        intermediateSize = bdd_size(product);
      }
      if (maxSize < intermediateSize) {
        maxSize = intermediateSize;
      }
    }
  }
  if (lazySiftFlag) {
    int nVars = bdd_num_vars(mddManager);
    for (i = 0; i < nVars; i++) {
      if (bdd_is_ps_var(mddManager, i))
        bdd_reset_var_to_be_grouped(mddManager, i);
      else
        bdd_set_var_to_be_grouped(mddManager, i);
    }
  }
  if (verbosity >= 2) {
    (void)fprintf(vis_stdout,
                  "Max. BDD size for intermediate product = %10d\n", maxSize);
  }/* for all clusters */
  
  return product;
} /* end of BddLinearAndSmooth */

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

  Synopsis    [Hashes bdd id to bdd in the table.]

  Description [The elements of the input array should be in one to one
  correpondence, i.e., idArray[i] should be the id of the variable
  bddArray[i]. Each element of idArray is hashed in the table with the value
  taken from the corresponding element from bddArray.]

  SideEffects []

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

static void
HashIdToBddTable(st_table *table, array_t *idArray,
                 array_t *bddArray)
{
  int i;
  for (i=0; i<array_n(idArray); i++) {
    int id;
    bdd_t *bdd;
    id = array_fetch(int, idArray, i);
    bdd = array_fetch(bdd_t*, bddArray, i);
    st_insert(table, (char*)(long)id, (char *)bdd);
  }
}


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

  Synopsis    [Prints information about the schedule of clusters and
  the the corresponding smooth cubes.]

  Description [This function is used to print information about the
  cluster sequence and the sequence of smooth cubes. For each cluster
  in sequence, it prints the number of variables quantified and the
  number of variables introduced.]

  SideEffects []

******************************************************************************/
static void
PrintSmoothIntroducedCount(array_t *clusterArray, array_t
                           **arraySmoothVarBddArrayPtr,
                           array_t *psBddIdArray, array_t
                           *nsBddIdArray)
{
  int i,j;
  st_table *nsBddIdTable, *supportTable, *psBddIdTable;
  bdd_t *trans;
  int introducedCount;
  bdd_t *tmp;
  st_generator *stgen;
  long varId;
  array_t *smoothVarBddArray, *arraySmoothVarBddArray;
  array_t *smoothVarBddIdArray = NIL(array_t);
  int psCount, piNdCount;

  if (arraySmoothVarBddArrayPtr == NIL(array_t*))
    arraySmoothVarBddArray = NIL(array_t);
  else
    arraySmoothVarBddArray = *arraySmoothVarBddArrayPtr;

  (void) fprintf(vis_stdout, "**********************************************\n");
  nsBddIdTable = st_init_table(st_numcmp, st_numhash);
  for (i=0; i<array_n(nsBddIdArray); i++)
    st_insert(nsBddIdTable, (char *)(long) array_fetch(int, nsBddIdArray, i),
              (char *) NULL);
  psBddIdTable = st_init_table(st_numcmp, st_numhash);
  for (i=0; i<array_n(psBddIdArray); i++)
    st_insert(psBddIdTable, (char *)(long) array_fetch(int, psBddIdArray, i),
              (char *) NULL);

  for (i=0; i<=array_n(clusterArray)-1; i++) {
    trans = array_fetch(bdd_t*, clusterArray, i);
    supportTable = ImgBddGetSupportIdTable(trans);
    psCount = 0;
    piNdCount = 0;
    if (arraySmoothVarBddArray != NIL(array_t)) {
      smoothVarBddArray = array_fetch(array_t*, arraySmoothVarBddArray, i);
      smoothVarBddIdArray = bdd_get_varids(smoothVarBddArray);
      for (j=0; j<array_n(smoothVarBddIdArray);j++) {
        if (st_is_member(psBddIdTable, (char *)(long)
                         array_fetch(int, smoothVarBddIdArray, j)))
          psCount++;
        else
          piNdCount++;
      }
    }
    introducedCount = 0;
    st_foreach_item(nsBddIdTable, stgen, &varId, &tmp) {
      if (st_is_member(supportTable, (char *)varId)) {
        introducedCount++;
        st_delete(nsBddIdTable,&varId, NULL);
      }
    }
    (void)fprintf(vis_stdout,
    "Tr no = %3d\tSmooth PS # = %3d\tSmooth PI # = %3d\tIntroduced # = %d\n", i,
                  psCount, piNdCount, introducedCount);
    st_free_table(supportTable);
    array_free(smoothVarBddIdArray);
  }
  st_free_table(nsBddIdTable);
  st_free_table(psBddIdTable);
}

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

  Synopsis    [Prints Ids of an array of BDDs.]

  Description [Prints Ids of an array of BDDs.]

  SideEffects []

******************************************************************************/
static void
PrintBddIdFromBddArray(array_t *bddArray)
{
  int i;
  array_t *idArray = bdd_get_varids(bddArray);
  fprintf(vis_stdout,
          "**************** printing bdd ids from the bdd array ***********\n");
  fprintf(vis_stdout,"%d::\t", array_n(idArray));
  for (i=0;i<array_n(idArray); i++) {
    fprintf(vis_stdout," %d ", (int)array_fetch(int, idArray, i));
  }
  fprintf(vis_stdout,"\n******************\n");
  array_free(idArray);
  return;
}

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

  Synopsis    [Prints the integers from a symbol table.]

  Description [Prints the integers from a symbol table.]

  SideEffects []

******************************************************************************/
static void
PrintBddIdTable(st_table *idTable)
{
  st_generator *stgen;
  long varId;
  fprintf(vis_stdout,
          "**************** printing bdd ids from the id table  ***********\n");
  fprintf(vis_stdout,"%d::\t", st_count(idTable));
  st_foreach_item(idTable, stgen, &varId, NIL(char *)) {
    fprintf(vis_stdout," %d ", (int) varId);
  }
  fprintf(vis_stdout,"\n******************\n");
  return;
}

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

  Synopsis    [Traverses the partition recursively, creating the
  relations for the intermediate vertices.]

  Description [Traverses the partition recursively, creating the
  relations for the intermediate vertices.]

  SideEffects []

******************************************************************************/
static void
PartitionTraverseRecursively(mdd_manager *mddManager, vertex_t *vertex,
                             int mddId, st_table *vertexTable,
                             array_t *relationArray,
                             st_table *domainQuantifyVarMddIdTable,
                             array_t *quantifyVarMddIdArray)
{
  Mvf_Function_t *mvf;
  lsList faninEdges;
  lsGen gen;
  vertex_t *faninVertex;
  edge_t *faninEdge;
  array_t *varBddRelationArray;

  if (st_is_member(vertexTable, (char *)vertex)) return;
  st_insert(vertexTable, (char *)vertex, NIL(char));
  if (mddId != -1) { /* This is not the next state function node */
    /* There is an mdd variable associated with this vertex */
    array_insert_last(int, quantifyVarMddIdArray, mddId);
    mvf = Part_VertexReadFunction(vertex);
    /*relation = Mvf_FunctionBuildRelationWithVariable(mvf, mddId);*/
    varBddRelationArray = mdd_fn_array_to_bdd_rel_array(mddManager, mddId, mvf);
    array_append(relationArray, varBddRelationArray);
    array_free(varBddRelationArray);
    /*array_insert_last(mdd_t *, relationArray, relation);*/
  }
  faninEdges = g_get_in_edges(vertex);
  if (lsLength(faninEdges) == 0) return;
  lsForEachItem(faninEdges, gen, faninEdge) {
    faninVertex = g_e_source(faninEdge);
    mddId = Part_VertexReadMddId(faninVertex);
    if (st_is_member(domainQuantifyVarMddIdTable, (char *)(long)mddId)) {
      /* This is either domain var or quantify var */
      /* continue */
      continue;
    }
    PartitionTraverseRecursively(mddManager, faninVertex, mddId, vertexTable,
                                 relationArray,
                                 domainQuantifyVarMddIdTable,
                                 quantifyVarMddIdArray);
  }
}


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

  Synopsis    [Prints the partition recursively.]

  Description [Prints the partition recursively.]

  SideEffects []

******************************************************************************/
static void
PrintPartitionRecursively(vertex_t *vertex, st_table *vertexTable, int indent)
{
  int i;
  lsList faninEdges;
  lsGen gen;
  vertex_t *faninVertex;
  edge_t *faninEdge;

  if (st_is_member(vertexTable, (char *)vertex)) return;
  st_insert(vertexTable, (char *)vertex, NIL(char));
  for(i=0; i<= indent; i++) fprintf(vis_stdout," ");
  fprintf(vis_stdout,"%s %d\n", Part_VertexReadName(vertex),
          Part_VertexReadMddId(vertex));
  faninEdges = g_get_in_edges(vertex);
  if (lsLength(faninEdges) == 0) return;
  lsForEachItem(faninEdges, gen, faninEdge) {
    faninVertex = g_e_source(faninEdge);
    PrintPartitionRecursively(faninVertex, vertexTable,indent+2);
  }
}

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

  Synopsis    [Recompute partial image by relaxing the parameters.]

  Description [Recompute partial image by relaxing the
  parameters. First the threshold to trigger approximation is
  increased and then the threshold to approximate is increased. If
  both fail to produce new states, the exact image is computed (no
  approximations).Partial flag comes in indicating whether partial
  images are allowed, and leave indicating whether the image was
  partial or not. ASSUME that toCareSet is in terms of
  function->domainVars.]

  SideEffects [partial flag is modified with each image
  computation. IMPORTANT NOTE: image may be in terms of domain
  variables even when the reverse is expected.]

******************************************************************************/
static bdd_t *
RecomputeImageIfNecessary(ImgFunctionData_t *functionData,
                          mdd_manager *mddManager,
                          bdd_t *domainSubset,
                          array_t *relationArray,
                          array_t *arraySmoothVarBddArray,
                          array_t *smoothVarCubeArray,
                          int verbosity,
                          ImgPartialImageOption_t *PIoption,
                          array_t *toCareSetArray,
                          boolean *partial, boolean lazySiftFlag)
{
  int oldPIThreshold = 0, oldPISize;
  double oldClippingFrac = IMG_PI_CLIP_DEPTH;
  boolean realRequired = FALSE;
  mdd_t *image = NIL(mdd_t);

  /* relax values to compute a non-trivial partial image */
  if (PIoption->partialImageMethod == Img_PIApprox_c) {
    oldPIThreshold = PIoption->partialImageThreshold;
    PIoption->partialImageThreshold *= 2;
    if (!PIoption->partialImageThreshold)
      PIoption->partialImageThreshold = IMG_PI_SP_THRESHOLD;
  } else if (PIoption->partialImageMethod == Img_PIClipping_c) {
    oldClippingFrac = PIoption->clippingDepthFrac;
    PIoption->clippingDepthFrac =
        (PIoption->clippingDepthFrac < IMG_PI_CLIP_DEPTH_FINAL) ?
        IMG_PI_CLIP_DEPTH_FINAL : PIoption->clippingDepthFrac;
  }
  if (((PIoption->partialImageMethod == Img_PIClipping_c) &&
       (oldClippingFrac != PIoption->clippingDepthFrac)) ||
      (PIoption->partialImageMethod == Img_PIApprox_c)) {
    if (verbosity >= 2) {
      fprintf(vis_stdout,
        "IMG: No new states, computing image again with relaxed thresholds.\n");
      if (PIoption->partialImageMethod == Img_PIApprox_c) {
        fprintf(vis_stdout, "IMG: Relaxed hd_partial_image_threshold to %d\n",
                PIoption->partialImageThreshold);
      } else if (PIoption->partialImageMethod == Img_PIClipping_c) {
        fprintf(vis_stdout,
                "IMG: Relaxed hd_partial_image_clipping_depth to %g\n",
                PIoption->clippingDepthFrac);
      }
    }

    *partial = TRUE;
    /* if values were relaxed, compute image again. */
    image = BddLinearAndSmooth(mddManager, domainSubset, relationArray,
                                arraySmoothVarBddArray, smoothVarCubeArray,
                                verbosity, PIoption, partial, lazySiftFlag);

    /* check if new states were produced. */
    if (*partial) {
      if (bdd_leq_array(image, toCareSetArray, 1, 0)) {
        bdd_free(image);
        if (PIoption->partialImageMethod == Img_PIApprox_c) {
          /* relax approximation values some more. */
          oldPISize = PIoption->partialImageSize;
          PIoption->partialImageSize *= 2;
          if (!PIoption->partialImageSize)
            PIoption->partialImageSize = IMG_PI_SP_THRESHOLD;

          if (verbosity >= 2) {
            fprintf(vis_stdout,
        "IMG: No new states, computing image again with relaxed thresholds.\n");
            fprintf(vis_stdout, "IMG: Relaxed hd_partial_image_size to %d\n",
                    PIoption->partialImageSize);
          }

          *partial = TRUE;
          image = BddLinearAndSmooth(mddManager, domainSubset, relationArray,
                                     arraySmoothVarBddArray, smoothVarCubeArray,
                                     verbosity, PIoption, partial,
                                     lazySiftFlag);
          PIoption->partialImageSize = oldPISize;

          /* check if new states were produced. */
          if (*partial) {
            if (bdd_leq_array(image, toCareSetArray, 1, 0)) {
              bdd_free(image);
              realRequired = TRUE;
            }
          }
        } else if (PIoption->partialImageMethod == Img_PIClipping_c) {
          /* compute real image since clipping depth cannot be increased. */
          realRequired = TRUE;
        }
      }
    } /* require recomputation */
  } else { /* clipping depth could not be increased. */
    /* compute real image since clipping depth cannot be increased. */
    realRequired = TRUE;
  }

  /* reset the values to the orginal values. */
  if (PIoption->partialImageMethod == Img_PIClipping_c) {
    PIoption->clippingDepthFrac = oldClippingFrac;
  } else if (PIoption->partialImageMethod == Img_PIApprox_c) {
    PIoption->partialImageThreshold = oldPIThreshold;
  }
  /* if the actual image is required compute it. */
  if (realRequired) {
    if (verbosity >= 2) {
      fprintf(vis_stdout, "IMG: No new states, computing exact image.\n");
    }
    *partial = FALSE;
    image = BddLinearAndSmooth(mddManager, domainSubset, relationArray,
                                arraySmoothVarBddArray, smoothVarCubeArray,
                                verbosity, PIoption, partial, lazySiftFlag);
    assert(!(*partial));
  }

  return (image);
} /* end of RecomputeImageIfNecessary */


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

  Synopsis    [Approximate intermediate product according to method.]

  Description [Approximate intermediate product according to method.]

  SideEffects []

******************************************************************************/
static bdd_t *
ComputeSubsetOfIntermediateProduct(bdd_t *product,
                                   ImgPartialImageOption_t *PIoption)
{
  int partialImageSize = PIoption->partialImageSize;
  int nvars = PIoption->nvars;
  bdd_t *newProduct;

  switch (PIoption->partialImageApproxMethod) {
  case BDD_APPROX_HB:
    {
      int tempSize;
      tempSize = (partialImageSize < IMG_PI_SP_THRESHOLD) ?
                IMG_PI_SP_THRESHOLD : partialImageSize;
      newProduct = bdd_approx_hb(product,
                                 (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                 nvars, tempSize);
      break;
    }
  case BDD_APPROX_SP:
    {
      int tempSize;
      tempSize = (partialImageSize < IMG_PI_SP_THRESHOLD) ?
                IMG_PI_SP_THRESHOLD : partialImageSize;
      newProduct = bdd_approx_sp(product,
                                 (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                 nvars, tempSize, 0);
      break;
    }
  case BDD_APPROX_UA:
    newProduct = bdd_approx_ua(product,
                               (bdd_approx_dir_t)BDD_UNDER_APPROX,
                               nvars, partialImageSize,
                               0, PIoption->quality);
    break;
  case BDD_APPROX_RUA:
    newProduct = bdd_approx_remap_ua(product,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     0, partialImageSize,
                                     PIoption->quality);
    break;
  case BDD_APPROX_COMP:
    {
      newProduct = bdd_approx_compress(product,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     nvars, partialImageSize);
    }
    break;
  case BDD_APPROX_BIASED_RUA:
    newProduct = bdd_approx_biased_rua(product,
                                       (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                       PIoption->bias,
                                       0, partialImageSize,
                                       PIoption->quality,
                                       PIoption->qualityBias);
    break;
  default :
    newProduct = bdd_approx_remap_ua(product,
                                     (bdd_approx_dir_t)BDD_UNDER_APPROX,
                                     0, partialImageSize,
                                     PIoption->quality);
    break;
  }
  return newProduct;
} /* end of ComputeSubsetOfIntermediateProduct */

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

  Synopsis    [Clipping and-abstract]

  Description [Compute clipping and-abstract. Try twice. If fail,
  compute exact and-abstract]

  SideEffects [partial flag is modified.]

******************************************************************************/
static bdd_t *
ComputeClippedAndAbstract(bdd_t *product,
                          bdd_t *relation,
                          array_t *smoothVarBddArray,
                          int nvars,
                          int clippingDepth,
                          boolean *partial,
                          int verbosity)
{
  bdd_t *tmpProduct;

  /* compute clipped image */
  tmpProduct = bdd_clipping_and_smooth(product, relation,
                                       smoothVarBddArray, clippingDepth, 0);
  /* if product is zero, compute clipped image with increased
   * clipping depth
   */
  if (bdd_is_tautology(tmpProduct, 0)) {
    bdd_free(tmpProduct);
    tmpProduct = bdd_clipping_and_smooth(product, relation,
                                         smoothVarBddArray,
                                         (int) (nvars*IMG_PI_CLIP_DEPTH_FINAL),
                                         0);

    /* if clipped image is zero, compute unclipped image */
    if (bdd_is_tautology(tmpProduct, 0)) {
      if (verbosity >= 2) {
        fprintf(vis_stdout,
                "IMG: Clipping failed, computing exact and-abstract\n");
      }
      bdd_free(tmpProduct);
      tmpProduct = bdd_and_smooth(product, relation,
                                  smoothVarBddArray);
    } else *partial = TRUE;

  } else *partial = TRUE;

  return tmpProduct;
} /* end of ComputeClippedAndAbstract */


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

  Synopsis    [Copy an array of bdd array.]

  Description [Copy an array of bdd array.]

  SideEffects []

******************************************************************************/
static array_t *
CopyArrayBddArray(array_t *arrayBddArray)
{
  int           i, j;
  array_t       *newArrayBddArray = array_alloc(array_t *, 0);
  array_t       *bddArray, *newBddArray;
  mdd_t         *bdd;

  for (i = 0; i < array_n(arrayBddArray); i++) {
    bddArray = array_fetch(array_t *, arrayBddArray, i);
    newBddArray = array_alloc(mdd_t *, 0);
    for (j = 0; j < array_n(bddArray); j++) {
      bdd = array_fetch(mdd_t *, bddArray, j);
      array_insert(mdd_t *, newBddArray, j, bdd_dup(bdd));
    }
    array_insert(array_t *, newArrayBddArray, i, newBddArray);
  }

  return(newArrayBddArray);
}


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

  Synopsis    [Prints info of the partitioned transition relation.]

  Description [Prints info of the partitioned transition relation.]

  SideEffects []

******************************************************************************/
static void
PrintPartitionedTransitionRelation(mdd_manager *mddManager,
                                   Iwls95Info_t *info,
                                   Img_DirectionType directionType)
{
  array_t *relationArray;
  array_t *arraySmoothVarBddArray;

  if (directionType == Img_Forward_c) {
    relationArray = info->fwdClusteredRelationArray;
    arraySmoothVarBddArray = info->fwdArraySmoothVarBddArray;
  } else if (directionType == Img_Backward_c) {
    relationArray = info->bwdClusteredRelationArray;
    arraySmoothVarBddArray = info->bwdArraySmoothVarBddArray;
  } else {
   return;
  }

  ImgPrintPartitionedTransitionRelation(mddManager, relationArray,
                                        arraySmoothVarBddArray);
}


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

  Synopsis    [Reorder partitioned transition relation.]

  Description [Reorder partitioned transition relation.]

  SideEffects []

******************************************************************************/
static void
ReorderPartitionedTransitionRelation(Iwls95Info_t *info,
                                     ImgFunctionData_t *functionData,
                                     Img_DirectionType directionType)
{
  array_t *domainVarMddIdArray = functionData->domainVars;
  array_t *rangeVarMddIdArray = functionData->rangeVars;
  array_t *quantifyVarMddIdArray = array_dup(functionData->quantifyVars);
  graph_t *mddNetwork = functionData->mddNetwork;
  mdd_manager *mddManager = Part_PartitionReadMddManager(mddNetwork);
  array_t *rangeVarBddArray, *domainVarBddArray, *quantifyVarBddArray;
  array_t *clusteredRelationArray;
  st_table *varsTable = st_init_table(st_numcmp, st_numhash);
  int i, j, k, mddId;
  array_t *smoothVarBddArray, *arraySmoothVarBddArray;
  mdd_t *mdd;
  array_t *supportIdArray;

  /* Finds the variables which are neither state variable nor primary input
   * from previous smooth variable array, and put those variable into
   * quantify variable array.
   */
  for (i =0 ; i<  array_n(domainVarMddIdArray); i++) {
    mddId = array_fetch(int, domainVarMddIdArray, i);
    st_insert(varsTable, (char *)(long)mddId, NIL(char));
  }
  for (i = 0; i < array_n(rangeVarMddIdArray); i++) {
    mddId = array_fetch(int, rangeVarMddIdArray, i);
    st_insert(varsTable, (char *)(long)mddId, NIL(char));
  }
    for (i = 0; i < array_n(quantifyVarMddIdArray); i++) {
    mddId = array_fetch(int, quantifyVarMddIdArray, i);
    st_insert(varsTable, (char *)(long)mddId, NIL(char));
  }
    if (directionType == Img_Forward_c || directionType == Img_Both_c)
    arraySmoothVarBddArray = info->fwdArraySmoothVarBddArray;
  else
    arraySmoothVarBddArray = info->bwdArraySmoothVarBddArray;
  for (i = 0; i < array_n(arraySmoothVarBddArray); i++) {
    smoothVarBddArray = array_fetch(array_t *, arraySmoothVarBddArray, i);
    for (j = 0; j < array_n(smoothVarBddArray); j++) {
      mdd = array_fetch(mdd_t *, smoothVarBddArray, j);
      supportIdArray = mdd_get_support(mddManager, mdd);
      for (k = 0; k < array_n(supportIdArray); k++) {
        mddId = array_fetch(int, supportIdArray, k);
        if (!st_lookup(varsTable, (char *)(long)mddId, NIL(char *))) {
          array_insert_last(int, quantifyVarMddIdArray, mddId);
          st_insert(varsTable, (char *)(long)mddId, NIL(char));
        }
      }
      array_free(supportIdArray);
    }
  }
  st_free_table(varsTable);

  /* Get the Bdds from Mdd Ids */
  rangeVarBddArray = functionData->rangeBddVars;
  domainVarBddArray = functionData->domainBddVars;
  quantifyVarBddArray = mdd_id_array_to_bdd_array(mddManager,
                                                  quantifyVarMddIdArray);
  array_free(quantifyVarMddIdArray);

  if (directionType == Img_Forward_c || directionType == Img_Both_c) {
    clusteredRelationArray = info->fwdClusteredRelationArray;
    mdd_array_array_free(info->fwdArraySmoothVarBddArray);
    info->fwdClusteredRelationArray = NIL(array_t);
    info->fwdArraySmoothVarBddArray = NIL(array_t);
    if (info->method == Img_Iwls95_c) {
      OrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, info->option,
                        &info->fwdClusteredRelationArray,
                        &info->fwdArraySmoothVarBddArray);
    } else if (info->method == Img_Mlp_c) {
      ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, Img_Forward_c,
                        &info->fwdClusteredRelationArray,
                        &info->fwdArraySmoothVarBddArray,
                        info->option);
    } else
      assert(0);
    mdd_array_free(clusteredRelationArray);

    if (functionData->createVarCubesFlag) {
      if (info->fwdSmoothVarCubeArray)
        mdd_array_free(info->fwdSmoothVarCubeArray);
      info->fwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->fwdArraySmoothVarBddArray);
    }

    if (info->option->verbosity > 2) {
      PrintPartitionedTransitionRelation(mddManager, info, Img_Forward_c);
    }
  }

  if (directionType == Img_Backward_c || directionType == Img_Both_c) {
    clusteredRelationArray = info->bwdClusteredRelationArray;
    mdd_array_array_free(info->bwdArraySmoothVarBddArray);
    info->bwdClusteredRelationArray = NIL(array_t);
    info->bwdArraySmoothVarBddArray = NIL(array_t);
    if (info->method == Img_Iwls95_c) {
      OrderRelationArray(mddManager, clusteredRelationArray,
                        rangeVarBddArray, quantifyVarBddArray,
                        domainVarBddArray, info->option,
                        &info->bwdClusteredRelationArray,
                        &info->bwdArraySmoothVarBddArray);
    } else if (info->method == Img_Mlp_c) {
      ImgMlpOrderRelationArray(mddManager, clusteredRelationArray,
                        domainVarBddArray, quantifyVarBddArray,
                        rangeVarBddArray, Img_Backward_c,
                        &info->bwdClusteredRelationArray,
                        &info->bwdArraySmoothVarBddArray, info->option);
    } else
      assert(0);
    mdd_array_free(clusteredRelationArray);

    if (functionData->createVarCubesFlag) {
      if (info->bwdSmoothVarCubeArray)
        mdd_array_free(info->bwdSmoothVarCubeArray);
      info->bwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->bwdArraySmoothVarBddArray);
    }

    if (info->option->verbosity > 2) {
      PrintPartitionedTransitionRelation(mddManager, info, Img_Backward_c);
    }
  }

  /* Free the Bdd arrays */
  mdd_array_free(quantifyVarBddArray);
}


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

  Synopsis    [Updates quantification schedule.]

  Description [Updates quantification schedule.]

  SideEffects []

******************************************************************************/
static void
UpdateQuantificationSchedule(Iwls95Info_t *info,
                             ImgFunctionData_t *functionData,
                             Img_DirectionType directionType)
{
  graph_t *mddNetwork = functionData->mddNetwork;
  mdd_manager *mddManager = Part_PartitionReadMddManager(mddNetwork);
  array_t *clusteredRelationArray;
  st_table *quantifyVarsTable;
  var_set_t *supportVarSet;
  int i, j, id, nRelations;
  int pos;
  long longid;
  array_t *smoothVarBddArray, **smoothVarBddArrayPtr;
  bdd_t *bdd, *relation;
  st_generator *gen;

  if (directionType == Img_Forward_c || directionType == Img_Both_c) {
    nRelations = array_n(info->fwdClusteredRelationArray);
    quantifyVarsTable = st_init_table(st_numcmp, st_numhash);
    for (i = 0; i <= nRelations; i++) {
      smoothVarBddArray = array_fetch(array_t *,
                      info->fwdArraySmoothVarBddArray, i);
      for (j = 0; j < array_n(smoothVarBddArray); j++) {
        bdd = array_fetch(mdd_t *, smoothVarBddArray, j);
        id = (int)bdd_top_var_id(bdd);
        st_insert(quantifyVarsTable, (char *)(long)id, NIL(char));
      }
    }

    clusteredRelationArray = info->fwdClusteredRelationArray;
    for (i = 0; i < nRelations; i++) {
      relation = array_fetch(bdd_t *, clusteredRelationArray, i);
      supportVarSet = bdd_get_support(relation);
      for (j = 0; j < supportVarSet->n_elts; j++) {
        if (var_set_get_elt(supportVarSet, j) == 1) {
          if (st_is_member(quantifyVarsTable, (char *)(long)j))
            st_insert(quantifyVarsTable, (char *)(long)j, (char *)(long)(i+1));
        }
      }
      var_set_free(supportVarSet);
    }

    mdd_array_array_free(info->fwdArraySmoothVarBddArray);
    info->fwdArraySmoothVarBddArray = array_alloc(array_t *, nRelations + 1);
    smoothVarBddArrayPtr = ALLOC(array_t *, sizeof(array_t *) * (nRelations+1));
    for (i = 0; i <= nRelations; i++) {
      smoothVarBddArray = array_alloc(bdd_t *, 0);
      smoothVarBddArrayPtr[i] = smoothVarBddArray;
      array_insert(array_t *, info->fwdArraySmoothVarBddArray, i,
                smoothVarBddArray);
    }

    st_foreach_item_int(quantifyVarsTable, gen, &longid, &pos) {
      id = (int) longid;
      smoothVarBddArray = array_fetch(array_t *,
                      info->fwdArraySmoothVarBddArray, pos);
      bdd = bdd_var_with_index(mddManager, id);
      array_insert_last(bdd_t *, smoothVarBddArray, bdd);
    }

    FREE(smoothVarBddArrayPtr);
    st_free_table(quantifyVarsTable);

    if (functionData->createVarCubesFlag) {
      if (info->fwdSmoothVarCubeArray)
        mdd_array_free(info->fwdSmoothVarCubeArray);
      info->fwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->fwdArraySmoothVarBddArray);
    }

    if (info->option->verbosity > 2) {
      PrintPartitionedTransitionRelation(mddManager, info, Img_Forward_c);
    }
  }

  if (directionType == Img_Backward_c || directionType == Img_Both_c) {
    nRelations = array_n(info->bwdClusteredRelationArray);
    quantifyVarsTable = st_init_table(st_numcmp, st_numhash);
    for (i = 0; i <= nRelations; i++) {
      smoothVarBddArray = array_fetch(array_t *,
                      info->bwdArraySmoothVarBddArray, i);
      for (j = 0; j < array_n(smoothVarBddArray); j++) {
        bdd = array_fetch(mdd_t *, smoothVarBddArray, j);
        id = (int)bdd_top_var_id(bdd);
        st_insert(quantifyVarsTable, (char *)(long)id, (char *)(long)0);
      }
    }

    clusteredRelationArray = info->bwdClusteredRelationArray;
    for (i = 0; i < nRelations; i++) {
      relation = array_fetch(bdd_t *, clusteredRelationArray, i);
      supportVarSet = bdd_get_support(relation);
      for (j = 0; j < supportVarSet->n_elts; j++) {
        if (var_set_get_elt(supportVarSet, j) == 1) {
          if (st_is_member(quantifyVarsTable, (char *)(long)j))
            st_insert(quantifyVarsTable, (char *)(long)j, (char *)(long)(i+1));
        }
      }
      var_set_free(supportVarSet);
    }

    mdd_array_array_free(info->bwdArraySmoothVarBddArray);
    info->bwdArraySmoothVarBddArray = array_alloc(array_t *, nRelations + 1);
    smoothVarBddArrayPtr = ALLOC(array_t *, sizeof(array_t *) * (nRelations+1));
    for (i = 0; i <= nRelations; i++) {
      smoothVarBddArray = array_alloc(bdd_t *, 0);
      smoothVarBddArrayPtr[i] = smoothVarBddArray;
      array_insert(array_t *, info->bwdArraySmoothVarBddArray, i,
                smoothVarBddArray);
    }

    st_foreach_item_int(quantifyVarsTable, gen, &longid, &pos) {
      id = (int) longid;
      smoothVarBddArray = array_fetch(array_t *,
                      info->bwdArraySmoothVarBddArray, pos);
      bdd = bdd_var_with_index(mddManager, id);
      array_insert_last(bdd_t *, smoothVarBddArray, bdd);
    }

    FREE(smoothVarBddArrayPtr);
    st_free_table(quantifyVarsTable);

    if (functionData->createVarCubesFlag) {
      if (info->bwdSmoothVarCubeArray)
        mdd_array_free(info->bwdSmoothVarCubeArray);
      info->bwdSmoothVarCubeArray = MakeSmoothVarCubeArray(mddManager,
                                        info->bwdArraySmoothVarBddArray);
    }

    if (info->option->verbosity > 2) {
      PrintPartitionedTransitionRelation(mddManager, info, Img_Backward_c);
    }
  }
}


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

  Synopsis    [Returns array of cubes of smoothing variables.]

  Description [Returns array of cubes of smoothing variables.]

  SideEffects []

******************************************************************************/
static array_t *
MakeSmoothVarCubeArray(mdd_manager *mddManager, array_t *arraySmoothVarBddArray)
{
  int           i, j;
  array_t       *smoothVarBddArray;
  array_t       *cubeArray = array_alloc(mdd_t *, 0);
  mdd_t         *cube, *var, *tmp;

  arrayForEachItem(array_t *, arraySmoothVarBddArray, i, smoothVarBddArray) {
    cube = mdd_one(mddManager);
    arrayForEachItem(mdd_t *, smoothVarBddArray, j, var) {
      tmp = mdd_and(cube, var, 1, 1);
      mdd_free(cube);
      cube = tmp;
    }
    array_insert_last(mdd_t *, cubeArray, cube);
  }
  return(cubeArray);
}


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

  Synopsis    [Eliminates dependent variables from transition relation.]

  Description [Eliminates dependent variables from a transition
  relation.  Returns a simplified copy of the given states if successful;
  NULL otherwise.]

  SideEffects [relationArray is also modified.]

  SeeAlso     []

******************************************************************************/
static mdd_t *
TrmEliminateDependVars(Img_ImageInfo_t *imageInfo, array_t *relationArray,
                        mdd_t *states, mdd_t **dependRelations)
{
  int           nDependVars;
  mdd_t         *newStates;
  Iwls95Info_t  *info = (Iwls95Info_t *) imageInfo->methodData;

  newStates = ImgTrmEliminateDependVars(&imageInfo->functionData, relationArray,
                                        states, dependRelations, &nDependVars);
  if (nDependVars) {
    if (imageInfo->verbosity > 0)
      (void)fprintf(vis_stdout, "Eliminated %d vars.\n", nDependVars);
    info->averageFoundDependVars = (info->averageFoundDependVars *
                        (float)info->nFoundDependVars + (float)nDependVars) /
                        (float)(info->nFoundDependVars + 1);
    info->nFoundDependVars++;
  }

  info->nPrevEliminatedFwd = nDependVars;
  return(newStates);
} /* end of TfmEliminateDependVars */


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

  Synopsis    [Comparison function used by qsort.]

  Description [Comparison function used by qsort to order the
  variables according to their signatures.]

  SideEffects [None]

******************************************************************************/
static int
TrmSignatureCompare(int *ptrX, int *ptrY)
{
  if (signatures[*ptrY] > signatures[*ptrX])
    return(1);
  if (signatures[*ptrY] < signatures[*ptrX])
    return(-1);
  return(0);
} /* end of TrmSignatureCompare */


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

  Synopsis    [Setups lazy sifting.]

  Description [Setups lazy sifting.]

  SideEffects [None]

******************************************************************************/
static void
SetupLazySifting(mdd_manager *mddManager, array_t *bddRelationArray,
                array_t *domainVarBddArray, array_t *quantifyVarBddArray,
                array_t *rangeVarBddArray, int verbosity)
{
  int i, id, psId, nsId, psLevel, nsLevel;
  int nVars = bdd_num_vars(mddManager);
  int nUngroup = 0;
  int nGroup = 0;
  int nConst = 0;
  int *nonLambda = ALLOC(int, nVars);
  int nLambda = 0;
  int nDepend1to1 = 0;
  int nIndepend1to1 = 0;
  int nPairPsOnce = 0;
  int pairPsId;
  int constFlag;
  int dependFlag;
  int nDependPs;
  int nDependPi;
  bdd_t* relation;
  var_set_t* supIds;
  bdd_t *bdd, *psVar, *nsVar;

  bdd_discard_all_var_groups(mddManager);
  for (i = 0; i < array_n(quantifyVarBddArray); i++) {
    bdd = array_fetch(bdd_t *, quantifyVarBddArray, i);
    id = bdd_top_var_id(bdd);
    bdd_set_pair_index(mddManager, id, id);
    bdd_set_pi_var(mddManager, id);
    bdd_reset_var_to_be_grouped(mddManager, id);
  }

  for (i = 0; i < array_n(rangeVarBddArray); i++) {
    psVar = array_fetch(bdd_t *, domainVarBddArray, i);
    nsVar = array_fetch(bdd_t *, rangeVarBddArray, i);
    psId = bdd_top_var_id(psVar);
    nsId = bdd_top_var_id(nsVar);
    bdd_set_pair_index(mddManager, psId, nsId);
    bdd_set_pair_index(mddManager, nsId, psId);
    bdd_set_ps_var(mddManager, psId);
    bdd_set_ns_var(mddManager, nsId);
    bdd_reset_var_to_be_grouped(mddManager, psId);
    bdd_set_var_to_be_grouped(mddManager, nsId);
  }

  if (array_n(domainVarBddArray) > array_n(rangeVarBddArray)) {
    for (i = array_n(rangeVarBddArray); i < array_n(domainVarBddArray); i++) {
      bdd = array_fetch(bdd_t *, domainVarBddArray, i);
      id = bdd_top_var_id(bdd);
      bdd_set_pair_index(mddManager, id, id);
      bdd_set_pi_var(mddManager, id);
      bdd_reset_var_to_be_grouped(mddManager, id);
    }
  }

  memset(nonLambda, 0, sizeof(int) * nVars);
  for (i = 0; i < array_n(bddRelationArray); i++) {
    nsId = -1;
    pairPsId = -1;
    constFlag = 1;
    dependFlag = 0;
    nDependPs = 0;
    nDependPi = 0;
    relation = array_fetch(bdd_t *, bddRelationArray, i);
    supIds = bdd_get_support(relation);

    for (id = 0; id < nVars; id++) {
      if (var_set_get_elt(supIds, id) == 1) {
        if (bdd_is_ns_var(mddManager, id)) {
          pairPsId = bdd_read_pair_index(mddManager, id);
          nsId = id;
        } else if (bdd_is_ps_var(mddManager, id)) {
          constFlag = 0;  /* nonconst next sate function */
          nonLambda[id]++;
          psId = id;
          /* At least one next state function depends on psvar id  */
          nDependPs++;
        } else {
          constFlag = 0;  /* nonconst next state function */
          nDependPi++;
        }
      }
    }
    if (nsId >= 0) {
      /* bitwise transition relation depends on some nsvar */
      for (id = 0; id < nVars; id++) {
        if (var_set_get_elt(supIds, id) == 1) {
          if (pairPsId == id) {
            /* dependendent state pair */
            nGroup++;
            dependFlag = 1;
          }
        }
      }
      if (dependFlag != 1) {
        /* independent state pair */
        nUngroup++;
        if (constFlag == 1) {
          /* next state function is constant */
          psLevel = bdd_get_level_from_id(mddManager, pairPsId);
          nsLevel = bdd_get_level_from_id(mddManager, nsId);
          if (psLevel == nsLevel - 1) {
            bdd = bdd_var_with_index(mddManager, pairPsId);
            bdd_new_var_block(bdd, 2);
            bdd_free(bdd);
            nConst++;
          } else if (psLevel == nsLevel + 1) {
            bdd = bdd_var_with_index(mddManager, nsId);
            bdd_new_var_block(bdd, 2);
            bdd_free(bdd);
            nConst++;
          }
        } else {
          bdd_set_var_to_be_ungrouped(mddManager, pairPsId);
          bdd_set_var_to_be_ungrouped(mddManager, nsId);
        }
      } else if (nDependPs == 1) {
        psLevel = bdd_get_level_from_id(mddManager, pairPsId);
        nsLevel = bdd_get_level_from_id(mddManager, nsId);
        if (psLevel == nsLevel - 1) {
          bdd = bdd_var_with_index(mddManager, pairPsId);
          bdd_new_var_block(bdd, 2);
          bdd_free(bdd);
          nDepend1to1++;
        } else if (psLevel == nsLevel + 1) {
          bdd = bdd_var_with_index(mddManager, nsId);
          bdd_new_var_block(bdd, 2);
          bdd_free(bdd);
          nDepend1to1++;
        }
      } else if (nonLambda[pairPsId] == 1) {
        psLevel = bdd_get_level_from_id(mddManager, pairPsId);
        nsLevel = bdd_get_level_from_id(mddManager, nsId);
        if (psLevel == nsLevel - 1) {
          bdd = bdd_var_with_index(mddManager, pairPsId);
          bdd_new_var_block(bdd, 2);
          bdd_free(bdd);
          nPairPsOnce++;
        } else if (psLevel == nsLevel + 1) {
          bdd = bdd_var_with_index(mddManager, nsId);
          bdd_new_var_block(bdd, 2);
          bdd_free(bdd);
          nPairPsOnce++;
        }
      }
    } else {
      /* bitwise transition relation does not depend on any nsvar */
      nUngroup++;
    }
    var_set_free(supIds);
  }
  for (i = 0; i < nVars; i++) {
    if (bdd_is_ps_var(mddManager, i)) {
      if (nonLambda[i] == 0) {
        /* no next state function depends on psvar i */
        psId = i;
        nsId = bdd_read_pair_index(mddManager, psId);
        psLevel = bdd_get_level_from_id(mddManager, psId);
        nsLevel = bdd_get_level_from_id(mddManager, nsId);
        nLambda++;
        if (psLevel == nsLevel - 1) {
          bdd = bdd_var_with_index(mddManager, psId);
          bdd_new_var_block(bdd, 2);
          bdd_free(bdd);
        } else if (psLevel == nsLevel + 1) {
          bdd = bdd_var_with_index(mddManager, nsId);
          bdd_new_var_block(bdd, 2);
          bdd_free(bdd);
        }
      }
    } else if (bdd_is_pi_var(mddManager, i)) {
      bdd_set_var_to_be_ungrouped(mddManager, i);
    }
  }
  if (verbosity) {
    fprintf(vis_stdout, "#grp=%d(#depend1to1=%d,#pairpsonce=%d)",
            nGroup, nDepend1to1, nPairPsOnce);
    fprintf(vis_stdout,
            " #ungrp=%d(#const=%d,#lambda=%d, independ1to1=%d)\n",
            nUngroup, nConst, nLambda, nIndepend1to1);
  }
  FREE(nonLambda);
}

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

  Synopsis    [Check whether the given transition relation is optimized.]

  Description [Check whether the given transition relation is optimized.]

  SideEffects [ ]

******************************************************************************/
int
Img_IsTransitionRelationOptimized(Img_ImageInfo_t * imageInfo)
{
  Iwls95Info_t *info = (Iwls95Info_t *)imageInfo->methodData;
  if(info->fwdOptClusteredRelationArray &&
     info->fwdClusteredRelationArray)   return(1);
  return(0);
}

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

  Synopsis    [Compare the size of two MDDs.]

  Description [Compare the size of two MDDs.]

  SideEffects [ ]

******************************************************************************/
static int
MddSizeCompare(const void *ptr1, const void *ptr2)
{
  mdd_t *mdd1 = *(mdd_t **)ptr1;
  mdd_t *mdd2 = *(mdd_t **)ptr2;
  int size1 = mdd_size(mdd1);
  int size2 = mdd_size(mdd2);

  if (size1 > size2)
    return 1;
  else if (size1 < size2)
    return -1;
  else
    return 0;
}

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

  Synopsis    [Restore original transition relations from restricted transition relations by conjoining  with winning strategy]

  Description [Restore original transition relations from restricted transition relations by conjoining  with winning strategy]

  SideEffects []

  SeeAlso []
******************************************************************************/
void
Img_ForwardImageInfoRecoverFromWinningStrategy(
        Img_ImageInfo_t * imageInfo
        )
{
  Iwls95Info_t *info;
  array_t *new_;

  info = (Iwls95Info_t *)imageInfo->methodData;

  if(info->fwdOriClusteredRelationArray) {
    new_ = info->fwdClusteredRelationArray;
    info->fwdClusteredRelationArray = info->fwdOriClusteredRelationArray;
    info->fwdOriClusteredRelationArray = 0;
    mdd_array_free(new_);
  }
}

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

  Synopsis    [Create restricted transition relations by conjoining  with winning strategy]

  Description [Create restricted transition relations by conjoining  with winning strategy.]

  SideEffects []

  SeeAlso []
******************************************************************************/
void
Img_ForwardImageInfoConjoinWithWinningStrategy(
        Img_ImageInfo_t * imageInfo,
        mdd_t *winningStrategy)
{
  Iwls95Info_t *info;
  array_t *old, *new_;
  int i;
  mdd_t *tr, *ntr, *winning;


  old = 0;
  winning = 0;
  new_ = array_alloc(mdd_t *, 0);
  info = (Iwls95Info_t *)imageInfo->methodData;
  old = info->fwdClusteredRelationArray;
  info->fwdClusteredRelationArray = new_;
  info->fwdOriClusteredRelationArray = old;
  winning = mdd_dup(winningStrategy);
/*    winning = ImgSubstitute(winningStrategy, &(imageInfo->functionData), Img_D2R_c);*/
  for(i=0; i<array_n(old); i++) {
    tr = array_fetch(mdd_t *, old, i);
    ntr = mdd_and(tr, winning, 1, 1);
    array_insert_last(mdd_t *, new_, ntr);
  }
  mdd_free(winning);
}