TravellingSalesman.cpp

// IE511: Integer Programming
// Hint for the TSP Solver
// C++ Program written by Prof. R.S. Sreenivas
//
// The input file is to have the following format:
// Line 1: #vertices
// Line 2: First row of the cost matrix ie. a (1 x #vertices) row
// Line 3: Second row of the cost matrix ie. a (1 x #vertices) row
// Line 4: Third row of the cost matrix ie. a (1 x #vertices) row
// etc. etc. etc.

#include <iostream>
#include <iomanip>
#include <cmath>
#include <fstream>
#include <cstdlib>
#include <vector>

#include "lp_lib.h"

using namespace std;

// Setting the appropriate pointers for lpsolve
// keep in mind the size of the solution vector will depend on the data
// it has to the allocated dynamically.
lprec *lp;
double *solution;

// some global variables
int number_of_vertices;

// This routine is used to generate the next mask, which in turn
// is used to generate all possible subsets of the set of
// vertices using the procedure from the following source
// http://compprog.wordpress.com/2007/10/10/generating-subsets/
//
int next(int mask[], int n)
{
	int i;
	for (i = 0; (i < n) && mask[i]; ++i)
		mask[i] = 0;
	
	if (i < n)
	{
		mask[i] = 1;
		return 1;
	}
	return 0;
}

// I need this to check how many vertices are in a subset that
// is generated by the mask-procedure from
// http://compprog.wordpress.com/2007/10/10/generating-subsets/
// I need this for the RHS of the subtour elimination constraints
int size_of_mask(int mask[], int n)
{
	int size = 0;
	for (int i = 0; i < n; i++)
		if (mask[i] == 1)
			size++;
	return (size);
}

// This routine sets the ILP up based on the input data
void initialize_lp_from_input_data(char* argv[])
{
	// reading the input filename from commandline
	ifstream input_filename(argv[1]);
	
	if (input_filename.is_open()) {
		
		cout << "Input File Name: " << argv[1] << endl;
		
		// first line contains the number of items
		input_filename >> number_of_vertices;
		cout << "Number of Vertices = " << number_of_vertices << endl;
		
		// initialize the lp
		// there is a variable for each pair of vertices
		lp = make_lp(0, number_of_vertices*number_of_vertices);
		
		// now that we know the number of items, we can
		// set the size of the solution arrays
		solution = new double[number_of_vertices*number_of_vertices];
		
		// the row that defines the objective function should have
		// 1 + number_of_vertices*number_of_vertices items in it...
		// the 0-th element of this row is ignored.
		double objective_function[1+number_of_vertices*number_of_vertices];
		objective_function[0] = 0;
		// initializing the rest of the objective function
		{
			// fill the necessary code here!
            for (int i=1; i<=number_of_vertices*number_of_vertices; i++)
            {
                int value_just_read;
				input_filename >> value_just_read;
                objective_function[i]=value_just_read;
                
            }
        }
		// set the objective function
		set_obj_fn(lp, objective_function);
		
		// This keeps the message reporting of lp_solve to a minimum
		set_verbose(lp, 3);
		
		// the interpretation for the decision variables is this x(i,j) is
		// 1 (0) if the optimal tour takes (does not take) the arc from i to j
		{
			// Fill the necessary code to make sure the tour leaves every vertex
            
            for (int i=1; i<=number_of_vertices; i++)
            {
                double row[number_of_vertices*number_of_vertices];
                for (int k=0;k<number_of_vertices*number_of_vertices+1;k++)
                    row[k]=0;
                for (int j=1; j<=number_of_vertices; j++)
                {
                    if (i!=j)
                        row[j+number_of_vertices*(i-1)]=1;
                    
                }
                add_constraint(lp, row, EQ, 1);
            }
            
		}
		
		{
			// Fill the necessary code to make sure the tour enters every vertex
            for (int i=1; i<=number_of_vertices; i++)
            {
                double row[number_of_vertices*number_of_vertices];
                for (int k=0;k<number_of_vertices*number_of_vertices+1;k++)
                    row[k]=0;
                for (int j=1; j<=number_of_vertices; j++)
                {
                    if (i!=j)
                        row[i+number_of_vertices*(j-1)]=1;
                    
                }
                add_constraint(lp, row, EQ, 1);
            }
            
            
		}
		
		// we do not have to worry about forcing all x(i,i) variables to zero
		// (why?) -- because this is automatically done when we generate the
		// subtour elimination constraints for all singleton sets
		
		// subtour elimination constraints
		{
			int mask[number_of_vertices];
			
			// initialize mask to all zeros
			for (int i = 0; i < number_of_vertices; i++)
				mask[i] = 0;
			
			// start adding constraints for each
			while ((next(mask, number_of_vertices)==1) &&
				   (size_of_mask(mask, number_of_vertices) < number_of_vertices))
			{
				// this part of the code will generate all sub-tours
                // you have to figure out a way of putting it into constraints
                // that are lp_solve friendly
                double row[number_of_vertices*number_of_vertices+1];
                for (int k=0;k<number_of_vertices*number_of_vertices+1;k++)
                    row[k]=0;
                
                for(int i=1;i<=number_of_vertices;i++)
                    
                {
                    for(int j=1;j<=number_of_vertices;j++)
                        if(i!=j && mask[i]==1 && mask[j]==1)
                        {
                            row[(i-1)*number_of_vertices+j]=1;
                        }
                    
                }
                
              add_constraint(lp, row, LE, size_of_mask(mask, number_of_vertices) - 1);
			}
		}
        
		// put the required code to make sure all variables are binary
        for(int k=1;k<number_of_vertices*number_of_vertices+1;k++)
            set_binary(lp,k,TRUE);
            }
	else {
		cout << "Input file missing" << endl;
		exit(0);
	}
}

// This routine solves the TSP instance.
// The procedure of extracting the optimal TSP path from the solution from
// lpsolve takes some effort to follow -- otherwise this is quite routine.
void solve_the_TSP()
{
	int ret;
	
	// solve the ILP
	ret = solve(lp);
	
	if (ret == 0)
	{
		// get the optimal assignment
		get_variables(lp, solution);
		
		cout << "Optimal Tour" << endl;
		int current_vertex = 1;
		cout <<  current_vertex;
		for (int j = 1; j <= number_of_vertices; j++)
		{
			if (solution[((current_vertex-1)*number_of_vertices)+j-1] == 1)
			{
				current_vertex = j;
				cout << " --> " << j;
			}
		}
		
		// we would have left vertex 1, and when we get back to vertex 1, we
		// are done with the tour.
		while (current_vertex != 1)
		{
			for (int j = 1; j <= number_of_vertices; j++)
			{
				if ((solution[((current_vertex-1)*number_of_vertices)+j-1] == 1) &&
					(current_vertex != 1))
				{
					current_vertex = j;
					cout << " --> " << j;
				}
			}
		}
		cout << endl;
		
		// optimal cost of the Hamiltonian Tour
		cout << "Optimal Tour Length = " << get_objective(lp) << endl;
	}
	else {
		cout << "Problem is infeasible, check input file... something went wrong" << endl;
	}
}

int main (int argc, char* argv[])
{
	// formulate the the TSP ILP from input data
	initialize_lp_from_input_data(argv);
    
	// solve the TSP instance
	solve_the_TSP();
	
	return(0);
}