// cop_ref.cxx

// 

// Copyright 2002 Alexei E. Zverovitch

//

// Permission to use, copy, modify, distribute and sell this software

// for any purpose is hereby granted without fee, provided that the

// above copyright notice appears in all copies and that both that

// copyright notice and this permission notice appear in supporting

// documentation.  The author makes no representations about the

// suitability of this software for any purpose. It is provided

// "as is" without express or implied warranty.

//

// The most up-to-date version of this file can be downloaded from

// http://www.cs.rhul.ac.uk/~zvero/thesis/cop_ref/

//

// This code requires the C++ code by Jonker and Volgenant for solving

// the Linear Assignment Problem. The code can be downloaded from:

// http://www.magiclogic.com/assignment/lap_cpp.zip

//

// This code has been tested with the version 3.0 of the g++

// compiler

//

// $Id: cop_ref.cxx,v 1.3 2002/04/23 08:14:00 alexei Exp $


#include <string>
#include <vector>
#include <fstream>
#include <stdexcept>
#include <algorithm>

#include "lap.h"

using namespace std;

// Version information

string revision("$Revision: 1.3 $");
string version_temp(revision.substr(revision.find(':') + 1));
string version(version_temp.substr(0, version_temp.length() - 1));

// Cost matrix

typedef vector<vector<cost> > cost_matrix_t;
// Cycle factor - cf[i] is the vertex that follows i in the cycle factor

typedef vector<int> cycle_factor_t;
// Tour - t[i] is the vertex that follows i in the tour

typedef vector<int> tour_t;

// Infinity

const int INF = 1000000;

// The threshold parameter

const int THRESHOLD = 3;

// Check that the cycle factor is valid

void validate_cycle_factor(const cycle_factor_t& cf)
{
   vector<bool> seen(cf.size(), false);
   for (int i = 0; i != cf.size(); ++i)
   {
      if (cf[i] < 0 || cf[i] >= cf.size() || seen[i])
      {
	 throw runtime_error("Invalid cycle factor");
      };
      seen[i] = true;
   };
};

// Check that the tour is valid

void validate_tour(const tour_t& tour)
{
   validate_cycle_factor(tour);
   int n(0);
   for (int i = 0, first = 1; i != 0 || first; i = tour[i], first = 0)
   {
      ++n;
   };
   if (n != tour.size()) throw runtime_error("Invalid tour");
};

int get_tour_cost(int size, const cost_matrix_t& C, const tour_t& tour)
{
   int cost(0);
   for (int i = 0; i != size; ++i)
   {
      cost += C[i][tour[i]];
   };
   return cost;
};

// Eliminate short cycles

cycle_factor_t eliminate_short_cycles(int size, const cost_matrix_t& C)
{
   // Solve the LAP

   vector<cost*> assigncost(size);
   vector<cost> u(size), v(size);
   vector<row> colsol(size);
   vector<col> rowsol(size);
   for (int i = 0; i != size; ++i)
   {
      assigncost[i] = const_cast<cost*>(&(C[i][0]));
   };
   cost lapcost(lap(size, &assigncost[0], &rowsol[0], &colsol[0], &u[0], &v[0]));
   checklap(size, &assigncost[0], &rowsol[0], &colsol[0], &u[0], &v[0]);

   cycle_factor_t next(rowsol);
   vector<bool> visited(size, false);
   int ncycles(0); // number of cycles in the factor

   vector<pair<int, int> > contracted_vertices;
   bool contracted(false);
   for (int i = 0; i != size; ++i)
   {
      if (!visited[i]) // not-yet-visited cycle

      {
	 int len(0); // cycle length

	 int most_expensive(i);

	 // find cycle length and most expensive arc

	 for (int j = i, first = 1; j != i || first; j = next[j], first = 0)
	 {
	    ++len;
	    visited[j] = true;
	    if (C[j][next[j]] > C[most_expensive][next[most_expensive]])
	    {
	       most_expensive = j;
	    };
	 };
	 
	 if (len < THRESHOLD) // short cycle

	 {
	    // delete most expensive arc and contract the path

	    contracted_vertices.push_back(make_pair(next[most_expensive],
						    most_expensive));
	    contracted = true;
	 }
	 else // long cycle - leave vertices intact

	 {
	    for (int j = i, first = 1; j != i || first; j = next[j], first = 0)
	    {
	       contracted_vertices.push_back(make_pair(j, j));
	    };
	 };

	 ++ncycles;
      };
   };

   if (ncycles == 1 || !contracted) return next;
   
   int contracted_size(contracted_vertices.size());
   cost_matrix_t contracted_C(contracted_size);
   
   for (int i = 0; i != contracted_size; ++i)
   {
      contracted_C[i].resize(contracted_size);
      
      for (int j = 0; j != contracted_size; ++j)
      {
	 contracted_C[i][j] = C[contracted_vertices[i].second][contracted_vertices[j].first];
      };
      
      contracted_C[i][i] = INF;
   };
   
   cycle_factor_t contracted_next(eliminate_short_cycles(contracted_size, contracted_C));
   cycle_factor_t new_next(size, -1);

   for (int i = 0; i != contracted_size; ++i)
   {
      new_next[contracted_vertices[i].second] = contracted_vertices[contracted_next[i]].first;
   };
   for (int i = 0; i != size; ++i)
   {
      if (new_next[i] == -1) new_next[i] = next[i]; // expand contracted paths

   };
   validate_cycle_factor(new_next);

   return new_next;
};

bool compare_cycle_info(const pair<int,int> c1, const pair<int,int> c2)
{
   return c1.first > c2.first;
};

// Transforms cycle factor next into a tour (in-place)

void KSP(int size, const cost_matrix_t& C, cycle_factor_t& next)
{
   vector<bool> visited(size, false);
   typedef vector<pair<int,int> > cycle_info_t;
   cycle_info_t cycle_info;
   
   for (int i = 0; i != size; ++i)
   {
      if (!visited[i]) // not-yet-visited cycle

      {
	 int len(0); // cycle length


	 // find cycle length

	 for (int j = i, first = 1; j != i || first; j = next[j], first = 0)
	 {
	    ++len;
	    visited[j] = true;
	 };

	 cycle_info.push_back(make_pair(len, i));
      }
   };

   sort(cycle_info.begin(), cycle_info.end(), compare_cycle_info);

   pair<int,int> c0(*cycle_info.begin());
   for (cycle_info_t::iterator ci = cycle_info.begin() + 1; ci != cycle_info.end(); ++ci)
   {
      int min_patching_cost(INF);
      int min_i;
      int min_j;

      for (int i = c0.second, first = 1; i != c0.second || first; i = next[i], first = 0)
      {
	 for (int j = ci->second, first = 1;
	      j != ci->second || first;
	      j = next[j], first = 0)
	 {
	    int patching_cost(C[i][next[j]] + C[j][next[i]] - C[i][next[i]] - C[j][next[j]]);

	    if (patching_cost < min_patching_cost)
	    {
	       min_i = i;
	       min_j = j;
	       min_patching_cost = patching_cost;
	    };
	 };
      };

      swap(next[min_i], next[min_j]);
   };

   validate_tour(next);
};

// Implementation of KSP/COP

tour_t COP(const string& filename)
{
   int size;
   cost_matrix_t C;

   // Read the ATSP instance

   cout << "Solving " << filename << endl;
   ifstream f(filename.c_str());
   if (!f) throw runtime_error("File not found: " + filename);
   f >> size;
   C.resize(size);
   for (int i = 0; i != size; ++i)
   {
      C[i].resize(size);
      for (int j = 0; j != size; ++j)
      {
	 f >> C[i][j];
      };
      C[i][i] = INF;
   };

   // Apply KSP/COP

   cycle_factor_t next(eliminate_short_cycles(size, C));
   KSP(size, C, next);
   cout << "Tour cost: " << get_tour_cost(size, C, next) << endl;
   return next;
};

// Format of input file:

//

// num_vertices

// c(1,1) c(1,2) ... c(1,n)

// . . .

// c(n,1) c(n,2) ... c(n,n)

int main(int argc, char** argv)
{
   try
   {
      cout << "cop_ref version" << version << "built on " __DATE__ << endl;
      if (argc < 2)
      {
	 cout << "Usage: cop_ref file ..." << endl;
      }
      else
      {
	 for (int i = 1; i < argc; ++i) COP(argv[i]);
      };
      return EXIT_SUCCESS;
   }
   catch (exception& e)
   {
      cerr << "Error: " << e.what() << endl;
      return EXIT_FAILURE;
   };
};