cup

aco.cc

@@ -13,19 +13,13 @@
 #include <iostream>
 #include <algorithm>
 #include <unistd.h>
+#include "opt.h"
+
+#define SINGLE_CORE 1
 
 using namespace std;
 typedef unsigned int uint;
 
-ostream& operator<< (ostream& out, vector<uint> a) {
-  out << "[";
-  for (int i = 0; i < a.size() - 1; i++) {
-    out << a[i] << ",";
-  }
-  out << a[a.size() - 1] << "]";
-  return out;
-}
-
 const uint MAX_NODES = 1577;
 
 // alpha >= 0
@@ -49,16 +43,11 @@ vector<uint> sh_route;
 bool first_pass;
 uint start;
 
-uint select_random(const set<uint> &s) {
-  auto n = rand() % s.size(); // not _really_ random
-  auto it = begin(s);
-  advance(it, n); // 'advance' the iterator n times
-  return *it;
-}
 
 struct ant {
-  uint other;
+  uint seed;
   uint idx;
+  uint iter;
   pthread_t t_handle;
   uint location;
   double distance_travelled;
@@ -66,8 +55,15 @@ struct ant {
   vector<uint> route;
   bool tour_complete = false;
 
+  uint select_random(const set<uint> &s) {
+    auto n = rand_r(&seed) % s.size(); // not _really_ random
+    auto it = begin(s);
+    advance(it, n); // 'advance' the iterator n times
+    return *it;
+  }
+
   void init() {
-    this->other = 0;
+    this->seed = this->idx * 1000 + this->iter;
 
     this->distance_travelled = 0.0;
 
@@ -78,7 +74,6 @@ struct ant {
       this->possible_locations.insert(i);
     }
 
-
     this->location = start;
 
     this->update_route(start);
@@ -95,6 +90,8 @@ struct ant {
     if (first_pass) {
       return select_random(this->possible_locations);
     } else {
+      double toss = (double) rand_r(&seed) / (double) RAND_MAX;
+
       double attractiveness[n_nodes];
       double sum_total = 0.0;
       vector<uint> idxs;
@@ -126,10 +123,6 @@ struct ant {
         sum_total = nextafter(sum_total, std::numeric_limits<double>::infinity());
       }
 
-      double toss = (double) rand() / (double) RAND_MAX;
-
-      // cerr << "ant " << idx << " sum_total is " << sum_total << " toss is " << toss << endl;
-
       double cumulative = 0.0;
       for (uint loc : idxs) {
         double weight = (attractiveness[loc] / sum_total);
@@ -140,11 +133,6 @@ struct ant {
         cumulative += weight;
       }
 
-      cerr << "cum " << cumulative << " n " << idxs.size() << " toss " << toss << endl;
-      sleep(1);
-
-      other++;
-      // cerr << "ant " << idx << " change statement" << endl;
       return idxs[idxs.size() - 1];
     }
   };
@@ -153,11 +141,9 @@ struct ant {
     this->update_route(next);
     this->distance_travelled += dist_matrix[this->location][next];
     this->location = next;
-    //cerr << "traversing " << next << " dist: " << this->distance_travelled << endl;
   }
 
   void run() {
-    // cerr << "started ant - # can pick: " << this->possible_locations.size() << endl;
     while (!this->possible_locations.empty()) {
       uint next = this->pick_path();
       this->traverse(next);
@@ -165,7 +151,6 @@ struct ant {
     this->distance_travelled += dist_matrix[this->route[this->route.size() - 1]]
       [this->route[0]];
 
-    // cerr << "stopped ant" << endl;
     this->tour_complete = true;
   }
 };
@@ -178,10 +163,11 @@ void *ant_thread (void *ant_ptr) {
 
 ant* ants;
 
-void init_ants() {
+void init_ants(size_t it) {
   for (size_t i = 0; i < n_ants; i++) {
-    ants[i].init();
     ants[i].idx = i;
+    ants[i].iter = it;
+    ants[i].init();
   }
 }
 
@@ -193,7 +179,7 @@ void init_aco() {
   first_pass = true;
   memset(pheromone_map, 0, sizeof(pheromone_map));
   memset(ant_updated_pheromone_map, 0, sizeof(ant_updated_pheromone_map));
-  init_ants();
+  init_ants(0);
   sh_dist = -1.0;
 }
 
@@ -220,8 +206,12 @@ void update_pheromone_map() {
 
 
 double mainloop(size_t i) {
-    srand(i);
     //cerr << "starting ants" << endl;
+#if SINGLE_CORE
+    for (uint j = 0; j < n_ants; j++) {
+      ant_thread(&(ants[j]));
+    }
+#else
     for (uint j = 0; j < n_ants; j++) {
       pthread_create(&(ants[j].t_handle), NULL, ant_thread, &(ants[j]));
     }
@@ -229,11 +219,10 @@ double mainloop(size_t i) {
     for (uint j = 0; j < n_ants; j++) {
       pthread_join(ants[j].t_handle, NULL);
     }
+#endif
     //cerr << "summing ants" << endl;
-    uint os = 0;
     unsigned long sum = 0L;
     for (uint j = 0; j < n_ants; j++) {
-      os += ants[j].other;
       populate_ant_updated_pheromone_map(ants[j]);
       sum += ants[j].distance_travelled;
 
@@ -250,7 +239,7 @@ double mainloop(size_t i) {
 
     first_pass = false;
 
-    init_ants();
+    init_ants(i+1);
 
     memset(ant_updated_pheromone_map, 0, sizeof(ant_updated_pheromone_map));
 
@@ -321,30 +310,22 @@ uint two_five_opt() {
         + dist_matrix[y1][x2] + dist_matrix[y3][y2];
 
       if (var_a < sh_dist && var_a < var_b && var_a < var_c) {
-        //cerr << "25opt(a) " << i << " " << j << ": " << var_a << endl;
         sh_dist = var_a;
         reverse(sh_route.begin() + i, sh_route.begin() + j);
         swaps++;
       } else if (var_b < sh_dist && var_b < var_a && var_b < var_c) {
-        //cerr << "25opt(b) " << i << " " << j << ": " << var_b << endl;
-        //cerr << vector<uint>(sh_route.begin() + i - 1, sh_route.begin() + j + 1) << endl;
         for (int k = i; k < j - 1; k++) {
           sh_route[k] = sh_route[k + 1];
         }
         sh_route[j - 1] = x2;
-        //cerr << vector<uint>(sh_route.begin() + i - 1, sh_route.begin() + j + 1) << endl;
 
         sh_dist = var_b;
         swaps++;
       } else if (var_c < sh_dist && var_c < var_a && var_c < var_b) {
-        //cerr << "25opt(c) " << i << " " << j << ": " << var_c << endl;
-        //cerr << vector<uint>(sh_route.begin() + i - 1, sh_route.begin() + j + 1) << endl;
         for (int k = j - 2; k >= i; k--) {
           sh_route[k + 1] = sh_route[k];
         }
         sh_route[i] = y1;
-        //cerr << vector<uint>(sh_route.begin() + i - 1, sh_route.begin() + j + 1) << endl;
-
         sh_dist = var_c;
         swaps++;
       }
@@ -357,8 +338,10 @@ uint two_five_opt() {
 
 int main(int argc, char** argv) {
   unsigned started = time(NULL);
-  if (argc < 8) {
+  int optruns;
-    cerr << argv[0] << " [n_nodes] [n_iter] [n_ants] [alpha] [beta] [evap] [weight]" << endl;
+
+  if (argc < 9) {
+    cerr << argv[0] << " [n_nodes] [n_iter] [n_ants] [alpha] [beta] [evap] [weight] [optruns]" << endl;
     return 1;
   }
 
@@ -369,6 +352,7 @@ int main(int argc, char** argv) {
   sscanf(argv[5], "%lf", &beta);
   sscanf(argv[6], "%lf", &pheromone_evaporation_coeff);
   sscanf(argv[7], "%lf", &pheromone_constant);
+  sscanf(argv[8], "%d", &optruns);
   ant ants_arr[n_ants];
   ants = ants_arr;
 
@@ -392,18 +376,7 @@ int main(int argc, char** argv) {
 
   cerr << "reading done" << endl;
 
-
-
-  /*uint nnz = 0;
-    for (int i = 0; i < n_nodes; i++) {
-    for (int j = 0; j < n_nodes - 1; j++) {
-    if (dist_matrix[i][j] != 0.0) nnz++;
-    }
-    }
-    cerr << "dist: nnz " << nnz << endl;*/
-
   init_aco();
-  // scan and parse
 
   unsigned limit = 175;
   unsigned last_time = started;
@@ -411,36 +384,25 @@ int main(int argc, char** argv) {
   for (size_t i = 0; i < n_iterations; i++) {
     double d = mainloop(i);
     unsigned now = time(NULL);
+    cerr << "iter: " << i << " - dist: " << sh_dist << " - avg: "
+         << d << " - elapsed: " << now - started << endl;
     if ((now - started) + last_delta > limit) {
         cerr << "out of time" << endl;
         break;
     }
-    cerr << "iter: " << i << " - dist: " << sh_dist << " - avg: "
-         << d << " - elapsed: " << now - started << endl;
     last_delta = now - last_time;
     last_time = now;
   }
 
-  //for (uint k = 0; k < n_nodes; k++) {
-  //  cerr << "p: " << sh_route[k] << ": " << dist_matrix[sh_route[k]][sh_route[(k + 1) % n_nodes]] << endl;
-  //}
-
   cerr << "pre-optimization length: " << sh_dist << endl;
 
-  uint i = 0;
+  cerr << sh_route << endl;
-  uint j = 0;
+
-  do {
+  for (uint j = 0; j < optruns; j++) {
-    //for (i = 0; i < 50; i++) {
+    two_five_opt();
-      //cerr << "2opt round " << i << endl;
+    three_opt(sh_route, dist_matrix, sh_dist);
-      //if(two_opt() == 0) break;
+    cerr << "optimized length " << j << " : " << sh_dist << endl;
-    //}
+  }
-    for (j = 0; j < 50; j++) {
-      //cerr << "25opt round " << j << endl;
-      if(two_five_opt() == 0) break;
-    }
-    cerr << "optimizing length: " << sh_dist << endl;
-  } while (i > 0 || j > 0);
 
-  cerr << "optimized length: " << sh_dist << endl;
   cout << sh_route << endl;
 }

opt.cc

@@ -0,0 +1,84 @@
+// vim: set ts=2 sw=2 et tw=80:
+#include "opt.h"
+
+using namespace std;
+
+typedef unsigned int uint;
+typedef tuple<uint, uint, uint> triple;
+
+ostream& operator<< (ostream& out, vector<uint> a) {
+  out << "[";
+  for (int i = 0; i < a.size() - 1; i++) {
+    out << a[i] << ",";
+  }
+  out << a[a.size() - 1] << "]";
+  return out;
+}
+
+ostream& operator<< (ostream& out, triple& a) {
+  out << "(" << get<0>(a) << "," << get<1>(a) << "," << get<2>(a) << ")";
+  return out;
+}
+
+double reverse_segment_if_better(vector<uint>& tour, const double (&dist_matrix) [1577][1577],
+    uint i, uint j, uint k) {
+
+  uint A = tour[i-1];
+  uint B = tour[i];
+  uint C = tour[j-1];
+  uint D = tour[j];
+  uint E = tour[k-1];
+  uint F = tour[k % (tour.size()-1)];
+
+  double d0 = (dist_matrix[A][B] + dist_matrix[C][D] + dist_matrix[E][F]);
+  double d1 = (dist_matrix[A][C] + dist_matrix[B][D] + dist_matrix[E][F]);
+  double d2 = (dist_matrix[A][B] + dist_matrix[C][E] + dist_matrix[D][F]);
+  double d3 = (dist_matrix[A][D] + dist_matrix[E][B] + dist_matrix[C][F]);
+  double d4 = (dist_matrix[F][B] + dist_matrix[C][D] + dist_matrix[E][A]);
+
+  if (d0 > d1) {
+    reverse(tour.begin() + i, tour.begin() + j);
+    return -d0 + d1;
+  }
+  else if (d0 > d2) {
+    reverse(tour.begin() + j, tour.begin() + k);
+    return -d0 + d2;
+  }
+  else if (d0 > d4) {
+    reverse(tour.begin() + i, tour.begin() + k);
+    return -d0 + d4;
+  }
+  else if (d0 > d3) {
+    vector<uint> tmp;
+    for (uint z = j; z < k; z++) {
+      tmp.push_back(tour[z]);
+    }
+    for (uint z = i; z < j; z++) {
+      tmp.push_back(tour[z]);
+    }
+
+
+    for (uint f = i; f < k; f++) {
+      tour[f] = tmp[f - i];
+    }
+    return -d0 + d3;
+  }
+  return 0;
+}
+
+void three_opt(vector<uint>& tour, const double (&dist_matrix) [1577][1577], double& sh_dist) {
+  double delta = 0;
+  uint n = tour.size() - 1;
+  for (uint i = 0; i < 1; i++) {
+    for (uint i = 1; i < n; i++)
+      for (uint j = i+2; j < n; j++)
+        for (uint k = j+2; k < n; k++) {
+          double d = reverse_segment_if_better(tour, dist_matrix, i, j, k);
+          if (d > 0) cerr << d << " " << i << "," << j << "," << k << endl;
+          delta += d;
+        }
+    if (delta >= 0)
+      break;
+  }
+  sh_dist += delta;
+}

opt.h

@@ -0,0 +1,11 @@
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <limits>
+#include <algorithm>
+#include <tuple>
+
+extern void three_opt(std::vector<uint>& tour,
+		const double (&dist_matrix) [1577][1577],
+		double& sh_dist);
+extern std::ostream& operator<< (std::ostream& out, std::vector<uint> a);

run.py

@@ -1,5 +1,5 @@
 import glob
-#import pandas as pd
+import pandas as pd
 from src.io_tsp import ProblemInstance
 from src.TSP_solver import TSPSolver, available_improvers, available_solvers
 import numpy as np
@@ -29,7 +29,7 @@ def use_solver_to_compute_solution(solver, improve, index, results, name, verbos
 def run(show_plots=False, verbose=False):
     problems = glob.glob('./problems/*.tsp')
 
-    problems = ["./problems/fl1577.tsp"]
+    # problems = ["./problems/fl1577.tsp"]
 
     solvers_names = available_solvers.keys()
     improvers_names = available_improvers.keys()
@@ -67,8 +67,8 @@ def run(show_plots=False, verbose=False):
 
     #index = pd.MultiIndex.from_tuples(index, names=['problem', 'method'])
 
-    return None
+    #return None
-    #return pd.DataFrame(results, index=index, columns=["tour length", "optimal solution", "gap", "time to solve"])
+    return pd.DataFrame(results, index=index, columns=["tour length", "optimal solution", "gap", "time to solve"])
 
 
 if __name__ == '__main__':

src/ant_colony.py

@@ -1,61 +1,41 @@
-from src.io_tsp import ProblemInstance
 import os
-# OOT = out-of-time
-
-# Run #0 unknown
-
-# Run #1
-# alpha = 0.5, beta = 7, evap = 0.5, weight = instance.best_sol
-# if instance.nPoints > 1000: ants = 550 loops = 15
-# elif instance.nPoints > 195: ants = 1000 loops = 55
-# else: ants = 2500 loops = 150
-#   pr439=112263 pcb442=53699 d198=16199 fl1577=<OOT> ch130=6332 u1060=247703
-#   kroA100=21737 eil76=554 rat783=<OOT> lin318=45580
-
-# Run #2
-# alpha, beta, evap, weight = (0.9, 8, 0.4, 100_000)
-# ants, loops = (900, 6) if instance.nPoints > 1100 \
-#     else (600, 15) if instance.nPoints > 1000 \
-#     else (750, 30) if instance.nPoints > 700 \
-#     else (975, 40) if instance.nPoints > 500 \
-#     else (1000, 50) if instance.nPoints > 300 \
-#     else (1100, 70) if instance.nPoints > 195 else (2700, 140)
-#   pr439=111322 pcb442=52176 d198=16177 ch130=6269 u1060=246085
-#   kroA100=21665 eil76=550 lin318=43675
-# Separate run: fl1577=24238 (132s) rat783=9389 (174s)
-
-# Run #3
-# alpha, beta, evap, weight = (0.1, 1, 0.1, 1 - 15/instance.nPoints)
-# ants, loops = (800, 75) if instance.nPoints > 1100 \
-#     else (800, 100) if instance.nPoints > 1000 \
-#     else (800, 300) if instance.nPoints > 700 \
-#     else (800, 750) if instance.nPoints > 500 \
-#     else (800, 1000) if instance.nPoints > 300 \
-#     else (800, 1250) if instance.nPoints > 195 else (800, 1500)
 
 dir = "/Users/maggicl/Git/AI2020BsC/c_prob/"
 def ant_colony_opt(instance):
     fname = dir + instance.name + ".txt"
 
     # write .mat file for C++
-    if not os.path.exists(fname):
+    with open(fname, "w") as f:
-        with open(fname, "w") as f:
+        for i in range(instance.nPoints):
-            for i in range(instance.nPoints):
+            print(" ".join(map(str, instance.dist_matrix[i])), file=f)
-                print(" ".join(map(str, instance.dist_matrix[i])), file=f)
 
-    alpha, beta, evap, weight = (0.9, 6, 0.6, 100_000)
+    # single core params
-    ants, loops = (400, 100_000) if instance.nPoints > 1100 \
+    params = {
-        else (600, 15) if instance.nPoints > 1000 \
+        "eil76": (0.9, 8, 0.4, 1000, 500, 200, 30), # 543
-        else (750, 30) if instance.nPoints > 700 \
+        "d198": (0.9, 8, 0.4, 1000, 500, 60, 30), # 16042
-        else (975, 40) if instance.nPoints > 500 \
+        "ch130": (0.9, 8, 0.4, 1000, 1750, 25, 30), # 6212
-        else (1000, 50) if instance.nPoints > 300 \
+        "kroA100": (0.9, 8, 0.4, 1000, 750, 100, 30), # 21378
-        else (1100, 70) if instance.nPoints > 195 else (2700, 140)
+        "lin318": (0.9, 8, 0.4, 1000, 500, 32, 30), # 43171
+        "pcb442": (0.9, 8, 0.4, 1000, 450, 24, 3), # 52466
+        "pr439": (0.9, 8, 0.4, 1000, 450, 24, 2), # 11734
+        "rat783": (0.9, 8, 0.4, 1000, 450, 24, 2), # 9232
+        "u1060": (0.9, 8, 0.4, 1000, 350, 6, 3), # 238025
+        "fl1577": (0.9, 8, 0.4, 1000, 50, 9, 3), # 24145
+    }
+    # params_multicore = {
+    #     "d198": (0.9, 8, 0.4, 1000, 750, 200, 30), # 15970
+    #     "ch130": (0.9, 8, 0.4, 1000, 1750, 200, 30), # 6173
+    # }
+
+    alpha, beta, evap, weight, ants, loops, optruns = params[instance.name]
 
     # Call C++ program
     cmd = dir + "aco " + str(instance.nPoints) + " " + str(loops) + " " + str(ants) + \
-        " " + str(alpha) + " " + str(beta) + " " + str(evap) + " " + str(weight) + \
+        " " + str(alpha) + " " + str(beta) + " " + str(evap) + " " + str(weight) + " " + str(optruns) + \
         " < " + fname + " &2>/dev/null"
     print(cmd)
     solution = eval(os.popen(cmd).read())
     solution.append(solution[0])
+    print(solution)
+    print(len(solution))
     return solution