convert_to_bdd.c

#include "mdb.h"
#include "scan.h"
#include "SDDS.h"
#include "SDDSaps.h"
#include <sys/types.h>
 
#define CLO_GOOD_ELEMENTS 0
#define CLO_BAD_ELEMENTS 1
#define CLO_FAILED_SUB_TREES 2
#define CLO_PIPE 3
#define CLO_VERBOSE 4
#define N_OPTIONS 5
char *option[N_OPTIONS] = {"goodElements", "badElements", "failedSubTrees", "pipe", "verbose"};
 
static char *USAGE = "covnert_to_bdd <database file> <output file> [-pipe=[input][,output]] \n\
[-goodElements=<list of base IDs>] [-badElements=<list of base IDs] [-verbose]\n\
[-failedSubTrees=<list of sub tree IDs>\n\
<database file>     fault tree database in SDDS format.\n\
<output file>       the output file contains the fault probablity of each base elements in the fault sub tree.\n\
-pipe=in|out        pipe option.\n\
-goodElements       provide list of base elements which are known to be good, i.e., their fault probability is 0.\n\
-badElements        provide list of base elements which are known (found to be bad), i.e., their fault probability is 1.\n\
-failedSubTrees      provide list of sub trees that failed. By default, compute each sub-tree in the database.\n\n";
 
typedef struct
{
  long id;
  double probability, ps, pes, MIF, DIF;
  char *label;
  char *guidance, *description;
} BASE;
 
typedef struct
{
  BASE *base;  /*first element is base */
  void *left;  /*second element is a pointer, could be base or ITE */
  void *right; /*first element is pointer,  couble be base or ITE */
  /*left=NULL, means 1; right = NULL, means 0 */
  short is_base;
  double probability, p1, p0, ps, pes;
  double MIF, DIF; /*pe1, pe0, pes for NO probability */
  char *description, *guidance, *label;
  long ID;
} ITE;
 
/*a base x can write as ITE(x, 1, 0) both left pointer and right point are NULLs*/
 
typedef struct
{
  ITE **ite_ptr;
  long ites;
  long *tree_ID;
  int32_t ID;
  char *description, *tree_name, *typeDesc;
  short type;     /*and for 0, or for 1 */
  short all_base; /*the member of sub-tree are all base elements */
  short calculated;
  double p1, p0, ps;
  ITE *CAL_ITE;
} SUB_TREE;
 
ITE *bdd_ite_cal(ITE *ite1, ITE *ite2, short type); /*type=0 for AND, and type=1 for or */
ITE *bdd_base_ite_cal(BASE *base, ITE *ite, short type);
ITE *bdd_base_base_cal(BASE *base1, BASE *base2, short type);
void load_data_base(char *filename);
void print_sub_tree(ITE *ite);
short is_base(ITE *ite);
void free_cal_ite(ITE *ite);
void locate_tree_id();
double compute_ps(ITE *ite);
void SetupOutputFile(char *outputFile, SDDS_DATASET *outData);
 
#define STACKSIZE 2000
long treeStack1[STACKSIZE];
long treeStackptr1 = 0;
long treeStack[STACKSIZE];
long treeStackptr = 0;
 
long push_node(long ID);
long pop_node(void);
void push_sub_tree(ITE *ite);
void push_sub_tree1(ITE *ite);
void compute_sub_tree_ps(ITE *ite, SDDS_DATASET *outData, long treeIndex);
long push_node1(long ID);
long pop_node1(void);
 
static BASE **base = NULL;
static long bases = 0;
static SUB_TREE *sub_tree = NULL;
static long sub_trees = 0;
static ITE **ite = NULL;
static long ites = 0;
 
static ITE **ite_ptr = NULL;
static long ite_ptrs = 0;
static long verbose = 0;
 
int main(int argc, char **argv) {
  SCANNED_ARG *s_arg;
  SDDS_DATASET outData;
  long i, j, alldone = 0, ready, tree_ID;
  char *inputFile = NULL, *outputFile = NULL;
  ITE *ite1, *ite2;
  char **badBase = NULL, **goodBase = NULL, **failedTree = NULL;
  long badBases = 0, goodBases = 0, failedTrees = 0, i_arg, tmpfile_used = 0, compute = 0;
  unsigned long pipeFlags = 0;
 
  argc = scanargs(&s_arg, argc, argv);
  if (argc < 3)
    bomb(NULL, USAGE);
 
  for (i_arg = 1; i_arg < argc; i_arg++) {
    if (s_arg[i_arg].arg_type == OPTION) {
      delete_chars(s_arg[i_arg].list[0], "_");
      switch (match_string(s_arg[i_arg].list[0], option, N_OPTIONS, 0)) {
      case CLO_GOOD_ELEMENTS:
        goodBases = s_arg[i_arg].n_items - 1;
        goodBase = malloc(sizeof(*goodBase) * goodBases);
        for (i = 0; i < goodBases; i++)
          goodBase[i] = s_arg[i_arg].list[i + 1];
        break;
      case CLO_BAD_ELEMENTS:
        badBases = s_arg[i_arg].n_items - 1;
        badBase = malloc(sizeof(*badBase) * badBases);
        for (i = 0; i < badBases; i++)
          badBase[i] = s_arg[i_arg].list[i + 1];
        break;
      case CLO_FAILED_SUB_TREES:
        failedTrees = s_arg[i_arg].n_items - 1;
        failedTree = malloc(sizeof(*failedTree) * failedTrees);
        for (i = 0; i < failedTrees; i++)
          failedTree[i] = s_arg[i_arg].list[i + 1];
        break;
      case CLO_PIPE:
        if (!processPipeOption(s_arg[i_arg].list + 1, s_arg[i_arg].n_items - 1, &pipeFlags))
          SDDS_Bomb("invalid -pipe syntax");
        break;
      case CLO_VERBOSE:
        verbose = 1;
        break;
      default:
        fprintf(stderr, "Uknown option %s provided.\n", s_arg[i_arg].list[0]);
        exit(1);
      }
    } else {
      if (!inputFile)
        inputFile = s_arg[i_arg].list[0];
      else if (!outputFile)
        outputFile = s_arg[i_arg].list[0];
      else {
        fprintf(stderr, "argument %s is invalid: too many filenames\n", s_arg[i_arg].list[0]);
        exit(1);
      }
    }
  }
  processFilenames("convert_to_bdd", &inputFile, &outputFile, pipeFlags, 1, &tmpfile_used);
  inputFile = argv[1];
  load_data_base(inputFile);
  locate_tree_id();
  if (goodBases) {
    for (i = 0; i < goodBases; i++) {
      for (j = 0; j < bases; j++) {
        if (strcmp(goodBase[i], base[j]->label) == 0) {
          base[j]->probability = 0;
          break;
        }
      }
    }
    free(goodBase);
  }
  if (badBases) {
    for (i = 0; i < badBases; i++) {
      for (j = 0; j < bases; j++) {
        if (strcmp(badBase[i], base[j]->label) == 0) {
          base[j]->probability = 1;
          break;
        }
      }
    }
    free(badBase);
  }
  SetupOutputFile(outputFile, &outData);
  for (i = 0; i < sub_trees; i++) {
    compute = 0;
    if (sub_tree[i].all_base) {
      if (verbose) {
        fprintf(stdout, "\nsub tree name %s, ID %d, the ite is\n", sub_tree[i].tree_name, sub_tree[i].ID);
        print_sub_tree(sub_tree[i].CAL_ITE);
      }
      if (failedTrees) {
        for (j = 0; j < failedTrees; j++) {
          if (strcmp(sub_tree[i].tree_name, failedTree[j]) == 0) {
            compute = 1;
            break;
          }
        }
      } else
        compute = 1;
      if (compute)
        compute_sub_tree_ps(sub_tree[i].CAL_ITE, &outData, i);
    }
  }
 
  while (1) {
    alldone = 1;
    for (i = 0; i < sub_trees; i++) {
      if (sub_tree[i].calculated)
        continue;
      else {
        alldone = 0;
        ready = 1;
        for (j = 0; j < sub_tree[i].ites; j++) {
          tree_ID = sub_tree[i].tree_ID[j];
          if (tree_ID >= 0 && sub_tree[tree_ID].calculated == 0) {
            ready = 0;
            break;
          }
        }
        if (ready) {
          if (sub_tree[i].ites == 1)
            sub_tree[i].CAL_ITE = sub_tree[i].ite_ptr[0];
          else {
            if (sub_tree[i].ite_ptr[0]->base && is_base(sub_tree[i].ite_ptr[0]))
              ite1 = sub_tree[i].ite_ptr[0];
            else
              ite1 = sub_tree[sub_tree[i].tree_ID[0]].CAL_ITE;
            if (sub_tree[i].ite_ptr[1]->base && is_base(sub_tree[i].ite_ptr[1]))
              ite2 = sub_tree[i].ite_ptr[1];
            else
              ite2 = sub_tree[sub_tree[i].tree_ID[1]].CAL_ITE;
            sub_tree[i].CAL_ITE = bdd_ite_cal(ite1, ite2, sub_tree[i].type);
            for (j = 2; j < sub_tree[i].ites; j++) {
              if (sub_tree[i].ite_ptr[j]->base && is_base(sub_tree[i].ite_ptr[j]))
                ite1 = sub_tree[i].ite_ptr[j];
              else
                ite1 = sub_tree[sub_tree[i].tree_ID[j]].CAL_ITE;
              sub_tree[i].CAL_ITE = bdd_ite_cal(sub_tree[i].CAL_ITE, ite1, sub_tree[i].type);
            }
            if (verbose) {
              fprintf(stdout, "\nsub tree name %s, ID %d, the ite is\n", sub_tree[i].tree_name, sub_tree[i].ID);
              print_sub_tree(sub_tree[i].CAL_ITE);
            }
            compute = 0;
            if (failedTrees) {
              for (j = 0; j < failedTrees; j++) {
                if (strcmp(sub_tree[i].tree_name, failedTree[j]) == 0) {
                  compute = 1;
                  break;
                }
              }
            } else
              compute = 1;
            if (compute)
              compute_sub_tree_ps(sub_tree[i].CAL_ITE, &outData, i);
          }
          sub_tree[i].calculated = 1;
        }
      }
    }
    if (alldone)
      break;
  }
  if (failedTrees)
    free(failedTree);
  if (!SDDS_Terminate(&outData))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  for (i = 0; i < bases; i++) {
    free(base[i]->guidance);
    base[i]->guidance = NULL;
    free(base[i]->label);
    base[i]->label = NULL;
    free(base[i]->description);
    base[i]->description = NULL;
    free(base[i]);
  }
  free(base);
  for (i = 0; i < ites; i++) {
    if (ite[i]->left)
      free((ITE *)(ite[i]->left));
    ite[i]->left = NULL;
    if (ite[i]->right)
      free((ITE *)(ite[i]->right));
    ite[i]->right = NULL;
    if (ite[i]->description)
      free(ite[i]->description);
    ite[i]->description = NULL;
    if (ite[i]->guidance)
      free(ite[i]->guidance);
    ite[i]->guidance = NULL;
    if (ite[i]->label)
      free(ite[i]->label);
    ite[i]->label = NULL;
    free(ite[i]);
    ite[i] = NULL;
  }
  free(ite);
 
  for (i = 0; i < sub_trees; i++) {
    free(sub_tree[i].description);
    free(sub_tree[i].tree_name);
    free(sub_tree[i].typeDesc);
    /* for (j=0; j<sub_tree[i].ites; j++) {
         if (sub_tree[i].ite_ptr[j]) free(sub_tree[i].ite_ptr[j]); sub_tree[i].ite_ptr[j]=NULL;
         } */
    free(sub_tree[i].ite_ptr);
    free(sub_tree[i].tree_ID);
    /* if (sub_tree[i].CAL_ITE) free(sub_tree[i].CAL_ITE); sub_tree[i].CAL_ITE = NULL; */
  }
  free(sub_tree);
  for (i = 0; i < ite_ptrs; i++) {
    if (ite_ptr[i]->guidance)
      free(ite_ptr[i]->guidance);
    if (ite_ptr[i]->description)
      free(ite_ptr[i]->description);
    if (ite_ptr[i]->label)
      free(ite_ptr[i]->label);
    free(ite_ptr[i]);
  }
  free(ite_ptr);
  free_scanargs(&s_arg, argc);
 
  return 0;
}
 
short is_base(ITE *ite) {
  if ((ite->left == NULL && ite->right == NULL) || (ite->is_base))
    return 1;
  else
    return 0;
}
 
ITE *bdd_base_base_cal(BASE *base1, BASE *base2, short type) {
  ITE *tmp = NULL;
  BASE *base_s, *base_b;
 
  tmp = calloc(sizeof(*tmp), 1);
  tmp->left = tmp->right = NULL;
  /*rememeber allocated ite pointers, in order to free them in the end */
  ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
  ite_ptr[ite_ptrs] = tmp;
  ite_ptrs++;
 
  if (base1->id == base2->id) {
    /* a+a = a; a*a = a */
    tmp->base = base1;
    tmp->is_base = 1;
    tmp->left = tmp->right = NULL;
    return tmp;
  }
  if (base1->id < base2->id) {
    base_s = base1;
    base_b = base2;
  } else {
    base_s = base2;
    base_b = base1;
  }
  if (type == 0) {
    /* and logical, returns ite(first, second, 0) */
    tmp->base = base_s;
    tmp->left = (ITE *)calloc(sizeof(ITE), 1);
    ((ITE *)tmp->left)->base = base_b;
    ((ITE *)tmp->left)->left = ((ITE *)tmp->left)->right = NULL;
 
    tmp->right = NULL;
    /*rememeber allocated ite pointers, in order to free them in the end */
    ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
    ite_ptr[ite_ptrs] = (ITE *)tmp->left;
    ite_ptrs++;
  } else {
    /* or logical ite(first, 1, second) */
    tmp->base = base_s;
    tmp->left = NULL;
    tmp->right = (ITE *)calloc(sizeof(ITE), 1);
    ((ITE *)tmp->right)->base = base_b;
    ((ITE *)tmp->right)->left = ((ITE *)tmp->right)->right = NULL;
    /*rememeber allocated ite pointers, in order to free them in the end */
    ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
    ite_ptr[ite_ptrs] = (ITE *)tmp->right;
    ite_ptrs++;
  }
  return tmp;
}
 
ITE *bdd_base_ite_cal(BASE *base, ITE *ite, short type) {
  ITE *tmp = NULL, *left, *right;
  BASE *base1 = ite->base;
 
  if (is_base(ite))
    return bdd_base_base_cal(base, ite->base, type);
 
  tmp = calloc(sizeof(*ite), 1);
  /*rememeber allocated ite pointers, in order to free them in the end */
  ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
  ite_ptr[ite_ptrs] = tmp;
  ite_ptrs++;
  left = (ITE *)ite->left;
  right = (ITE *)ite->right;
 
  if (base->id < base1->id) {
    tmp->base = base;
    if (type == 1) {
      /*or, ITE(base, 1, ite) */
      tmp->left = NULL;
      tmp->right = ite;
    } else {
      /*and  ITE(base, ite, 0) */
      tmp->left = ite;
      tmp->right = NULL;
    }
  } else if (base->id > base1->id) {
    /* ITE(ite->base, ite->left type base, ite->right type base */
    tmp->base = ite->base;
    /* ite(b,c,0) + d = ite(b,c+d, d) = ite(b,c,d) */
    /*if (!right && type==1) {
        tmp->right = (ITE*)calloc(sizeof(ITE), 1);
        ((ITE*)tmp->right)->base = base;
        ((ITE*)tmp->right)->left = ((ITE*)tmp->right)->right = NULL;
        tmp->left = left;
 
        ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr)*(ite_ptrs+1));
        ite_ptr[ite_ptrs] = (ITE*)tmp->right;
        ite_ptrs ++;
        return tmp;
        } */
    if (!left) {
      if (type == 0) {
        tmp->left = (ITE *)calloc(sizeof(ITE), 1);
        ((ITE *)tmp->left)->base = base; /* left is one 1*base = base */
        ((ITE *)tmp->left)->left = ((ITE *)tmp->left)->right = NULL;
 
        /*rememeber allocated ite pointers, in order to free them in the end */
        ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
        ite_ptr[ite_ptrs] = (ITE *)tmp->left;
        ite_ptrs++;
      } else {
        tmp->left = NULL; /* 1+base = 1 , left pointer is NULL */
      }
    } else {
      tmp->left = bdd_base_ite_cal(base, (ITE *)(ite->left), type);
    }
    if (!right) {
      /*right is 0 */
      if (type == 0)
        /* 0 * base = 0 */
        tmp->right = NULL;
      else {
        /* 0 + base = base */
        tmp->right = (ITE *)calloc(sizeof(ITE), 1);
        ((ITE *)tmp->right)->base = base;
        ((ITE *)tmp->right)->left = ((ITE *)tmp->right)->right = NULL;
 
        /*rememeber allocated ite pointers, in order to free them in the end */
        ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
        ite_ptr[ite_ptrs] = (ITE *)tmp->right;
        ite_ptrs++;
      }
    } else {
      tmp->right = bdd_base_ite_cal(base, (ITE *)ite->right, type);
    }
  } else {
    /*same */
    tmp->base = base;
    /* base, is ite(base, 1, 0) */
    if (type == 0) {
      /*and logical */
      tmp->left = left;  /* 1 * left */
      tmp->right = NULL; /* 0 * right */
    } else {
      tmp->left = NULL;   /* 1 + left = 1 */
      tmp->right = right; /* 0 + right = right */
    }
  }
  return tmp;
}
 
ITE *bdd_ite_cal(ITE *ite1, ITE *ite2, short type) {
  ITE *ite = NULL, *ite_l, *ite_r;
  BASE *base1 = NULL, *base2 = NULL;
 
  if (!ite1 || !ite2)
    SDDS_Bomb("Null pointer provided to bdd_ite_cal().");
  if (is_base(ite1) && is_base(ite2))
    return bdd_base_base_cal(ite1->base, ite2->base, type);
 
  if (is_base(ite1))
    return bdd_base_ite_cal(ite1->base, ite2, type);
  if (is_base(ite2))
    return bdd_base_ite_cal(ite2->base, ite1, type);
 
  ite = calloc(sizeof(*ite), 1);
  base1 = ite1->base;
  base2 = ite2->base;
  /*rememeber allocated ite pointers, in order to free them in the end */
  ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
  ite_ptr[ite_ptrs] = ite;
  ite_ptrs++;
  if (base1->id == base2->id) {
    ite->base = ite1->base;
    if (ite1->left && ite2->left)
      ite->left = bdd_ite_cal((ITE *)ite1->left, (ITE *)ite2->left, type);
    else {
      if (type == 1)
        ite->left = NULL; /* one left is 1, 1+any = 1 */
      else {
        /* 1 * any = any */
        if (ite1->left)
          ite->left = ite1->left;
        else
          ite->left = ite2->left;
      }
    }
    if (ite1->right && ite2->right)
      ite->right = bdd_ite_cal((ITE *)ite1->right, (ITE *)ite2->right, type);
    else {
      if (type == 0)
        ite->right = NULL; /* one left is 0,  0*any = 0 */
      else {
        /* 0 + any = any */
        if (ite1->right)
          ite->right = ite1->right;
        else
          ite->right = ite2->right;
      }
    }
  } else {
    if (base1->id < base2->id) {
      ite_l = ite1;
      ite_r = ite2;
    } else {
      ite_l = ite2;
      ite_r = ite1;
    }
    ite->base = ite_l->base;
    /* ite(b,c,0) + d = ite(b,c+d, d) = ite(b,c,d) */
    /* if (!ite_l->right && type==1) {
         ite->right = ite_r;
         ite->left = ite_l->left;
         return ite;
         } */
    if (ite_l->left)
      ite->left = bdd_ite_cal(ite_l->left, ite_r, type);
    else {
      /*ite_l->left is NULL, i.e, 1 */
      if (type == 0)
        ite->left = ite_r;
      else
        ite->left = NULL;
    }
    if (ite_l->right)
      ite->right = bdd_ite_cal(ite_l->right, ite_r, type);
    else {
      /*ite_l->left is NULL, i.e, 0 */
      if (type == 0)
        ite->right = NULL;
      else
        ite->right = ite_r;
    }
  }
  return ite;
}
 
void load_data_base(char *filename) {
  int32_t i, rows, pages, j, allbase, index;
  int32_t *id = NULL;
  double *prob = NULL;
  char **guid = NULL, **desc = NULL, **label = NULL;
  SDDS_DATASET table;
 
  pages = 0;
  if (!SDDS_InitializeInput(&table, filename))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  pages = 0;
  sub_tree = NULL;
  ite = NULL;
  while (SDDS_ReadPage(&table) > 0) {
    if (!(rows = SDDS_CountRowsOfInterest(&table)))
      continue;
    sub_tree = SDDS_Realloc(sub_tree, sizeof(*sub_tree) * (pages + 1));
    sub_tree[pages].ites = rows;
    sub_tree[pages].ite_ptr = malloc(sizeof(ITE *) * rows);
    sub_tree[pages].calculated = 0;
    if (!SDDS_GetParameter(&table, "Description", &sub_tree[pages].description) ||
        !SDDS_GetParameter(&table, "ID", &sub_tree[pages].ID) ||
        !SDDS_GetParameter(&table, "LogicalType", &sub_tree[pages].type) ||
        !SDDS_GetParameter(&table, "LogicalTypeDesc", &sub_tree[pages].typeDesc) ||
        !SDDS_GetParameter(&table, "TreeName", &sub_tree[pages].tree_name))
      SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    if (!(id = SDDS_GetColumn(&table, "ID")) ||
        !(prob = SDDS_GetColumn(&table, "Probability")) ||
        !(desc = SDDS_GetColumn(&table, "Description")) ||
        !(guid = SDDS_GetColumn(&table, "Guidance")) ||
        !(label = SDDS_GetColumn(&table, "Label")))
      SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
    allbase = 1;
    for (i = 0; i < rows; i++) {
      if (id[i] > 1000) {
        /*base element */
        index = -1;
        if (bases) {
          /*search if base already exist */
          index = -1;
          for (j = 0; j < bases; j++) {
            if (id[i] == base[j]->id) {
              index = j;
              break;
            }
          }
        }
        if (index < 0) {
          base = SDDS_Realloc(base, sizeof(*base) * (bases + 1));
          base[bases] = calloc(sizeof(**base), 1);
          base[bases]->id = id[i];
          base[bases]->probability = prob[i];
          SDDS_CopyString(&base[bases]->guidance, guid[i]);
          SDDS_CopyString(&base[bases]->label, label[i]);
          SDDS_CopyString(&base[bases]->description, desc[i]);
          index = bases;
          bases++;
        }
        sub_tree[pages].ite_ptr[i] = calloc(sizeof(ITE), 1);
        sub_tree[pages].ite_ptr[i]->base = base[index];
        sub_tree[pages].ite_ptr[i]->left = sub_tree[pages].ite_ptr[i]->right = NULL;
        sub_tree[pages].ite_ptr[i]->ID = id[i];
        sub_tree[pages].ite_ptr[i]->probability = prob[i];
        sub_tree[pages].ite_ptr[i]->is_base = 1;
        SDDS_CopyString(&sub_tree[pages].ite_ptr[i]->guidance, guid[i]);
        SDDS_CopyString(&sub_tree[pages].ite_ptr[i]->label, label[i]);
        SDDS_CopyString(&sub_tree[pages].ite_ptr[i]->description, desc[i]);
        /*rememeber allocated ite pointers, in order to free them in the end */
        ite_ptr = SDDS_Realloc(ite_ptr, sizeof(*ite_ptr) * (ite_ptrs + 1));
        ite_ptr[ite_ptrs] = sub_tree[pages].ite_ptr[i];
        ite_ptrs++;
      } else {
        /* non base -- ITE */
        index = -1;
        if (ites) {
          for (j = 0; j < ites; j++) {
            if (id[i] == ite[j]->ID) {
              index = j;
              break;
            }
          }
        }
        if (index < 0) {
          ite = SDDS_Realloc(ite, sizeof(*ite) * (ites + 1));
          ite[ites] = calloc(sizeof(ITE), 1);
          ite[ites]->ID = id[i];
          ite[ites]->base = NULL;
          ite[ites]->probability = prob[i];
          ite[ites]->is_base = 0;
          SDDS_CopyString(&ite[ites]->guidance, guid[i]);
          SDDS_CopyString(&ite[ites]->description, desc[i]);
          index = ites;
          ites++;
        }
        sub_tree[pages].ite_ptr[i] = ite[index];
        allbase = 0;
      }
    }
    sub_tree[pages].ites = rows;
    sub_tree[pages].all_base = allbase;
    /*compute base sub_tree */
    if (allbase) {
      if (rows == 1) {
        /*only have one base element, no need to do computation */
        sub_tree[pages].CAL_ITE = sub_tree[pages].ite_ptr[0];
      } else {
        sub_tree[pages].CAL_ITE = bdd_ite_cal(sub_tree[pages].ite_ptr[0], sub_tree[pages].ite_ptr[1], sub_tree[pages].type);
        for (j = 2; j < rows; j++) {
          sub_tree[pages].CAL_ITE = bdd_ite_cal(sub_tree[pages].ite_ptr[j], sub_tree[pages].CAL_ITE, sub_tree[pages].type);
        }
      }
      sub_tree[pages].calculated = 1;
    }
    free(id);
    id = NULL;
    free(prob);
    prob = NULL;
    SDDS_FreeStringArray(guid, rows);
    SDDS_FreeStringArray(desc, rows);
    SDDS_FreeStringArray(label, rows);
    free(label);
    free(guid);
    free(desc);
    pages++;
  }
  if (verbose)
    fprintf(stdout, "total bases %ld, total ites %ld\n", bases, ites);
  sub_trees = pages;
  if (!SDDS_Terminate(&table))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
}
 
long push_node(long ID) {
  if (treeStackptr >= STACKSIZE) {
    fprintf(stderr, "stack overflow.\n");
    exit(1);
  }
  treeStack[treeStackptr++] = ID;
  return 1;
}
 
long pop_node(void) {
  if (treeStackptr < 1) {
    fprintf(stderr, "too few items on stack\n");
    exit(1);
  }
  return (treeStack[--treeStackptr]);
}
 
/*pos = 0, top, pos=1, left; pos=2 right */
void print_sub_tree(ITE *ite) {
  ITE *left, *right;
  if (!ite->base)
    return;
  push_node(ite->base->id);
  fprintf(stdout, "%ld\n", pop_node());
  left = (ITE *)ite->left;
  right = (ITE *)ite->right;
  if (!left)
    push_node(1);
  else
    push_node(left->base->id);
  fprintf(stdout, "%ld ", pop_node());
  if (!right)
    push_node(0);
  else
    push_node(right->base->id);
  fprintf(stdout, "%ld \n", pop_node());
  if (left)
    print_sub_tree(left);
  if (right)
    print_sub_tree(right);
}
 
void free_cal_ite(ITE *ite) {
  ITE *left, *right;
  if (!ite)
    return;
  left = (ITE *)(ite->left);
  right = (ITE *)(ite->right);
  if (!left && !right) {
    free(ite);
    ite = NULL;
    return;
  }
  if (left) {
    if (is_base(left)) {
      free(left);
      left = NULL;
    } else {
      free_cal_ite(left);
    }
  }
  if (right) {
    if (is_base(right)) {
      free(right);
      right = NULL;
    } else {
      free_cal_ite(right);
    }
  }
}
 
void locate_tree_id() {
  long i, j, k;
  for (i = 0; i < sub_trees; i++) {
    sub_tree[i].tree_ID = NULL;
    sub_tree[i].tree_ID = malloc(sizeof(long) * sub_tree[i].ites);
    for (j = 0; j < sub_tree[i].ites; j++) {
      sub_tree[i].tree_ID[j] = -1;
      if (!(sub_tree[i].ite_ptr[j]->base)) {
        for (k = 0; k < sub_trees; k++) {
          if (sub_tree[i].ite_ptr[j]->ID == sub_tree[k].ID) {
            sub_tree[i].tree_ID[j] = k;
            break;
          }
        }
        if (sub_tree[i].tree_ID[j] < 0) {
          fprintf(stderr, "Error: no tree_ID found for ite %ld of tree %ld %s\n", j, i, sub_tree[i].tree_name);
          exit(1);
        }
      }
    }
  }
}
 
void push_sub_tree(ITE *ite) {
  ITE *left, *right;
  push_node((long)ite);
  left = ite->left;
  right = ite->right;
  if (!left)
    push_node(1);
  else
    push_node((long)left);
  if (!right)
    push_node(0);
  else
    push_node((long)right);
  if (left)
    push_sub_tree(left);
  if (right)
    push_sub_tree(right);
}
 
double compute_ps(ITE *ite) {
  ITE *left, *right, *c_ite = NULL;
  long left_adr, right_adr;
  double p;
  push_sub_tree1(ite);
  while (treeStackptr1 > 0) {
    right_adr = pop_node1();
    left_adr = pop_node1();
    c_ite = (ITE *)pop_node1();
    p = c_ite->base->probability;
    if (left_adr == 1) {
      c_ite->p1 = 1;
    } else {
      left = (ITE *)left_adr;
      c_ite->p1 = left->ps;
    }
    if (right_adr == 0) {
      c_ite->p0 = 0;
    } else {
      right = (ITE *)right_adr;
      c_ite->p0 = right->ps;
    }
    c_ite->ps = p * c_ite->p1 + (1 - p) * c_ite->p0;
  }
  return c_ite->ps;
}
 
void compute_sub_tree_ps(ITE *ite, SDDS_DATASET *outData, long treeIndex) {
  ITE *t_ite;
  long this_adr;
  long t_bases = 0, i, index;
  BASE **t_base = NULL, *base1;
  double prob, PS;
  /*first compute system PS */
  PS = compute_ps(ite);
 
  push_sub_tree(ite);
  while (treeStackptr > 0) {
    this_adr = pop_node();
    if (this_adr != 1 && this_adr != 0) {
      t_ite = (ITE *)this_adr;
      base1 = t_ite->base;
      index = -1;
      if (!t_bases) {
        t_base = SDDS_Realloc(t_base, sizeof(*t_base) * (t_bases + 1));
        t_base[t_bases] = base1;
        t_bases++;
      } else {
        for (i = 0; i < t_bases; i++) {
          if (base1->id == t_base[i]->id) {
            index = i;
            break;
          }
        }
        if (index < 0) {
          t_base = SDDS_Realloc(t_base, sizeof(*t_base) * (t_bases + 1));
          t_base[t_bases] = base1;
          t_bases++;
        }
      }
      if (index >= 0)
        continue;
      prob = base1->probability; /* remember its original probability */
      /*first compute ps when base1->probability=1 */
      base1->probability = 1.0;
      base1->ps = compute_ps(ite);
      /*then compute pes when base1->probability=0 */
      base1->probability = 0.0;
      base1->pes = compute_ps(ite);
 
      base1->probability = prob; /*restore base1->probability */
      base1->MIF = base1->ps - base1->pes;
      base1->DIF = prob + prob * (1 - prob) * base1->MIF / PS;
      if (verbose)
        fprintf(stdout, "base %ld, prob%f, ps %f, pes %f, mif %f, dif %f\n",
                base1->id, base1->probability, base1->ps, base1->pes, base1->MIF, base1->DIF);
    }
  }
  if (!SDDS_StartPage(outData, t_bases))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (!SDDS_SetParameters(outData, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE,
                          "TreeName", sub_tree[treeIndex].tree_name,
                          "Description", sub_tree[treeIndex].description,
                          "LogicalType", sub_tree[treeIndex].type,
                          "LogicalTypeDesc", sub_tree[treeIndex].typeDesc,
                          "ID", sub_tree[treeIndex].ID, NULL))
    SDDS_PrintErrors(stderr, SDDS_EXIT_PrintErrors | SDDS_VERBOSE_PrintErrors);
  for (i = 0; i < t_bases; i++) {
    if (!SDDS_SetRowValues(outData, SDDS_SET_BY_NAME | SDDS_PASS_BY_VALUE, i,
                           "BaseID", t_base[i]->id,
                           "Label", t_base[i]->label,
                           "Probability", t_base[i]->probability,
                           "DIF", t_base[i]->DIF,
                           "MIF", t_base[i]->MIF,
                           "PS", t_base[i]->ps,
                           "PES", t_base[i]->pes,
                           "Description", t_base[i]->description,
                           "Guidance", t_base[i]->guidance, NULL))
      SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  }
  if (!SDDS_WritePage(outData))
    SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (t_base)
    free(t_base);
}
 
long push_node1(long ID) {
  if (treeStackptr1 >= STACKSIZE) {
    fprintf(stderr, "stack overflow.\n");
    exit(1);
  }
  treeStack1[treeStackptr1++] = ID;
  return 1;
}
 
long pop_node1(void) {
  if (treeStackptr1 < 1) {
    fprintf(stderr, "too few items on stack\n");
    exit(1);
  }
  return (treeStack1[--treeStackptr1]);
}
 
void push_sub_tree1(ITE *ite) {
  ITE *left, *right;
  push_node1((long)ite);
  left = ite->left;
  right = ite->right;
  if (!left)
    push_node1(1);
  else
    push_node1((long)left);
  if (!right)
    push_node1(0);
  else
    push_node1((long)right);
  if (left)
    push_sub_tree1(left);
  if (right)
    push_sub_tree1(right);
}
 
void SetupOutputFile(char *outputFile, SDDS_DATASET *outData) {
  if (!SDDS_InitializeOutput(outData, SDDS_ASCII, 1, NULL, NULL, outputFile))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (!SDDS_DefineSimpleParameter(outData, "TreeName", NULL, SDDS_STRING) ||
      !SDDS_DefineSimpleParameter(outData, "Description", NULL, SDDS_STRING) ||
      !SDDS_DefineSimpleParameter(outData, "LogicalType", NULL, SDDS_SHORT) ||
      !SDDS_DefineSimpleParameter(outData, "LogicalTypeDesc", NULL, SDDS_STRING) ||
      !SDDS_DefineSimpleParameter(outData, "ID", NULL, SDDS_LONG) ||
      !SDDS_DefineSimpleColumn(outData, "BaseID", NULL, SDDS_LONG) ||
      !SDDS_DefineSimpleColumn(outData, "Label", NULL, SDDS_STRING) ||
      !SDDS_DefineSimpleColumn(outData, "Probability", NULL, SDDS_DOUBLE) ||
      !SDDS_DefineSimpleColumn(outData, "DIF", NULL, SDDS_DOUBLE) ||
      !SDDS_DefineSimpleColumn(outData, "MIF", NULL, SDDS_DOUBLE) ||
      !SDDS_DefineSimpleColumn(outData, "PS", NULL, SDDS_DOUBLE) ||
      !SDDS_DefineSimpleColumn(outData, "PES", NULL, SDDS_DOUBLE) ||
      !SDDS_DefineSimpleColumn(outData, "Description", NULL, SDDS_STRING) ||
      !SDDS_DefineSimpleColumn(outData, "Guidance", NULL, SDDS_STRING))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
  if (!SDDS_WriteLayout(outData))
    SDDS_PrintErrors(stdout, SDDS_VERBOSE_PrintErrors | SDDS_EXIT_PrintErrors);
}