Traveling Salesman

Brandon Blatter

CS 312

View this report in your browser

Let's be honest. PDFs were not meant for displaying code, and code was not meant to be displayed in PDFs. Please view this report in your browser so it's easy to read.

Code

State.cs

using System;
using System.Collections;
using System.Collections.Generic;
using System.Drawing;

namespace TSP {

    // custom class that manages data associated with each node
    public class State {

        public State(int cityIndex, ArrayList path, double lb, double[,] matrix) {
            cityListIndex = cityIndex;
            partialPath = path;
            lowerBound = lb;
            reducedCostMatrix = matrix;
            weight = lb / (path.Count + 1);
        }

        public int cityListIndex;
        public int listIndex;
        public ArrayList partialPath;
        public double lowerBound;
        public double[,] reducedCostMatrix;
        public double weight;

    }
}

BinaryHeapPQ.cs

using System;
using System.Collections.Generic;
using System.Drawing;

namespace TSP {

    public class BinaryHeapPQ {

        public BinaryHeapPQ() {
            // reset list of states
            heapList = new List<State>();
            heapLocations = new List<int>();
        }

        public BinaryHeapPQ(List<State> states) {
            MakeQueue(states);
        }

        // locations of nodes on the heap (access from DeleteMin)
        private List<int> heapLocations;
        private List<State> heapList;

        public bool IsEmpty() {
            return heapList.Count == 0;
        }

        // gets node's parent from binary heap by index
        private int GetParentIndex(int nodeIndex) {
            return nodeIndex / 2;
        }

        private int GetSmallerChildIndex(int nodeIndex) {
            int childIndex1 = nodeIndex * 2;
            // if left child index is greater than size of heap, use nodeIndex
            if(childIndex1 >= heapList.Count) return nodeIndex;
            int childIndex2 = childIndex1 + 1;
            // if right child index is greater than size of heap, use left child index
            if(childIndex2 >= heapList.Count) return childIndex1;
            State child1 = heapList[childIndex1];
            State child2 = heapList[childIndex2];
            // if right child doesn't exist or has greater weight, use left child
            if(child1.weight < child2.weight) return childIndex1;
            // otherwise use right child
            else return childIndex2;
        }

        // swap nodes by index in heap list
        private void SwapNodes(int nodeIndex1, int nodeIndex2) {
            // get nodes to switch
            State node1 = heapList[nodeIndex1];
            State node2 = heapList[nodeIndex2];
            // swap nodes actual nodes in heapList
            heapList[nodeIndex2] = node1;
            heapList[nodeIndex1] = node2;
            // swap indices in heapLocations
            heapLocations[node1.listIndex] = nodeIndex2;
            heapLocations[node2.listIndex] = nodeIndex1;
        }

        private void BubbleUp(int nodeIndex) {
            int parentIndex = GetParentIndex(nodeIndex);
            // while parent value is higher
            while(heapList[parentIndex].weight > heapList[nodeIndex].weight) {
                // swap nodes
                SwapNodes(parentIndex, nodeIndex);
                // update current node and parent indices
                nodeIndex = parentIndex;
                parentIndex = GetParentIndex(parentIndex);
            }
        }

        private void BubbleDown(int nodeIndex) {
            int smallerChildIndex = GetSmallerChildIndex(nodeIndex);
            if(smallerChildIndex == 0) return;
            // while child value is lower
            while(heapList[smallerChildIndex].weight < heapList[nodeIndex].weight) {
                // swap nodes
                SwapNodes(smallerChildIndex, nodeIndex);
                // update current node and parent indices
                nodeIndex = smallerChildIndex;
                smallerChildIndex = GetSmallerChildIndex(smallerChildIndex);
                if(smallerChildIndex == 0) return;
            }
        }

        // MakeQueue: O(n log n)
        public void MakeQueue(List<State> states) {
            // reset list of states
            heapList = new List<State>();
            heapLocations = new List<int>();
            // add each state given; order doesn't matter
            for (int i = 0; i < states.Count; i++) {
                Insert(states[i]);
            }
        }

        // Insert: O(log n)
        public void Insert(State state) {
            // add state to end of list (bottom right of tree)
            heapList.Add(state);
            // add location as the last element in the heap list
            heapLocations.Add(heapList.Count - 1);
            state.listIndex = heapLocations.Count - 1;
            // bubble up
            BubbleUp(heapList.Count - 1);
        }

        // DeleteMin: O(log n)
        public State DeleteMin() {
            // get top node
            State min = heapList[0];
            heapLocations[min.listIndex] = -1;
            // replace top node with bottommost rightmost
            heapList[0] = heapList[heapList.Count - 1];
            heapList.RemoveAt(heapList.Count - 1);
            if(heapList.Count > 0) {
                heapLocations[heapList[0].listIndex] = 0;
            }
            // bubble down
            BubbleDown(0);
            return min;
        }

        // DecreaseKey: O(log n)
        public void DecreaseKey(int listIndex) {
            // bubble up actual location in heap
            BubbleUp(heapLocations[listIndex]);
        }
    }
}

BranchAndBound.cs

using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;

namespace TSP {

    class BranchAndBound {

        public BranchAndBound() {

        }

        // for storing the result of reducing a matrix
        public class MatrixAndBound {

            public MatrixAndBound(double[,] costMatrix, double oldLowerBound) {
                reducedCostMatrix = costMatrix;
                lowerBound = oldLowerBound;
            }

            public double[,] reducedCostMatrix;
            public double lowerBound;

        }

        public static double[,] CreateCostMatrix(City[] cities) {
            int length = cities.Length;
            double[,] costMatrix = new double[length, length];
            // in every cell
            for(int i = 0; i < length; i++) {
                for(int j = 0; j < length; j++) {
                    // store cost between cities
                    costMatrix[i, j] = cities[i].costToGetTo(cities[j]);
                }
            }
            // turn middle distance zeros to infinities
            for(int i = 0; i < length; i++) {
                costMatrix[i, i] = double.PositiveInfinity;
            }
            return costMatrix;
        }

        public static MatrixAndBound PrepareMatrix(    double[,] intialMatrix,
                                                                                                int rowIndex,
                                                                                                int columnIndex,
                                                                                                double oldLowerBound) {
            double[,] matrix = intialMatrix.Clone() as double[,];
            // add cost of this path
            double lowerBound = oldLowerBound + matrix[rowIndex, columnIndex];
            // infinity out row
            for(int i = 0; i < matrix.GetLength(0); i++) {
                matrix[rowIndex, i] = double.PositiveInfinity;
            }
            // infinity out column
            for(int i = 0; i < matrix.GetLength(1); i++) {
                matrix[i, columnIndex] = double.PositiveInfinity;
            }
            // infinity out reverse path
            matrix[columnIndex, rowIndex] = double.PositiveInfinity;
            // return prepared result
            MatrixAndBound result = new MatrixAndBound(matrix, lowerBound);
            return result;
        }

        public static MatrixAndBound ReduceMatrix(MatrixAndBound start) {
            double [,] matrix = start.reducedCostMatrix.Clone() as double[,];
            double lowerBound = start.lowerBound;
            // reduce each row
            double lowestRowValue;
            for(int i = 0; i < matrix.GetLength(0); i++) {
                // find the lowest non-infinity value
                lowestRowValue = double.PositiveInfinity;
                for(int j = 0; j < matrix.GetLength(1); j++) {
                    if(!double.IsPositiveInfinity(matrix[i, j]) && matrix[i, j] < lowestRowValue) {
                        lowestRowValue = matrix[i, j];
                    }
                }
                // if lowest non infinite value > 0, subtract this value from every non-infinite value
                if(lowestRowValue > 0 && !double.IsPositiveInfinity(lowestRowValue)) {
                    // add cost difference to lower bound
                    lowerBound += lowestRowValue;
                    for(int j = 0; j < matrix.GetLength(1); j++) {
                        //if(!double.IsPositiveInfinity(matrix[i, j]))
                        matrix[i, j] -= lowestRowValue;
                    }
                }
            }
            // reduce each column
            double lowestColumnValue;
            for(int i = 0; i < matrix.GetLength(1); i++) {
                // find the lowest non-infinity value
                lowestColumnValue = double.PositiveInfinity;
                for(int j = 0; j < matrix.GetLength(0); j++) {
                    if(!double.IsPositiveInfinity(matrix[j, i]) && matrix[j, i] < lowestColumnValue) {
                        lowestColumnValue = matrix[j, i];
                    }
                }
                // if lowest value > 0, subtract this value from every non-infinite value
                if(lowestColumnValue > 0 && !double.IsPositiveInfinity(lowestColumnValue)) {
                    // add cost difference to lower bound
                    lowerBound += lowestColumnValue;
                    for(int j = 0; j < matrix.GetLength(0); j++) {
                        if(!double.IsPositiveInfinity(matrix[j, i])) matrix[j, i] -= lowestColumnValue;
                    }
                }
            }

            return new MatrixAndBound(matrix, lowerBound);
        }

        public static void PrintMatrix(double[,] matrix) {
            Console.Write("\n");
            Console.Write("\n");

            for(int i = 0; i < matrix.GetLength(0); i++) {
                for(int j = 0; j < matrix.GetLength(1); j++) {
                    Console.Write((double.IsPositiveInfinity(matrix[i, j]) ? "∞" : matrix[i, j].ToString()) + "\t");
                }
                Console.Write("\n");
            }
        }

        public static bool IsValidPath(ArrayList path, int numCities) {
            // not valid if count differs from number of cities
            if(path.Count != numCities) return false;

            // go through each edge in the route and add up the cost
            int x;
            City here;
            double cost = 0D;

            for (x = 0; x < path.Count - 1; x++) {
                here = path[x] as City;
                cost += here.costToGetTo(path[x + 1] as City);
            }
            // go from the last city to the first
            here = path[path.Count - 1] as City;
            cost += here.costToGetTo(path[0] as City);
            // not valid if cost is infinite
            if(double.IsPositiveInfinity(cost)) return false;

            // only valid if every item in the list is unique
            return path.Cast<City>().Distinct().Count() == path.Count;
        }

    }
}

ProblemAndSolver.cs

The following is the most relevant excerpt.

/// <summary>
/// performs a Branch and Bound search of the state space of partial tours
/// stops when time limit expires and uses BSSF as solution
/// </summary>
/// <returns>results array for GUI that contains three ints: cost of solution, time spent to find solution, number of solutions found during search (not counting initial BSSF estimate)</returns>
public string[] bBSolveProblem() {

    string[] results = new string[3];
    Stopwatch timer = new Stopwatch();
    BinaryHeapPQ branchAndBoundPQ = new BinaryHeapPQ();
    City[] bbCities = GetCities();
    ArrayList newPartialPath;
    MatrixAndBound newResult;
    TSPSolution possibleNewSolution;
    int count = 0;
    double currentCostOfBssf = double.PositiveInfinity;

    timer.Start();

    // create a matrix of costs
    double[,] costMatrix = CreateCostMatrix(GetCities());

    // produce a reduced costs matrix
    MatrixAndBound originalResult = ReduceMatrix(new MatrixAndBound(costMatrix, 0));

    // pick a random city to start
    Random r = new Random();
    int firstCityIndex = r.Next(bbCities.Length);

    // add to queue
    branchAndBoundPQ.Insert(new State(    firstCityIndex,
                                                                            new ArrayList(),
                                                                            originalResult.lowerBound,
                                                                            originalResult.reducedCostMatrix));

    State currentState;
    // while there is still stuff in the queue and time hasn't run out
    while(!branchAndBoundPQ.IsEmpty()) {
        // halt if over a certain time limit
        // if(timer.ElapsedMilliseconds > 30000) break;
        // store current cost of bssf instead of calculating it all the time
        currentCostOfBssf = bssf == null ? double.PositiveInfinity : costOfBssf();
        // get the current state for operations; delete the first item in the queue
        currentState = branchAndBoundPQ.DeleteMin();

        if(currentState.lowerBound > currentCostOfBssf) continue;

        // make a new partial path to add to children
        newPartialPath = new ArrayList(currentState.partialPath);
        newPartialPath.Add(bbCities[currentState.cityListIndex]);

        // make a new TSPSolution for testing
        possibleNewSolution = new TSPSolution(newPartialPath);
        // forget about it if it costs more than the best one
        double costOfNewSolution = possibleNewSolution.costOfRoute();
        if(costOfNewSolution > currentCostOfBssf) continue;

        // update the best solution so far (bssf) if possible
        if(IsValidPath(newPartialPath, bbCities.Length)) {
            bssf = possibleNewSolution;
            count++;
            // path is valid; don't process children
            continue;
        }

        // iterate over every city
        for(int i = 0; i < bbCities.Length; i++) {
            // if there is viable path to the next city
            if(            !double.IsPositiveInfinity(currentState.reducedCostMatrix[currentState.cityListIndex, i])
                    &&    i != firstCityIndex) {
                // calculate the cost matrix and make a little object
                newResult = ReduceMatrix(    PrepareMatrix(    currentState.reducedCostMatrix,
                                                                                                    currentState.cityListIndex,
                                                                                                    i,
                                                                                                    currentState.lowerBound));
                // add to the queue
                branchAndBoundPQ.Insert(new State(    i,
                                                                                        newPartialPath,
                                                                                        newResult.lowerBound,
                                                                                        newResult.reducedCostMatrix));
            }
        }
    }

    timer.Stop();

    results[COST] = costOfBssf().ToString();
    results[TIME] = timer.Elapsed.ToString();
    results[COUNT] = count.ToString();

    return results;
}

Complexity

Time Complexity

The overall Big O for this algorithm is O(n! n²) and consists of two main parts: the time taken to add each state object to the queue and the matrix operations performed on each state object.

Overall, as many as O(n!) objects may need to be added to the queue. This represents the total number of possible paths that may be taken, each which could have a corresponding

Each object on the queue needs to go under a O(n²) operation so it can pass over a row and a column in the matrix. This is required for changing a row and a column's values to infinity as well as subtracting the smallest value from each cell in a row and a column.

My binary heap has a Big O of O(log n) but because this n is different from the number of cities, it is not considered in the final Big O calculation.

Space Complexity

Since every state that is created needs to be stored along with its matrix, the space complexity has the same big O as the time complexity: O(n! n²).

State Data Structure

I created a very simple state data structure (State.cs) to hold information about each state on the queue. The state's data members are as follows:

• int cityListIndex: the index of the "from" city as found in the array of cities
• int listIndex: the index of the state in the binary heap queue
• ArrayList partialPath: the cities in the path this state is a part of
• double lowerBound: the lower bound of the final cost found by following this path
• double[,] reducedCostMatrix: the reduced cost matrix for this state
• double weight: the weight in the priority queue (smaller has higher priority), which is calculated using the lower bound and the number of cities in the partial path

I create and add a state to the queue each time a path has a possible solution with a better lower bound than the best solution so far. States are removed from the priority queue in order of weight and processed individually.

Priority Queue Data Structure

The priority queue data structure is a simple binary heap. I repurposed the code I wrote for the network routing lab here. All operations are O(log n).

Whenever a node is added to or removed from the binary heap, the heap rebalances itself so the nodes stay sorted. Nodes are added to the end and then they "bubble up" to the top (swapping with parents if the new node has a lower weight). Nodes that are removed are replaced with the last node (the node with the highest weight) which then "bubbles down" to its correct position.

Approach for Initial BSSF

Before I have any solution to a TSP problem I continue to explore and process states in order of their weight in the queue. Any potential states are added to the queue, regardless of their upper bound. Once an initial solution is found, I check its cost against the lower bound of states being pulled off the queue. If the states coming from the queue end up having a lower bound higher than the initial solution, I don't process them.

After that, I simply update the initial solution with a new solution every time a full path is found with a smaller cost.

Outcomes

All times are measured in seconds.

Raw

Binary Heap

Discussion

Generally problems tend to get much much harder to solve with every city added. While good pruning algorithms may be able to keep the maximum number of stored states low, the number of total states calculated is almost unavoidably high.

As explained earlier, the time complexity and space complexity grow very quickly. Pruning algorithms are effective at reducing the space complexity, but the time complexity is unavoidable as every state in the queue must be calculated at some point.