adjacencywdigraph.h

function to compute an expression using int value parameters throw an exception of type illegalPara

               {
                  // distanceFromSource[j] decreases
                  distanceFromSource[j] = distanceFromSource[v] + a[v][j];
                  // add j to newReachableVertices
                  if (predecessor[j] == -1)
                     // not reached before
                     newReachableVertices.insert(0, j);
                  predecessor[j] = v;
               }
            }
         }
      }

      template <class T>
      void allPairs(T **c, int **kay)
      {// Dynamic programming all pairs shortest paths algorithm.
       // Compute c[i][j] and kay[i][j] for all i and j.
         if (!weighted())
            throw undefinedMethod
            ("adjacencyWDigraph::allPairs() not defined for unweighted graphs");
   
         // initialize c[i][j] = c(i,j,0)
         for (int i = 1; i <= n; i++)
            for (int j = 1; j <= n; j++)
            {
               c[i][j] = a[i][j];
               kay[i][j] = 0;
            }
         for (int i = 1; i <= n; i++)
            c[i][i] = 0;
         
         // compute c[i][j] = c(i,j,k)
         for (int k = 1; k <= n; k++)
            for (int i = 1; i <= n; i++)
               for (int j = 1; j <= n; j++)
                  if (c[i][k] != noEdge && c[k][j] != noEdge &&
                     (c[i][j] == noEdge || c[i][j] > c[i][k] + c[k][j]))
                  {// smaller value for c[i][j] found
                       c[i][j] = c[i][k] + c[k][j];
                       kay[i][j] = k;
                  }
      }

      T btSalesperson(int *bestTour)
      {// Traveling salesperson by backtracking.
       // bestTour[1:n] is set to best tour.
       // Return cost of best tour.
         if (!weighted())
            throw undefinedMethod
            ("adjacencyWDigraph::btSalesperson() not defined for unweighted graphs");
         // set partialTour to identity permutation
         partialTour = new int [n + 1];
         for (int i = 1; i <= n; i++)
            partialTour[i] = i;
   
         costOfBestTourSoFar = noEdge;
         bestTourSoFar = bestTour;
         costOfPartialTour = 0;
      
         // search permutations of partialTour[2:n]
         rTSP(2);
      
         return costOfBestTourSoFar;
      }

      // struct used by least-cost branch-and-bound traveling salesperson
      struct heapNode
      {
         // data members
         T lowerCost;            // lower bound on cost of tours in subtree
         T costOfPartialTour;    // cost of partial tour
         T minAdditionalCost;    // min additional cost to complete tour
         int sizeOfPartialTour;  // partial tour is
         int *partialTour;       // partialTour[sizeOfPartialTour+1:n-1]
                                 // gives remaining vertices to be added
                                 // to partialTour[0:sizeOfPartialTour]
         // constructors
         heapNode() {}

         heapNode(T lC, T cOPT, T mAC, int sOPT, int* pT)
         {
            lowerCost = lC;
            costOfPartialTour = cOPT;
            minAdditionalCost = mAC;
            sizeOfPartialTour = sOPT;
            partialTour = pT;
         }
   
         operator int() {return lowerCost;}

         operator>(const heapNode right)
            {return lowerCost > right.lowerCost;}
      };
      
      T leastCostBBSalesperson(int *bestTour)
      {// least-cost branch-and-bound code to find a shortest tour
       // bestTour[i] set to i'th vertex on shortest tour
       // Return cost of shortest tour.
         if (!weighted())
            throw undefinedMethod
            ("adjacencyWDigraph::leastCostBBSalesperson() not defined for unweighted graphs");

         minHeap<heapNode> liveNodeMinHeap;
      
         // costOfMinOutEdge[i] = cost of least-cost edge leaving vertex i
         T *costOfMinOutEdge = new T [n + 1];
         T sumOfMinCostOutEdges = 0;
         for (int i = 1; i <= n; i++)
         {// compute costOfMinOutEdge[i] and sumOfMinCostOutEdges
            T minCost = noEdge;
            for (int j = 1; j <= n; j++)
               if (a[i][j] != noEdge && (minCost == noEdge ||
                   minCost > a[i][j]))
                  minCost = a[i][j];
   
            if (minCost == noEdge) return noEdge; // no route
            costOfMinOutEdge[i] = minCost;
            sumOfMinCostOutEdges += minCost;
         }
      
         // initial E-node is tree root
         heapNode eNode(0, 0, sumOfMinCostOutEdges, 0, new int [n]);
         for (int i = 0; i < n; i++)
            eNode.partialTour[i] = i + 1;
         T costOfBestTourSoFar = noEdge;            // no tour found so far
         int *partialTour = eNode.partialTour;      // shorthand for
                                                    // eNode.partialTour
         // search permutation tree
         while (eNode.sizeOfPartialTour < n - 1)
         {// not at leaf
            partialTour = eNode.partialTour;
            if (eNode.sizeOfPartialTour == n - 2)
            {// parent of leaf
               // complete tour by adding two edges
               // see whether new tour is better
               if (a[partialTour[n - 2]][partialTour[n - 1]] != noEdge
                   && a[partialTour[n - 1]][1] != noEdge
                   && (costOfBestTourSoFar == noEdge ||
                      eNode.costOfPartialTour
                      + a[partialTour[n - 2]][partialTour[n - 1]]
                      + a[partialTour[n - 1]][1]
                      < costOfBestTourSoFar))
               {// better tour found
                  costOfBestTourSoFar = eNode.costOfPartialTour
                          + a[partialTour[n - 2]][partialTour[n - 1]]
                          + a[partialTour[n - 1]][1];
                   eNode.costOfPartialTour = costOfBestTourSoFar;
                   eNode.lowerCost = costOfBestTourSoFar;
                   eNode.sizeOfPartialTour++;
                   liveNodeMinHeap.push(eNode);
               }
            } 
            else
            {// generate children
               for (int i = eNode.sizeOfPartialTour + 1; i < n; i++)
                  if (a[partialTour[eNode.sizeOfPartialTour]]
                      [partialTour[i]] != noEdge)
                  {
                     // feasible child, bound path cost
                     T costOfPartialTour = eNode.costOfPartialTour
                         + a[partialTour[eNode.sizeOfPartialTour]]
                            [partialTour[i]];
                     T minAdditionalCost = eNode.minAdditionalCost
                               -  costOfMinOutEdge[partialTour
                                  [eNode.sizeOfPartialTour]];
                     T leastCostPossible = costOfPartialTour
                                           + minAdditionalCost;
                     if (costOfBestTourSoFar == noEdge ||
                       leastCostPossible < costOfBestTourSoFar)
                     {// subtree may have better leaf, put root in min heap
                         heapNode hNode(leastCostPossible,
                                        costOfPartialTour,
                                        minAdditionalCost,
                                        eNode.sizeOfPartialTour + 1,
                                        new int [n]);
                         for (int j = 0; j < n; j++)
                            hNode.partialTour[j] = partialTour[j];
                         hNode.partialTour[eNode.sizeOfPartialTour + 1] =
                                  partialTour[i];
                         hNode.partialTour[i] =
                                  partialTour[eNode.sizeOfPartialTour + 1];
                         liveNodeMinHeap.push(hNode);
                     }
                  }
               }
      
            // get next E-node
            delete [] eNode.partialTour;
            if (liveNodeMinHeap.empty()) break;
            eNode = liveNodeMinHeap.top();
            liveNodeMinHeap.pop();
         }
      
         if (costOfBestTourSoFar == noEdge)
            return NULL; // no route
   
         // copy best route into bestTour[1:n]
         for (int i = 0; i < n; i++)
            bestTour[i + 1] = partialTour[i];
      
         return costOfBestTourSoFar;
      }
      
   protected:
      void rTSP(int currentLevel)
      {// Recursive backtracking code for traveling salesperson.
       // Search the permutation tree for best tour. Start  at a node
       // at currentLevel.
         if (currentLevel == n)
         {// at parent of a leaf
            // complete tour by adding last two edges
            if (a[partialTour[n - 1]][partialTour[n]] != noEdge &&
                a[partialTour[n]][1] != noEdge &&
                (costOfBestTourSoFar == noEdge ||
                 costOfPartialTour + a[partialTour[n - 1]][partialTour[n]]
                 + a[partialTour[n]][1] < costOfBestTourSoFar))
            {// better tour found
               copy(partialTour + 1, partialTour + n + 1, bestTourSoFar + 1);
               costOfBestTourSoFar = costOfPartialTour
                                + a[partialTour[n - 1]][partialTour[n]]
                                + a[partialTour[n]][1];
            }
         }
         else
         {// try out subtrees
            for (int j = currentLevel; j <= n; j++)
               // is move to subtree labeled partialTour[j] possible?
               if (a[partialTour[currentLevel - 1]][partialTour[j]] != noEdge
                   && (costOfBestTourSoFar == noEdge ||
                    costOfPartialTour +
                      a[partialTour[currentLevel - 1]][partialTour[j]]
                      < costOfBestTourSoFar))
               {// search this subtree
                  swap(partialTour[currentLevel], partialTour[j]);
                  costOfPartialTour += a[partialTour[currentLevel - 1]]
                                        [partialTour[currentLevel]];
                  rTSP(currentLevel + 1);
                  costOfPartialTour -= a[partialTour[currentLevel - 1]]
                                        [partialTour[currentLevel]];
                  swap(partialTour[currentLevel], partialTour[j]);
               }
         }
      }

      // class data members used by btSalesperson
      static int *partialTour;
      static int *bestTourSoFar;
      static T costOfBestTourSoFar;
      static T costOfPartialTour;
   
};

int* adjacencyWDigraph<int>::partialTour;
int* adjacencyWDigraph<int>::bestTourSoFar;
int adjacencyWDigraph<int>::costOfBestTourSoFar;
int adjacencyWDigraph<int>::costOfPartialTour;
int* adjacencyWDigraph<float>::partialTour;
int* adjacencyWDigraph<float>::bestTourSoFar;
float adjacencyWDigraph<float>::costOfBestTourSoFar;
float adjacencyWDigraph<float>::costOfPartialTour;
int* adjacencyWDigraph<double>::partialTour;
int* adjacencyWDigraph<double>::bestTourSoFar;
double adjacencyWDigraph<double>::costOfBestTourSoFar;
double adjacencyWDigraph<double>::costOfPartialTour;
      
// overload <<
template <class T>
ostream& operator<<(ostream& out, const adjacencyWDigraph<T>& x)
   {x.output(out); return out;}
#endif