TSP branch bound with Held-Karp lower bound

branch_bound_HK.cc

//
// Held and Karp bound
//
// T が 1-tree iff T は {2,...,n} 上の全域木 + 1 から枝が 2 本．
// 1-tree かつ全部の頂点の誘導次数が 2 であればそれはサイクル．
//
// 巡回セールス人を
//   minimize     c(T) 
//   subject to   T is a 1-tree
//                deg_T(i) = 2 for all i
// と定式化する．次数制約が辛いので目的関数に入れて
//   min_T max_p c(T) + \sum p_i (deg_T(i) - 2)
// とする（ここまで同値変形）．
//
// min と max の順序を入れ替えて双対問題を作る（双対ギャップあり）
//   max_p min_T c(T) + \sum p_i (deg_T(i) - 2)
// これを解くと min max >= max min より，元の問題の下界が得られる．
// この下界を Held-Karp 下界という．
//
// 双対問題を劣勾配法で解く．
// 内側の min_T 以下は重みを c'(i,j) = c(i,j) + p_i + p_j とした
// グラフで最小 1-tree を計算すればいい．最小 1-tree は最小全域木と
// 同じようにして計算可能．あとは p_i <- p_i + eps (deg_T(i) - 2)
// の更新式で p に関する最大化ができる．
//
// これを branch-bound するときの下界に使って tsp を解いてみる．
//

#include <iostream>
#include <vector>
#include <cstdio>
#include <algorithm>
#include <queue>
#include <functional>

using namespace std;

#define fst first
#define snd second
#define all(c) ((c).begin()), ((c).end())
#define dout if(0)cout


#include <string>
string bitstr(unsigned long long x, int n = 6) {
  char s[n+1];
  for (int i = 0; i < n; ++i) 
    s[n-i-1] = ((x >> i) & 1) + '0';
  s[n] = 0;
  return s;
}


struct union_find {
  vector<int> p; 
  union_find(int n) : p(n, -1) { };
  bool unite(int u, int v) { 
    if ((u = root(u)) == (v = root(v))) return false;
    if (p[u] > p[v]) swap(u, v);
    p[u] += p[v]; p[v] = u;
    return true;
  }
  bool find(int u, int v) { return root(u) == root(v); }
  int root(int u) { return p[u] < 0 ? u : p[u] = root(p[u]); }
  int size(int u) { return -p[root(u)]; }
};

template <class T>
struct graph {
  struct edge {
    int src, dst;
    T weight;
  };
  int n, m;
  vector<edge> edges;
  vector<vector<int>> adj;
  graph(int n) : n(n), adj(n), m(0) { }
  void add_edge(int src, int dst, T weight) {
    adj[src].push_back(edges.size());
    adj[dst].push_back(edges.size());
    edges.push_back({src, dst, weight});
    ++m;
  }

  struct state {
    vector<double> p; // dual ver
    vector<int> b;    // fixed edges; 1 for forbiddedn
    vector<int> x;    // current solution
    vector<int> d;    // degree
    double score;
    bool operator<(const state &s) const {
      return score > s.score;
    };
  };
  bool evaluate(state &s) {
    double eta = 1.0;
    vector<int> id(m); iota(all(id), 0);
    for (int iter = 0; iter < 20; ++iter) {
      //dout << "iter " << iter << endl;
      auto residue = [&](int i) {
        if (s.b[i] == 1) return 1.0/0.0;
        return edges[i].weight + s.p[edges[i].src] + s.p[edges[i].dst];
      };
      sort(all(id), [&](int i, int j) {
        if (s.b[i] != s.b[j]) return s.b[i] < s.b[j]; // -1 < 0 < 1
        return residue(i) < residue(j);
      });
      union_find uf(n);
      s.d.assign(n, 0);
      s.score = 0;
      s.x.assign(m, 0);
      for (int i: id) {
        const edge &e = edges[i];
        if (((e.src == 0 || e.dst == 0) && s.d[0] < 2) || 
            ((e.src != 0 && e.dst != 0) && uf.unite(e.src, e.dst))) {
          ++s.d[e.src];
          ++s.d[e.dst];
          s.score += residue(i);
          s.x[i] = 1;
        }
      }
      eta = 1.0 / (1.0 + iter);
      bool found = true;
      for (int i = 1; i < n; ++i) {
        if (s.d[i] != 2) found = false;
        s.p[i] += eta * (s.d[i] - 2);
      }
      if (found) return true;
    }
    return false;
  }

  double held_karp() {
    if (adj[0].size() < 2) return 1.0/0.0;
    priority_queue<state> que;
    state s = {vector<double>(n), vector<int>(m)};

    int expanded_nodes = 0;

    double UB = 1.0 / 0.0;
    if (evaluate(s)) {
      //printf("%.8lf\n", s.score);
      return round(s.score);
    }
    que.push(s);
    while (!que.empty()) {
      s = que.top(); que.pop();
      ++expanded_nodes;
      if (s.score >= UB) continue;

      dout << "LB = " << s.score << " / UB = " << UB << endl;
      dout << "expand new state" << endl;
      for (int e = 0; e < m; ++e) {
        if (s.b[e] == 1)  
          dout << "  " << edges[e].src << " " << edges[e].dst << " is excluded" << endl;
        if (s.b[e] ==-1)  
          dout << "  " << edges[e].src << " " << edges[e].dst << " is included" << endl;
        if (s.x[e]) 
          dout << "  " << edges[e].src << " " << edges[e].dst << " is used" << endl;
      }
      dout << endl;

      int u = 1;
      for (; u < n; ++u) 
        if (s.d[u] > 2) break;
      dout << "  bounding vertex " << u << endl;

      int freedom = 2;
      for (int e: adj[u]) 
        if (s.b[e] ==-1) --freedom;
      vector<int> x = s.x;
      for (int e: adj[u]) {
        if (x[e] == 1 && s.b[e] == 0) {
          if (freedom-- > 0) {
            dout << "... try to exclude " << edges[e].src << " " << edges[e].dst << endl;
            s.b[e] = 1;
            if (evaluate(s)) {
              dout << "  ==> found a feasible solution with score " << s.score << endl;
              UB = min(UB, s.score);
            } else if (s.score < UB) {
              dout << "  ==> insert queue" << endl;
              que.push(s);
            }
            s.b[e] = -1;
          } else {
            s.b[e] =  1;
          }
        }
      }
      dout << "... use first twos " << endl;
      if (evaluate(s)) {
        dout << "  ==> found a feasible solution with score " << s.score << endl;
        UB = min(UB, s.score);
      } else if (s.score < UB) {
        dout << "  ==> insert queue" << endl;
        que.push(s);
      }
    }
    cout << "expanded nodes: " << expanded_nodes << endl;
    //printf("%.8lf\n", UB);
    return round(UB);
  }


  // Held and Karp's dynamic programming
  int dynamic_programming() {

    /*
      vector<int> ord(n);
      iota(all(ord), 0);
      int opt = 99999999;
      do {
        int ans = 0;
        for (int e = 0; e < m; ++e) {
          int u = edges[e].src;
          int v = edges[e].dst;
          int w = edges[e].weight;
          if ((ord[u]+1)%n == ord[v] || (ord[v]+1)%n == ord[u]) ans += w;
        }
        opt = min(opt, ans);

        if (ans == 191) {
          int ans = 0;
          for (int e = 0; e < m; ++e) {
            int u = edges[e].src;
            int v = edges[e].dst;
            int w = edges[e].weight;
            if ((ord[u]+1)%n == ord[v] || (ord[v]+1)%n == ord[u]) {
              ans += w;
              dout << "(" << u << ", " << v << "; " << w << ") ";
            }
          }
          dout << endl;
        }
      } while (next_permutation(ord.begin()+1, ord.end()));
      return opt;
    }
    */


    vector<vector<int>> f(n, vector<int>(1<<n, 99999999));
    f[0][1] = 0;
    for (int S = 0; S < (1 << n); ++S) {
      for (int e = 0; e < m; ++e) {
        int u = edges[e].src;
        int v = edges[e].dst;
        int w = edges[e].weight;
        if (!(S & (1 << u))) f[u][S|(1<<u)] = min(f[u][S|(1<<u)], w + f[v][S]);
        if (!(S & (1 << v))) f[v][S|(1<<v)] = min(f[v][S|(1<<v)], w + f[u][S]);
      }
    }
    int opt = 99999999;
    for (int e = 0; e < m; ++e) {
      int u = edges[e].src;
      int v = edges[e].dst;
      int w = edges[e].weight;
      if (u == 0) opt = min(opt, f[v][(1<<n)-1] + w);
      if (v == 0) opt = min(opt, f[u][(1<<n)-1] + w);
    }
    return opt;
  }
};


// === tick a time ===
#include <ctime>
double tick() {
  static clock_t oldtick;
  clock_t newtick = clock();
  double diff = 1.0*(newtick - oldtick) / CLOCKS_PER_SEC;
  oldtick = newtick;
  return diff;
}

int main() {

  double ta = 0, tb = 0;
//  for (int seed = 38772; ; ++seed) {
  for (int seed = 1; seed <= 7; ++seed) {
    dout << "seed = " << seed << endl;
    srand(seed);

    int n = 40;
    graph<int> g(n);
    dout << "input: " << endl;
    for (int i = 0; i < n; ++i) {
      for (int j = i+1; j < n; ++j) {
        g.add_edge(i, j, 1 + rand() % 1000) ;
        dout << i << " " << j << " " << g.edges.back().weight << endl;
      }
    }
    /*
    for (int i = 0; i < n; ++i) {
      for (int j = i+1; j < n; ++j) {
        g.add_edge(i, j, 1 + rand() % 1000) ;
        dout << i << " " << j << " " << g.edges.back().weight << endl;
      }
    }
    */
    tick();
    int a = 0; //g.dynamic_programming();
    ta += tick();
    int b = g.held_karp();
    tb += tick();
    cout <<  a << " " <<  b << endl;
    cout << ta << " " << tb << endl;
    if (a != b) {
      if (a >= 99999999) continue;
      cout << "seed = " << seed << endl;
      cout << "DP = " << a << " / branch-bound = " << b << endl;
      break;
    }
  }

}