#define not_in_tree 0
#define temp_in_tree 1
#define perm_in_tree 2

double TheNodeList::routerOPT(int alg, Node *source, int addr,  double &d, 
			     double &maxd, double &mind, 
			     double &h, double &nodes) {

  //An optimal branch and bound algorithm for the constrained Steiner tree 
  //problem. I believe it is optimal only for undirected graphs;

  //First delete any previous routing for the this group and source set.
  removeTree(source, addr);

  //Locate the source to get the peak rate
  SourceList *ss = source->sourceList();
  int found = False;
  while ((ss != NULL) && (found == False)) {
    if ((ss->source()->type() != Background) &&
	(ss->source()->address() == addr)) found = True;
    else ss = ss->next();
  };
  double pk = ss->source()->peak();
  double avg = ss->source()->average();
  
  //Locate the MC group that contains the destination set
  MCGroup *group = groupsHd;
  while ((group !=  NULL) && (group->address() != addr))   
    group = group->next();
  
  if (group != NULL)  {
    int grCount = group->count();
    if ((grCount == 1) && 
	(group->headm()->nodePtr()->name() == source->name())) {
      source->addRoutingEntry(addr, source);
      h = d = maxd = mind = 0;
      nodes = 1;
      return(0);
    };
    
    int found = False;
    NodeListEntry *tmp = group->headm();
    while ((tmp != NULL) && (found == False)) {
      if (tmp->nodePtr()->name() == source->name()) {
	found = True;
	grCount--;
      };
      tmp = tmp->next();
    };
    int *gr_matrix = new int[grCount];
    int i = 0;
    tmp = group->headm();
    while (tmp != NULL) {
      if (tmp->nodePtr()->name() != source->name())
	*(gr_matrix + i++) = tmp->nodePtr()->name();
      tmp = tmp->next();
    };

    int links = 0;
    tmp = nodeListHd;
    while (tmp != NULL) {
      AdjacencyListEntry *adj = tmp->nodePtr()->adjacentNodes();
      while (adj != NULL) {
	if ((((fn == PEAK) && (adj->peak() <= 
			       ((adj->linkCapacity() * ADMITRATIO) - pk))) ||
	    ((fn == AVERAGE) && (adj->average() <= 
				 ((adj->linkCapacity() * ADMITRATIO) - avg))))
	    && (adj->delay() < DELAYBOUND)) { 
	  links++;
	};
	adj = adj->next();
      };
      tmp = tmp->next();
    };
    
    int *head = new int[links];
    int *tail = new int[links];
    int *deg = new int[num];
    for (i = 0; i < num; i++) *(deg + i) = 0;
    double *delay = new double[links];
    double *cost = new double[links];

    int j = 0;
    tmp = nodeListHd;
    while (tmp != NULL) {
      i = tmp->nodePtr()->name();
      AdjacencyListEntry *adj1 = tmp->nodePtr()->adjacentNodes();
      while (adj1 != NULL) {

	if ((fn == PEAK) && (adj1->peak() <= 
			     ((adj1->linkCapacity() * ADMITRATIO) - pk))) { 
	                       //in order to eliminate links saturated links
	  if (adj1->delay() < DELAYBOUND) {
	    *(tail + j) = i;
	    *(head + j) = adj1->nodePtr()->name();
	    *(delay + j) = adj1->delay();
	    *(cost + j) = adj1->peak();
	    (*(deg + *(head + j)))++;
	    j++;
	  };
	}
	else if ((fn == AVERAGE) && (adj1->average() <= 
				((adj1->linkCapacity() * ADMITRATIO) - avg))) {
	                       //in order to eliminate links saturated links
	  if (adj1->delay() < DELAYBOUND) {				     
	    *(tail + j) = i;
	    *(head + j) = adj1->nodePtr()->name();
	    *(delay + j) = adj1->delay();
	    *(cost + j) = adj1->average();
	    (*(deg + *(head + j)))++;
	    j++;
	  };
	};
	adj1 = adj1->next();
      };
      tmp = tmp->next();
    };

    double *minwt_at_node = new double[num];
    int **adj = new  (int *)[num];
    int k;
    for (i = 0; i < num; i++) {
      k = 0;
      if (*(deg + i) > 0) *(adj + i) = new int[*(deg + i)];
      else *(adj + i) = NULL;	      
      for (j = 0; j < links; j++) {
        if (*(head + j) == i) {
	  *(*(adj + i) + k) = j;
	  k++;
	};
      };
      if (k > 1) sort(i, k, cost, *(adj + i));
      *(minwt_at_node + i) = *(cost + *(*(adj + i))); 
    };

    int *child = new int[num * grCount];
    for (i = 0; i < num; i++)
      for (j = 0; j < grCount; j++)  *(child + (i * grCount) + j) = INT_MAX;;

    double *sp = new double[num*num];
    int i0;
    double  spp;
    for (i = 0; i < num; i++)
      for (j = 0; j < num; j++) *(sp + (i * num) + j) = DBL_MAX;
    for (i = 0; i < num; i++) *(sp +(i * num) + i) = 0.0;
    
    for (j = 0; j < links; j++)
      *(sp + (*(tail + j) * num) + *(head + j)) = *(delay + j);
    
    for (i0 = 0; i0 < num; i0++) {
      for (i = 0; i < num; i++) {
	for (j = 0; j < num; j++) {
	  spp = *(sp + (i * num) + i0) + *(sp + (i0 * num) + j);
	  if (*(sp + (i * num) + j) > spp) *(sp + (i * num) + j) = spp;
	};
      };
    };

    double ub_now = DBL_MAX;
    int destn = *gr_matrix;
    int pos = 0;
    double *node_delay = new double[num];
    double *perm_delay = new double[num];
    int *tree = new int[num];
    for (i = 0; i < num; i++) {
      *(tree + i) = not_in_tree;
      *(node_delay + i) = 0.0;
    };
    node_delay[destn] = 0.0;
    int nodes_visited = 0;
    int *from = new int[num];
    int *to = new int[num];
    for (i = 0; i < num; i++) *(from + i) = *(to + i) = INT_MAX;
    bandb(alg, source->name(), destn, destn, pos, adj, node_delay, 
	  perm_delay, deg, cost, delay, head, tail, sp, nodes_visited,
	  ub_now, grCount, gr_matrix, tree, from, to, links, child,
	  minwt_at_node);
    delete [] gr_matrix;
    delete [] head;
    delete [] tail;
    delete [] deg;
    delete [] cost;
    delete [] delay;
    delete [] sp;
    delete [] child;
    delete [] minwt_at_node;
    for (i = 0; i < num; i++) delete [] *(adj + i);
    delete [] adj;
    delete [] node_delay;
    delete [] perm_delay;
    delete [] tree;
    if (ub_now == DBL_MAX) {
      delete [] from;
      delete [] to;
      return(SATORDB);
    }
    else {
      //Create a routing table entry for each MC group member;
      tmp = group->headm();
      while (tmp != NULL) {
	tmp->nodePtr()->addRoutingEntry(addr, source);
	tmp = tmp->next();
      };

      //Create the actual MC tree;
      i = 0;
      nodes = 1;
      while (*(from + i) != INT_MAX) {
	nodes++;
	Node *nd = nodeOf(*(to + i));
	Node *nd2 = nodeOf(*(from + i));
	AdjacencyListEntry *adj = nd2->adjacentNodes();
	while (adj->nodePtr() != nd) adj = adj->next();
	//update the link cost;
	double wght = adj->peak() + pk;
	adj->peak(wght);
	double average = adj->average() + avg;
	adj->average(average);

	//and add it to the routing table of the best connection node;
	nd2->addChild(addr, source, nd);
	i++;
      };
      delete [] from;
      delete [] to;

      maxd = 0;
      mind = DBL_MAX;
      int dbViolation = False;
      //calculate the expected end-to-end delay and the average number of hops;
      //and the cost per destination;
      NodeListEntry *tmp2 = group->headm();
      double avgDelay = 0;
      double avgHops = 0;
      while (tmp2 != NULL) {
	int gotit = False;
	int hops = 0;
	double delay = 0;
	if (source != tmp2->nodePtr()) 
	  results(source, addr, source, tmp2->nodePtr(), delay, hops, gotit);   
	if (delay >= DELAYBOUND) dbViolation = True;
	if (delay > maxd) maxd = delay;
	if (delay < mind) mind = delay;
	avgDelay += delay;
	avgHops += hops;
	tmp2 = tmp2->next();
      };
      avgDelay /= group->count();
      avgHops /= group->count();
      d = avgDelay;
      h = avgHops;
      
      if ((DBV == True) && (dbViolation == True)) {
	removeTree(source, addr);
	return(DBVIOL);
      }
      else return(ub_now);
    };
  }
  else return(NOGROUP);
};

void TheNodeList::sort(int node, int degree, double *cost, int *adj) {

  int i, j, temp;
  for (i = 1; i < degree; i++) {
    for (j = i; j != 0 ; j--) {
      if (*(cost + *(adj + j)) < *(cost + *(adj + j - 1))) {
	temp = *(adj + j);
	*(adj + j) = *(adj + j - 1);
	*(adj + j - 1) = temp;
      };
    };
  };
};

void TheNodeList::bandb(int alg, int srce, int destn, 
			int node, int pos, int **adj,
			double *node_delay, double *perm_delay, int *deg, 
			double *cost, double *delay, int *head, int *tail, 
			double *sp, int &nodes_visited, double &ub_now,
			int grCount, int *gr_matrix, int *tree,
			int *from, int *to, int links, int *child,
			double *minwt_at_node) {

  int arc,count,child_v,r,r1,r2,vertx,j,deg_d,parent,t_parent;
  double lb_here,delay_temp1,delay_temp2;

  if(node == srce) *(perm_delay + srce) = 0.0;
  if (pos == 0) *(tree + srce) = not_in_tree;
  deg_d = *(deg + node);

  if ((node != srce) && (*(tree + node) == not_in_tree)) {
    nodes_visited++;
    for (r = 0; r < deg_d; r++) {
      arc = *(*(adj + node) + r);
      *(tree + node) = temp_in_tree;
      parent = *(tail + arc);

      t_parent = *(tree + parent);
      //Find the reverse arc
      lb_here = lb(pos, grCount, head, tail, child, tree, cost, 
		   minwt_at_node, gr_matrix, links) + *(cost + arc);

      if (lb_here < ub_now) {  

	delay_temp1 = *(node_delay + node) + *(delay + arc);
	if (t_parent == perm_in_tree) 
	  delay_temp2 = *(perm_delay + parent) + delay_temp1;
	if (t_parent == not_in_tree) 
	  delay_temp2 = *(sp + (srce * num) + parent) + delay_temp1;

	r1 = 2;
	if (alg == copt) {
	  if (t_parent==not_in_tree && delay_temp2<DELAYBOUND) r1=1;
	  if (t_parent==not_in_tree && delay_temp2>=DELAYBOUND) r1=2;
	  if (t_parent==temp_in_tree) r1=2;
	  if (t_parent==perm_in_tree && delay_temp2<DELAYBOUND) r1=3;
	  if (t_parent==perm_in_tree && delay_temp2>=DELAYBOUND) r1=2;
	}
	else {
	  if (t_parent==not_in_tree) r1=1;
	  if (t_parent==temp_in_tree) r1=2;
	  if (t_parent==perm_in_tree) r1=3;
	};

	switch(r1){
	case 1: /*  incoming_edge = arc  */
	  *(node_delay + parent) = delay_temp1;    
	  *(child + (parent * grCount) + pos) = node;
	  bandb(alg, srce, destn, parent, pos, adj, node_delay, 
		perm_delay, deg, cost, delay, head, tail, sp, nodes_visited,
		ub_now, grCount, gr_matrix, tree, from, to, links,
		child, minwt_at_node);
	  *(child + (parent * grCount)  + pos) = INT_MAX;
	  break;
	case 2: break;
	case 3: /*  incoming_edge = arc  */
	  *(child + (parent * grCount) + pos) = node;
	  bandb(alg, srce, destn, parent, pos, adj, node_delay, 
		perm_delay, deg, cost, delay, head, tail, sp, nodes_visited,
		ub_now, grCount, gr_matrix, tree, from, to, links,
		child, minwt_at_node);
	  *(child + (parent * grCount)  + pos) = INT_MAX;
	  break;
	}; /* end switch */
      }  /* end if */
      else break;  	/*  break from for loop  */
    }; /* end for */
    *(tree + node) = not_in_tree;
    *(node_delay + node) = 0.0; 
  } /* end if */
  r2 = 3;
  if ((*(tree + node) == perm_in_tree) && (node == destn)) r2 = 1;
  if (((node == srce) || (*(tree + node) == perm_in_tree)) &&
      (node != destn)) r2 = 2;

  switch(r2){
  case 2:
    *(tree + node) = perm_in_tree;
    vertx = node;
    
    while (*(child + (vertx * grCount) + pos) != INT_MAX) {
      child_v = *(child + (vertx * grCount) + pos);
      deg_d = *(deg + child_v);
      for (r1 = 0; r1 < deg_d; r1++)
	if (*(tail + *(*(adj + child_v) + r1)) == vertx) 
	  j = *(*(adj + child_v) + r1);
      *(perm_delay + child_v) = *(perm_delay + vertx) + *(delay + j);
      vertx = child_v;
      *(tree + vertx) = perm_in_tree;
    };

    if (pos < (grCount - 1)) {
      *(child + (destn * grCount) + pos) = INT_MAX;
      pos++;
      destn = *(gr_matrix + pos);
      bandb(alg, srce, destn, destn, pos, adj, node_delay, 
	    perm_delay, deg, cost, delay, head, tail, sp, nodes_visited,
	    ub_now, grCount, gr_matrix, tree, from, to, links,
	    child, minwt_at_node);
      pos--; 
      destn = *(gr_matrix + pos);
    }
    else {
      lb_here = lb(pos, grCount, head, tail, child, tree, cost, 
		   minwt_at_node, gr_matrix, links);
      if (lb_here < ub_now) {
	ub_now = lb_here;  /* update upper bd */
	int i, k = 0;
	for (j = 0; j < grCount; j++) {
	  for (i = 0; i < num; i++) {
	    if (*(child + (i * grCount) + j) != INT_MAX) {
	      *(from + k) = i;
	      *(to +k++) = *(child + (i * grCount) + j);
	    };
	  };
	};
	*(from + k) = *(to + k) = INT_MAX;
      };
    };
    vertx = node;
    /* count = 0; */
    while (*(child + (vertx * grCount) + pos) != INT_MAX) {
      child_v = *(child + (vertx * grCount) + pos);
      *(tree + child_v) = temp_in_tree;
      vertx = child_v;
    };
    break;
    
  case 1:
    if (pos < (grCount - 1)) {
      *(child + (destn * grCount) + pos) = INT_MAX;
      pos++;
      destn = *(gr_matrix + pos);
      bandb(alg, srce, destn, destn, pos, adj, node_delay, 
	    perm_delay, deg, cost, delay, head, tail, sp, nodes_visited,
	    ub_now, grCount, gr_matrix, tree, from, to, links,
	    child, minwt_at_node);
      pos--; 
      destn = *(gr_matrix + pos);
    }
    else {
      lb_here = lb(pos, grCount, head, tail, child, tree, cost, 
		   minwt_at_node, gr_matrix, links);
      if (lb_here < ub_now) {
	ub_now = lb_here;  /* update upper bd */
	ub_now = lb_here;  /* update upper bd */
	int i, k = 0;
	for (j = 0; j < grCount; j++) {
	  for (i = 0; i < num; i++) {
	    if (*(child + (i * grCount) + j) != INT_MAX) {
	      *(from + k) = i;
	      *(to + k++) = *(child + (i * grCount) + j);
	    };
	  };
	};
	*(from + k) = *(to + k) = INT_MAX;
      };
    };
    break;

  case 3: break;
  };
};

double TheNodeList::lb(int curr_pos, int grCount, int *head, int *tail,
		       int *child, int *tree, double *cost,
		       double *minwt_at_node, int *gr_matrix, int links) {
  int tj, hj, j, pos;
  double lbd;
  
  lbd = 0.0;
  for (j = 0; j < links; j++) {
    hj = *(head + j);
    tj = *(tail + j);
    int child_hj_tj = False;
    for (pos = 0; pos < grCount; pos++) {
      if (*(child + (tj * grCount) + pos) == hj)
	child_hj_tj = True;
    };
    if ((*(tree + hj) == temp_in_tree) && (*(tree + tj) == temp_in_tree) &&
	(child_hj_tj == True)) lbd += *(cost + j);
    if ((*(tree + hj) == perm_in_tree) && (*(tree + tj) == perm_in_tree) &&
	(child_hj_tj == True)) lbd += *(cost + j);
  }; /* end for */

  for (pos = curr_pos; pos <  grCount; pos++)
    if (*(tree + *(gr_matrix + pos)) == not_in_tree) 
      lbd += *(minwt_at_node + *(gr_matrix + pos));
  return lbd;
};