/*****************************************************************/
/*                                                               */
/* Program to Solve the                                          */
/* Multi-level Network Optimization Problem                      */
/*                                                               */
/* Branch and Bound Retated Routines                             */
/*                                                               */
/* Author: Frederico R. B. da Cruz                               */
/*         Department of Industrial Engineering and              */
/*         Operations Research                                   */
/*         University of Massachusetts at Amherst                */
/*         USA                                                   */
/*                                                               */
/* E-mail: fcruz@sml1.ecs.umass.edu                              */
/*                                                               */
/* Version: 3.00                                                 */
/*                                                               */
/* Date: May/98                                                  */
/*                                                               */
/*****************************************************************/

#include <stdio.h>
#include <stdlib.h>

#include "lr.c"

void Branch_Bound    (int                n_Nodes,
                      int                n_Levels,
                      ListArcPtrType     *AdjaNodInd,
                      ListNodPtrType     *DemaNodLev,
                      ListNodPtrType     *RootNodLev,
                      MlType             *Ml,
                      slType             *sl,
                      PredLevPtrType     **BestPredLev,
                      PredLevPtrType     **GoodPredLev,
                      PredLevPtrType     **CurrPredLev,
                      float              *BestUpper,
                      float              *UpperBound,
                      float              *LowerBound,
                      int                *Visited);

void SetNode         (NodPtrType         NodPtr,
                      int                Value);

void SetArc          (ArcPtrType         ArcPtr,
                      int                Value);

void Reduce          (int                n_Nodes,
                      int                n_Levels,
                      ListArcPtrType     *AdjaNodInd,
                      ListNodPtrType     *DemaNodLev,
                      ListNodPtrType     *RootNodLev,
                      MlType             *Ml,
                      slType             *sl,
                      PredLevPtrType     **CurrPredLev,
                      float              *BestUpper,
                      ListArcPtrType     *ListArcFix,
                      ListNodPtrType     *ListNodFix);

void UnReduce        (ListArcPtrType     *ListArcFix,
                      ListNodPtrType     *ListNodFix);

int ChooseVar        (int                n_Nodes,
                      int                n_Levels,
                      ListArcPtrType     *AdjaNodInd,
                      ListNodPtrType     *DemaNodLev,
                      ListNodPtrType     *RootNodLev,
                      MlType             *Ml,
                      slType             *sl,
                      PredLevPtrType     **CurrPredLev,
                      float              *BestUpper,
                      void               **Chosen);



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


void Branch_Bound    (int                n_Nodes,
                      int                n_Levels,
                      ListArcPtrType     *AdjaNodInd,
                      ListNodPtrType     *DemaNodLev,
                      ListNodPtrType     *RootNodLev,
                      MlType             *Ml,
                      slType             *sl,
                      PredLevPtrType     **BestPredLev,
                      PredLevPtrType     **GoodPredLev,
                      PredLevPtrType     **CurrPredLev,
                      float              *BestUpper,
                      float              *UpperBound,
                      float              *LowerBound,
                      int                *Visited)
{
  float             Gap;
  PredLevPtrType    *AuxPredLev;
  ArcPtrType        ArcPtr;
  NodPtrType        NodPtr;
  int               Variable;
  void              *Chosen;

  /*   interrupt explorations?  */
#if   _B_B_NODES_ == _ONE_
  if ( *Visited >= 1 ) return;
#elif _B_B_NODES_ == _TEN_THOUSAND_
  if ( *Visited >= 10000 ) return;
#elif _B_B_NODES_ == _FIFTY_THOUSAND_
  if ( *Visited >= 50000 ) return;
#else
  static float      MaxTime = 8000.0;
  static int        FirstTime = YES;
  if ((*Visited%100==0)&&(ElapsedTime()>MaxTime)) {
    if (FirstTime==YES) {
      FirstTime=NO;
      fprintf(OUTFILE, "\nWARNING: Time Overflow (%1.0f seconds)\n", MaxTime);
    }
    return;
  }
#endif

  /*   bound  */
  (*Visited)++;
#if _LR_USE_ == _ONCE_
  if ( *Visited == 1 )
     Lagr_Relax          (n_Nodes, n_Levels,
                          AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                          GoodPredLev, CurrPredLev,
                          BestUpper, UpperBound, LowerBound);
  else
     Bounds              (n_Nodes, n_Levels,
                          AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                          GoodPredLev, CurrPredLev,
                          UpperBound, LowerBound);
#elif _LR_USE_ == _ALL_THE_TIME_
     Lagr_Relax          (n_Nodes, n_Levels,
                          AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                          GoodPredLev, CurrPredLev,
                          BestUpper, UpperBound, LowerBound);
#endif

  /*  update best upper bound  */
  if ( *UpperBound < *BestUpper ) {
     *BestUpper   = *UpperBound;
     AuxPredLev   = *BestPredLev;
     *BestPredLev = *GoodPredLev;
     *GoodPredLev = AuxPredLev;
  }

  /*  compute gap  */
/*
  if ( *LowerBound > 0 )
     Gap = *UpperBound / *LowerBound - 1;
  else if ( *LowerBound < 0 )
     Gap = -( *UpperBound / *LowerBound - 1);
  else
     Gap = INFINITY;
*/
  Gap = 1 - ( *LowerBound / *UpperBound );

  /*  print results */
  if ( *Visited == 1 ) {
     fprintf(OUTFILE, "  (First Node)");
     fprintf(OUTFILE, "\n");
     fprintf(OUTFILE, "Upper bound  %f    \n", *UpperBound);
     fprintf(OUTFILE, "Lower bound  %f    \n", *LowerBound);
     fprintf(OUTFILE, "Gap          %f %% \n", 100 * Gap);
     fprintf(OUTFILE, "Time (watch) %10.5f\n", ElapsedTime());
     StartCronograph();
/*
*/
  }
/*
  else {
     fprintf(OUTFILE, "\n");
     fprintf(OUTFILE, "Upper bound  %f    \n", *UpperBound);
     fprintf(OUTFILE, "Lower bound  %f    \n", *LowerBound);
     fprintf(OUTFILE, "Gap          %f %% \n", 100 * Gap);
  }
*/

  /* branch if necessary */
  if ((*LowerBound >= INFINITY) || (*BestUpper / *LowerBound - 1 < -ZERO)) {
/* 
     fprintf(OUTFILE,"Infeasible node reached\n");
*/
  }
  else if (Gap < EPSILON) {
/* 
     fprintf(OUTFILE, "An optimum solution reached\n"); 
*/
  }
  else {
/* 
     fprintf(OUTFILE, "Uh-oh! It is really necessary to branch\n"); 
*/
#if _REDUCTION_ == _IN_USE_
     ListArcPtrType    ListArcFix;
     ListNodPtrType    ListNodFix;

     Reduce(n_Nodes, n_Levels,
            AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
            CurrPredLev,
            BestUpper,
            &ListArcFix,
            &ListNodFix);
#endif
     Variable = ChooseVar(n_Nodes, n_Levels,
                          AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                          CurrPredLev, BestUpper, &Chosen);
     if ( Variable == IS_NODE ) {
        NodPtr = (NodPtrType) Chosen;
        SetNode( NodPtr,  SET_TO_1 );
        Branch_Bound(n_Nodes, n_Levels,
                     AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                     BestPredLev, GoodPredLev, CurrPredLev,
                     BestUpper, UpperBound, LowerBound, Visited);
        SetNode( NodPtr,  SET_TO_0 );
        Branch_Bound(n_Nodes, n_Levels,
                     AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                     BestPredLev, GoodPredLev, CurrPredLev,
                     BestUpper, UpperBound, LowerBound, Visited);
        SetNode( NodPtr,  SET_FREE );
     }
     else if ( Variable == IS_ARC ) {
        ArcPtr = (ArcPtrType) Chosen;
        SetArc( ArcPtr, SET_TO_1 );
        Branch_Bound(n_Nodes, n_Levels,
                     AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                     BestPredLev, GoodPredLev, CurrPredLev,
                     BestUpper, UpperBound, LowerBound, Visited);
        SetArc( ArcPtr, SET_TO_0 );
        Branch_Bound(n_Nodes, n_Levels,
                     AdjaNodInd, DemaNodLev, RootNodLev, Ml, sl,
                     BestPredLev, GoodPredLev, CurrPredLev,
                     BestUpper, UpperBound, LowerBound, Visited);
        SetArc( ArcPtr, SET_FREE );
     }
     else {
/* 
        fprintf(OUTFILE, "No more free variables\n"); 
*/
     }
#if _REDUCTION_ == _IN_USE_
     UnReduce(&ListArcFix, &ListNodFix);
#endif
  }
}

void SetNode         (NodPtrType         NodPtr,
                      int                Value)
{
  NodPtr->Status = Value;
/*
  if ( Value == SET_FREE )
      fprintf(OUTFILE, "\nNode (%d) set free\n", NodPtr->i+1);
  else
      fprintf(OUTFILE, "\nNode (%d) set to %d\n", NodPtr->i+1, Value);
*/
}

void SetArc          (ArcPtrType         ArcPtr,
                      int                Value)
{
  ArcPtr->Status = Value;
/*
  if ( Value == SET_FREE )
      fprintf(OUTFILE, "\nArc (%d,%d) set free\n", ArcPtr->i+1, ArcPtr->j+1 );
  else
      fprintf(OUTFILE, "\nArc (%d,%d) set to %d\n", 
              ArcPtr->i+1, ArcPtr->j+1, Value);
*/
}

void Reduce          (int                n_Nodes,
                      int                n_Levels,
                      ListArcPtrType     *AdjaNodInd,
                      ListNodPtrType     *DemaNodLev,
                      ListNodPtrType     *RootNodLev,
                      MlType             *Ml,
                      slType             *sl,
                      PredLevPtrType     **CurrPredLev,
                      float              *BestUpper,
                      ListArcPtrType     *ListArcFix,
                      ListNodPtrType     *ListNodFix)
{
  ListArcPtrType    ListArcPtr, ListArcAux;
  ListNodPtrType    ListNodPtr, ListNodAux;
/*
  ArcPtrType        *ArcPtr;
  NodPtrType        NodPtr;
*/
  float             Solution_L;
  int               l, i;
  int               Reduced;
  
  (*ListArcFix) = (ListArcPtrType) malloc( sizeof(ListArcType) );
  (*ListArcFix)->APtr = (ArcPtrType *) calloc(1,sizeof(ArcPtrType));
  (*ListArcFix)->APtr = NULL;
  (*ListArcFix)->Next = NULL;
  (*ListNodFix) = (ListNodPtrType) malloc( sizeof(ListNodType) );
  (*ListNodFix)->NPtr = NULL;
  (*ListNodFix)->Next = NULL;

  do {
    Reduced = NO;
    for ( l = 0; l < n_Levels; l++ ) {
      /*                    */
      /*  working on nodes  */
      /*                    */
      ListNodPtr = RootNodLev[l];
      while ( ListNodPtr != NULL ) {
         if ( ListNodPtr->NPtr->Status == SET_FREE ) {
/*
            fprintf(OUTFILE, "\nNow working on node %2d\n",
                    ListNodPtr->NPtr->i+1 );
*/
            SetNode( ListNodPtr->NPtr, SET_TO_1 );
            Solution_L = Solve_L(n_Nodes, n_Levels,
                                 AdjaNodInd, DemaNodLev, RootNodLev,
                                 *CurrPredLev, Ml, sl); 
            if ( Solution_L > *BestUpper ) {
               Reduced = YES;
               SetNode(  ListNodPtr->NPtr, SET_TO_0 );
               /*  annotate node just fixed  */
               ListNodAux = (ListNodPtrType) malloc(sizeof(ListNodType));
               ListNodAux->NPtr = ListNodPtr->NPtr;
               ListNodAux->Next = (*ListNodFix)->Next;
               (*ListNodFix)->Next = ListNodAux;
            }
            else {
               SetNode(  ListNodPtr->NPtr, SET_TO_0 );
               Solution_L = Solve_L(n_Nodes, n_Levels,
                                    AdjaNodInd, DemaNodLev, RootNodLev,
                                    *CurrPredLev, Ml, sl); 
               if ( Solution_L > *BestUpper ) {
                  Reduced = YES;
                  SetNode(  ListNodPtr->NPtr, SET_TO_1 );
                  /*  annotate node just fixed  */
                  ListNodAux = (ListNodPtrType) malloc(sizeof(ListNodType));
                  ListNodAux->NPtr = ListNodPtr->NPtr;
                  ListNodAux->Next = (*ListNodFix)->Next;
                  (*ListNodFix)->Next = ListNodAux;
               }
               else {
                  SetNode(  ListNodPtr->NPtr, SET_FREE );
               }
            }
         }
         ListNodPtr = ListNodPtr->Next;
      }
      /*                   */
      /*  working on arcs  */
      /*                   */
      for ( i = 0; i < n_Nodes; i++ ) {
         ListArcPtr = AdjaNodInd[i];
         while ( ListArcPtr != NULL ) {
            if ( ListArcPtr->APtr[l]->Status == SET_FREE ) {
/*
               fprintf(OUTFILE, "\nNow working on arc (%2d,%2d)^%d\n",
                       ListArcPtr->APtr[l]->i+1,
                       ListArcPtr->APtr[l]->j+1,
                       l+1 );
*/
               SetArc(  ListArcPtr->APtr[l], SET_TO_1 );
               Solution_L = Solve_L(n_Nodes, n_Levels,
                                    AdjaNodInd, DemaNodLev, RootNodLev,
                                    *CurrPredLev, Ml, sl); 
               if ( Solution_L > *BestUpper ) {
                  SetArc(  ListArcPtr->APtr[l], SET_TO_0 );
                  Reduced = YES;
                  /*  annotate arcs just fixed  */
                  ListArcAux = (ListArcPtrType) malloc(sizeof(ListArcType));
                  ListArcAux->APtr=(ArcPtrType *)calloc(1,sizeof(ArcPtrType));
                  ListArcAux->APtr[0] = ListArcPtr->APtr[l];
                  ListArcAux->Next = (*ListArcFix)->Next;
                  (*ListArcFix)->Next = ListArcAux;
               }
               else {
                  SetArc(  ListArcPtr->APtr[l], SET_TO_0 );
                  Solution_L = Solve_L(n_Nodes, n_Levels,
                                       AdjaNodInd, DemaNodLev, RootNodLev,
                                       *CurrPredLev, Ml, sl); 
                  if ( Solution_L > *BestUpper ) {
                      SetArc(  ListArcPtr->APtr[l], SET_TO_1 );
                      Reduced = YES;
                      /*  annotate arcs just fixed  */
                      ListArcAux = (ListArcPtrType) malloc(sizeof(ListArcType));
                      ListArcAux->APtr = (ArcPtrType *)
                                         calloc(1,sizeof(ArcPtrType));
                      ListArcAux->APtr[0] = ListArcPtr->APtr[l];
                      ListArcAux->Next = (*ListArcFix)->Next;
                      (*ListArcFix)->Next = ListArcAux;
                   }
                   else {
                      SetArc(  ListArcPtr->APtr[l], SET_FREE );
                   }
               }
            }
            ListArcPtr = ListArcPtr->Next;
        }    
      }
    }
  } while ( Reduced == YES );
/* 
  do { 
  } while (0);
  do {
  } while ( Reduced == YES );
*/

  /* testing routine */
/*
  fprintf(OUTFILE, "\nNew nodes/arcs in K0 = {" );
  ListNodPtr = (*ListNodFix)->Next;
  while ( ListNodPtr != NULL ) {
      if ( ListNodPtr->NPtr->Status == SET_TO_0 )
         fprintf(OUTFILE, " (%2d)",
                 ListNodPtr->NPtr->i+1  );
      ListNodPtr = ListNodPtr->Next;
  }
  ListArcPtr = (*ListArcFix)->Next;
  while ( ListArcPtr != NULL ) {
      if ( ListArcPtr->APtr[0]->Status == SET_TO_0 )
         fprintf(OUTFILE, " (%2d,%2d)",
                          ListArcPtr->APtr[0]->i+1,
                          ListArcPtr->APtr[0]->j+1  );
      ListArcPtr = ListArcPtr->Next;
  }
  fprintf(OUTFILE, " }\n");

  fprintf(OUTFILE, "\nNew nodes/arcs in K1 = {" );
  ListNodPtr = (*ListNodFix)->Next;
  while ( ListNodPtr != NULL ) {
      if ( ListNodPtr->NPtr->Status == SET_TO_1 )
         fprintf(OUTFILE, " (%2d)",
                 ListNodPtr->NPtr->i+1  );
      ListNodPtr = ListNodPtr->Next;
  }
  ListArcPtr = (*ListArcFix)->Next;
  while ( ListArcPtr != NULL ) {
      if ( ListArcPtr->APtr[0]->Status == SET_TO_1 )
         fprintf(OUTFILE, " (%2d,%2d)",
                          ListArcPtr->APtr[0]->i+1,
                          ListArcPtr->APtr[0]->j+1  );
      ListArcPtr = ListArcPtr->Next;
  }
  fprintf(OUTFILE, " }\n");

  for ( l = 0; l < n_Levels; l++ ) {
    fprintf(OUTFILE, "\nK0^%d = {", l+1 );
    ListNodPtr = RootNodLev[l];
    while ( ListNodPtr != NULL ) {
        if ( ListNodPtr->NPtr->Status == SET_TO_0 )
           fprintf(OUTFILE, " (%2d)",
                   ListNodPtr->NPtr->i+1  );
        ListNodPtr = ListNodPtr->Next;
    }
    for ( i = 0; i < n_Nodes; i++ ) {
        ListArcPtr = AdjaNodInd[i];
        while ( ListArcPtr != NULL ) {
          if ( ListArcPtr->APtr[l]->Status == SET_TO_0 )
             fprintf(OUTFILE, " (%2d,%2d)",
                              ListArcPtr->APtr[l]->i+1,
                              ListArcPtr->APtr[l]->j+1  );
          ListArcPtr = ListArcPtr->Next;
        }
    }
    fprintf(OUTFILE, " }\n");
  }

  for ( l = 0; l < n_Levels; l++ ) {
    fprintf(OUTFILE, "\nK1^%d = { ", l+1 );
    ListNodPtr = RootNodLev[l];
    while ( ListNodPtr != NULL ) {
        if ( ListNodPtr->NPtr->Status == SET_TO_1 )
           fprintf(OUTFILE, " (%2d)",
                   ListNodPtr->NPtr->i+1  );
        ListNodPtr = ListNodPtr->Next;
    }
    for ( i = 0; i < n_Nodes; i++ ) {
        ListArcPtr = AdjaNodInd[i];
        while ( ListArcPtr != NULL ) {
          if ( ListArcPtr->APtr[l]->Status == SET_TO_1 )
             fprintf(OUTFILE, " (%2d,%2d)",
                              ListArcPtr->APtr[l]->i+1,
                              ListArcPtr->APtr[l]->j+1  );
          ListArcPtr = ListArcPtr->Next;
        }
    }
    fprintf(OUTFILE, " }\n");
  }
*/
}

void UnReduce        (ListArcPtrType     *ListArcFix,
                      ListNodPtrType     *ListNodFix)
{
  ListArcPtrType  ListArcPtr, ListArcAux;
  ListNodPtrType  ListNodPtr, ListNodAux;

  ListArcPtr = (*ListArcFix)->Next;
  free( (*ListArcFix)->APtr );
  free( (*ListArcFix) );
  while ( ListArcPtr != NULL ) {
        ListArcAux = ListArcPtr;
        ListArcPtr = ListArcPtr->Next;

        SetArc(  ListArcAux->APtr[0], SET_FREE );
        free( ListArcAux->APtr );
        free( ListArcAux );
  }
  ListNodPtr = (*ListNodFix)->Next;
  free( (*ListNodFix) );
  while ( ListNodPtr != NULL ) {
        ListNodAux = ListNodPtr;
        ListNodPtr = ListNodPtr->Next;

        SetNode(  ListNodAux->NPtr, SET_FREE );
        free( ListNodAux );
  }
}

int ChooseVar        (int                n_Nodes,
                      int                n_Levels,
                      ListArcPtrType     *AdjaNodInd,
                      ListNodPtrType     *DemaNodLev,
                      ListNodPtrType     *RootNodLev,
                      MlType             *Ml,
                      slType             *sl,
                      PredLevPtrType     **CurrPredLev,
                      float              *BestUpper,
                      void               **Chosen)
{
  int            i,l;
  ListNodPtrType ListNodPtr;
  ListArcPtrType ListArcPtr;
  NodPtrType     NodPtr;
  ArcPtrType     ArcPtr;
  int            ArcFound;
  int            NodFound;

#if _CHOICE_STRATEGY_ == _FIRST_FREE_

  ArcFound = NO;
  NodFound = NO;
  l = 0;
  while ( ( l < n_Levels ) && ( ArcFound == NO ) && ( NodFound == NO ) ) {
      /*  try to choose a node  */
     ListNodPtr = RootNodLev[l];
     while ( ( ListNodPtr != NULL ) && 
             ( ListNodPtr->NPtr->Status != SET_FREE ) )
           ListNodPtr = ListNodPtr->Next;
     if ( ListNodPtr != NULL ) {
        (*Chosen) = (void *) ListNodPtr->NPtr;
        NodFound = YES;
     }
      /*  otherwise try to choose an arc  */
     else {
         i = 0;
         while ( ( i < n_Nodes ) && ( ArcFound == NO ) ) {
            ListArcPtr = AdjaNodInd[i];
            while ( ( ListArcPtr != NULL ) && 
                    ( ListArcPtr->APtr[l]->Status != SET_FREE ) )
                  ListArcPtr = ListArcPtr->Next;
            if ( ListArcPtr != NULL ) {
               (*Chosen) = (void *) ListArcPtr->APtr[l];
               ArcFound = YES;
            }
            i++;
         }
     }
     l++;
  }
#elif _CHOICE_STRATEGY_ == _GOOD_INCREASE_
  float          Delta, Max_Delta_y_ij, Max_Delta_z_k;;

  ArcFound = NO;
  NodFound = NO;
  /* compute Max_Delta_z_k */
  Max_Delta_z_k = - INFINITY;
  for (l = 0; l < n_Levels; l++ ) {
    ListNodPtr = RootNodLev[l];
    while (ListNodPtr != NULL) {
      if ( ListNodPtr->NPtr->Status == SET_FREE ) {
        Delta = ListNodPtr->NPtr->vi * ListNodPtr->NPtr->Gamma_vi;
/*
        fprintf(OUTFILE," vi * Gamma = %f * %f = %f\n",
                ListNodPtr->NPtr->vi, ListNodPtr->NPtr->Gamma_vi, Delta);
*/
        if ( Max_Delta_z_k < Delta ) {
          Max_Delta_z_k = Delta;
          NodPtr = ListNodPtr->NPtr;
        }
      }
      ListNodPtr = ListNodPtr->Next;
    }
  }

  /* compute Max_Delta_y_ij */
  Max_Delta_y_ij = -INFINITY;
  for (i = 0; i < n_Nodes; i++ ) {
    ListArcPtr = AdjaNodInd[i];
    while (ListArcPtr != NULL) {
      for (l = 0; l < n_Levels; l++) {
        if ( ListArcPtr->APtr[l]->Status == SET_FREE ) {
          Delta = ListArcPtr->APtr[l]->wij * ListArcPtr->APtr[l]->Gamma_wij;
/*
          fprintf(OUTFILE,"wij * Gamma = %f * %f = %f\n",
                  ListArcPtr->APtr[l]->wij, 
                  ListArcPtr->APtr[l]->Gamma_wij, Delta);
*/
          if ( Max_Delta_y_ij < Delta ) {
            Max_Delta_y_ij = Delta;
            ArcPtr = ListArcPtr->APtr[l];
          }
        }
      }
      ListArcPtr = ListArcPtr->Next;
    }
  }

  /* chose the maximum */
  if ( (Max_Delta_y_ij > -INFINITY) || (Max_Delta_z_k > -INFINITY) ) {
    if (Max_Delta_y_ij > Max_Delta_z_k) {
       ArcFound = YES;
       (*Chosen) = (void *) ArcPtr;
     }
     else {
       NodFound = YES;
       (*Chosen) = (void *) NodPtr;
     }
  }
#elif _CHOICE_STRATEGY_ == _NO_CYCLES_
  int          *Reached, Disconnected;

  ArcFound = NO;
  NodFound = NO;
  l = 0;
   /*  try to choose a first-level node */
  ListNodPtr = RootNodLev[l];
  while ( ( ListNodPtr != NULL ) &&
          ( ListNodPtr->NPtr->Status != SET_FREE ) )
       ListNodPtr = ListNodPtr->Next;
  if ( ListNodPtr != NULL ) {
     (*Chosen) = (void *) ListNodPtr->NPtr;
     NodFound = YES;
  }
  /* otherwise go throughout all levels */
  while ( ( l < n_Levels ) && ( ArcFound == NO ) && ( NodFound == NO ) ) {
      /*  try to choose a (l+1)th-level node */
     ListNodPtr = RootNodLev[l+1];
     while ( ( ListNodPtr != NULL ) &&
             ( ListNodPtr->NPtr->Status != SET_FREE ) )
          ListNodPtr = ListNodPtr->Next;
     if ( ListNodPtr != NULL ) {
        (*Chosen) = (void *) ListNodPtr->NPtr;
        NodFound = YES;
     }
      /*  otherwise try to choose a l-th-level arc  */
     else {
          /*  compute set of reached nodes  */
              /*  initialize  */
         Reached = (int *) calloc ( n_Nodes, sizeof(int) );
         for ( i = 0; i < n_Nodes; i++ ) 
             Reached[i] = NO;
             /*  nodes definetly reached  */
         ListNodPtr = RootNodLev[l];
         while ( ListNodPtr != NULL ) {
            if ( ListNodPtr->NPtr->Status == SET_TO_1 )
               Reached[ListNodPtr->NPtr->i] = YES;
            ListNodPtr = ListNodPtr->Next;
         }
             /*  others  */
         for ( i = 0; i < n_Nodes; i++ ) {
             ListArcPtr = AdjaNodInd[i];
             while ( ListArcPtr != NULL ) {
                   if ( ListArcPtr->APtr[l]->Status == SET_TO_1 ) {
                      Reached[ListArcPtr->APtr[l]->i] = YES;
                      Reached[ListArcPtr->APtr[l]->j] = YES;
                   }
                   ListArcPtr = ListArcPtr->Next;
             }
         }
          /*  is there a not reached node in D^l \cup (R^{l+1} \cap J_1)?  */
         Disconnected = NO;
         ListNodPtr = DemaNodLev[l];
         while ( (ListNodPtr != NULL ) && (Disconnected == NO) ) {
               if ( Reached[ListNodPtr->NPtr->i] == NO )
                  Disconnected = YES;
               ListNodPtr = ListNodPtr->Next;
         }
         ListNodPtr = RootNodLev[l+1];
         while ( (ListNodPtr != NULL ) && (Disconnected == NO) ) {
               if ( ( ListNodPtr->NPtr->Status == SET_TO_1 ) &&
                    ( Reached[ListNodPtr->NPtr->i] == NO ) )
                  Disconnected = YES;
               ListNodPtr = ListNodPtr->Next;
         }
          /*  if yes choose a free arc emanating from a reached   */
          /*  node whose other extremity has not been reached yet */
         if ( Disconnected == YES ) {
            i = 0;
            while ( ( i < n_Nodes ) && ( ArcFound == NO ) ) {
               ListArcPtr = AdjaNodInd[i];
               while ( ( ListArcPtr != NULL ) && ( ArcFound == NO ) ) {
                     if ( ( ListArcPtr->APtr[l]->Status == SET_FREE ) &&
                          ( Reached[ListArcPtr->APtr[l]->i] == YES )  &&
                          ( Reached[ListArcPtr->APtr[l]->j] == NO  )     ) {
                          (*Chosen) = (void *) ListArcPtr->APtr[l];
                          ArcFound = YES;
                     }
                     ListArcPtr = ListArcPtr->Next;
               }
               i++;
            }
         }
          /*  free set of reached nodes  */
         free( Reached );
     }
     l++;
  }
#endif
  if ( NodFound == YES ) {
    NodPtr = (NodPtrType)(*Chosen);
/*
    fprintf(OUTFILE,"\nNode %d chosen\n", NodPtr->i+1);
*/
    return(IS_NODE);
  }
  else if ( ArcFound == YES ) {
    ArcPtr = (ArcPtrType)(*Chosen);
/*
    fprintf(OUTFILE,"\nArc (%d,%d) chosen\n", ArcPtr->i+1, ArcPtr->j+1);
*/
    return(IS_ARC);
  }
  else {
/*
    fprintf(OUTFILE,"\nNobody chosen\n");
*/
    return(NONE);
  }
}