In [1]:
# load libraries and set plot parameters
import os
import numpy as np
import pandas as pd
import numpy as np
import scipy.spatial as sp
import scipy.stats as st
from pathlib import Path
import logging
import datetime
import pickle
logger = logging.getLogger()
def setup_file_logger(log_file):
hdlr = logging.FileHandler(log_file)
formatter = logging.Formatter('%(asctime)s %(levelname)s %(message)s')
hdlr.setFormatter(formatter)
logger.addHandler(hdlr)
logger.setLevel(logging.INFO)
def log(message):
#outputs to Jupyter console
print('{} {}'.format(datetime.datetime.now(), message))
#outputs to file
logger.info(message)
setup_file_logger('out.log')
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set_style('ticks')
sns.axes_style({'spines.right': False,
'axes.spines.top': False})
sns.set_palette(sns.color_palette("husl", 8))
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('pdf', 'png')
plt.rcParams['savefig.dpi'] = 75
plt.rcParams['figure.autolayout'] = False
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['axes.labelsize'] = 22
plt.rcParams['axes.titlesize'] = 22
plt.rcParams['xtick.labelsize'] = 20
plt.rcParams['ytick.labelsize'] = 22
plt.rcParams['font.size'] = 22
plt.rcParams['lines.linewidth'] = 1.6
plt.rcParams['lines.markersize'] = 8
plt.rcParams['legend.fontsize'] = 16
# figures:
# Where to save the figures
PROJECT_ROOT_DIR = "../"
IMAGES_PATH = os.path.join(PROJECT_ROOT_DIR, "figures")
if not os.path.exists(IMAGES_PATH):
os.makedirs(IMAGES_PATH)
def save_fig(fig_id, tight_layout=True, fig_extension="png", resolution=300):
path = os.path.join(IMAGES_PATH, fig_id + "." + fig_extension)
if tight_layout:
plt.tight_layout()
plt.savefig(path, format=fig_extension, dpi=resolution)
%load_ext autoreload
%autoreload 2
In [2]:
from tsp.tsp import TravelingSalesPersonProblem
import tsp.tsp_heuristic as th
In [3]:
def build_plot(heuristic_result):
start = heuristic_result.get_starting_node_for_plotting()
x, y = np.array(heuristic_result.get_cycle_for_plotting()).T
plt.rcParams['figure.figsize'] = (12, 8)
plt.plot(x, y)
plt.plot(start[0], start[1], "ro")
plt.title("Distance travelled: {0:.1f}".format(
heuristic_result.loss()))
plt.axis('off')
plt.show()
Traveling Salesman Problem¶
Introduction¶
The traveling salesman problem is an optimisation problem for which to date no analytical solution has been found. It describes a network of nodes that have to be connected in a tour so that the path connecting all nodes is as short as possible. Each node can only be visited once and thus has exactly two connecting edges.
There are two main families of heuristics for solving the problem. The construction heuristics start with one or a few nodes drawn at random and build the tour from there. The algorithm stops when all nodes are connected in the tour. The family of improvement heuristics start with a randomly chosen tour and then look for possible improvements using several different strategies. The algorithm stops after a certain number of executions.
In the scope of this project, the following heuristics will be showcased:
- Construction heuristics
- Best insertion
- Best-best insertion
- Shortest Edge
- Improvement heuristics
- Greedy local search with the moves:
- Swaping
- Translation
- Inversion
- Mixed move
- Simulated Annealing with the criteria
- Metropolis
- Heat bath
- Greedy local search with the moves:
Implementation¶
Data structures¶
Three classes represent the network on which heuristics operate. These are:
- Node: A vertex in the network with x- and y-coordinates and a number
- Edge: A line connecting two nodes with a length (cost/distance)
- TravelingSalesPersonProblem: Has a list of all nodes in the network, holds the distance matrix describing the distance between any two nodes in the network and has a distance metric registered used to compute the distance matrix.
In [4]:
nodes = pd.read_csv("TSP_411.txt", sep='\s+', names=['node', 'x', 'y'])
The network consists of {{len(nodes.index)}} nodes. The nodes are read in from a csv file. The three columns {{list(nodes.columns.values)}} in the file represent the node number and the coordinates of a given node. The unconnected nodes are shown in Figure .
In [5]:
nodes.plot(kind='scatter', x='x', y='y')
save_fig("network_nodes")
In [6]:
tsp = TravelingSalesPersonProblem(nodes)
Construction heuristics¶
Best Insertion¶
The best insertion heuristic starts with 3 randomly chosen nodes. It then randomly selects the next node not included in the tour and finds the pair of nodes between which insertion results in the smallest increase in tour length.
This strategy can be seen as a global optimization.
Best-Best Insertion¶
Shortest Edge¶
The shortest edge works on a list of all possible $n(n-1)$ edges between the nodes in the network. The edges are first sorted by length. Subsequently, the shortest edge is selected as the starting point. From there, the sorted list is searched for the next node meeting the two constraints:
- Inserted edge does not produce a look that is shorter than the one connecting all nodes
- The two nodes at the endpoints of the edge do not have more than 2 connected edges after insertion
The algorithm is repeated until the number of edges is equal to the number of edges connecting all nodes in the network.
Improvement heuristics¶
The family of improvement heuristics starts with a randomly constructed tour and seeks to improve the loss by performing various operations. The stopping criterion here is a fixed number of improvement operations.
Greedy Local Search¶
Starting with the random tour, randomly selected nodes are subjected to one of several possible moves. The loss with the move applied is compared to the previous loss and if smaller, the new tour is accepted.
Simulated Annealing¶
This heuristic mimics the physical process of annealing metals. At a higher temperature, locally modified cycles have a higher probability of being accepted, so we deliberately accept solutions that increase loss, albeit with a relatively small probability. The temperature is then gradually decreased, raising the bar for accepting new cycles. When the low temperature is reached, the algorithm stops.
There are two conditions for accepting new cycles:
- Metropolis: This condition accepts all cycles with an improved loss and those with a larger loss at probability [INSERT FORMULA]
- Heatbath: The heatbath condition does not automatically accept cycles with a lower loss.
Results¶
In [7]:
runs_file = Path('./runs.pkl')
if runs_file.is_file():
with open(runs_file, 'rb') as input:
runs = pickle.load(input)
else:
n_runs = 30
runs = {
'BestInsertion': [],
'BestBestInsertion': [],
'ShortestEdge': [],
'GreedyLocalSearchSwap': [],
'GreedyLocalSearchTranslate': [],
'GreedyLocalSearchInvert': [],
'GreedyLocalSearchMixed': [],
'SimulatedAnnealingMetropolis': [],
'SimulatedAnnealingHeatBath': []
}
for i in range(n_runs):
heuristics = [
th.BestInsertion(tsp),
th.BestBestInsertion(tsp),
th.ShortestEdge(tsp),
th.GreedyLocalSearchSwap(tsp),
th.GreedyLocalSearchTranslate(tsp),
th.GreedyLocalSearchInvert(tsp),
th.GreedyLocalSearchMixed(tsp),
th.SimulatedAnnealingMetropolis(tsp),
th.SimulatedAnnealingHeatBath(tsp)
]
log("Starting calculations for statistical comparison of heuristics. Doing {} runs...".format(n_runs))
for h in heuristics:
h.calculate_cycle()
runs[h.__class__.__name__].append([h.loss(), h.get_cycle()])
log("Iteration {}: Heuristic {} done".format(i+1, h.__class__.__name__))
with open('./runs.pkl', 'wb') as output:
pickle.dump(runs, output)
In [8]:
losses = {}
for (key, value) in runs.items():
l = []
for i in range(len(value)):
l.append(value[i][0])
losses[key] = l
losses_pd = pd.DataFrame.from_dict(losses)
Summary Table¶
In [9]:
stats = losses_pd.describe().drop(['count','25%','50%','75%']).apply(lambda x: round(x, 1))
In [10]:
def calc_ci(x, conf_int=0.95):
interval = st.t.interval(0.95, len(x)-1, loc=np.mean(x), scale=np.std(x, ddof=1))
return (round(interval[0], 1), round(interval[1], 1))
In [11]:
ci = pd.DataFrame(losses_pd.apply(calc_ci),columns=["CI"]).T
In [12]:
stats = pd.concat([stats,ci])
In [13]:
stats
Out[13]:
Best Solutions¶
In [14]:
hr = th.TspHeuristic(tsp)
In [15]:
import pandas as pd
In [16]:
runs_pd = pd.DataFrame.from_dict(runs)
In [17]:
min_results = runs_pd.min()
In [18]:
min_results
Out[18]:
In [19]:
def plot_resulting_cycle(heuristic, results):
x, y = np.array(hr.get_cycle_for_plotting(results[1])).T
plt.rcParams['figure.figsize'] = (12, 8)
plt.plot(x, y)
plt.rcParams['axes.titlesize'] = 42
plt.title("{algo}: $L = {loss:.1f}$".format(algo=heuristic, loss=results[0]))
plt.axis('off')
save_fig("min_"+heuristic)
plt.show()
In [20]:
for i, heuristic in enumerate(min_results.index):
plot_resulting_cycle(heuristic, min_results[i])
Pairwise Comparison¶
For the two-sample T-test, we assume equal variance in- and idependence between- the two samples.
$H_0$: The sample means are equal. $H_1$: The sample means are not equal.
The confidence level $\alpha$ is set at 5 %.
A negative statistic indicates that the 'left' sample has a lower mean while a positive statistic indicates that the 'right' sample has a lower mean. When we observe a p-value $\leq \alpha$, we reject the null hypothesis $H_0$ at $\alpha = 5\%$ and thus conclude that the difference in means is significant.
Thus, We can find out which heuristic yields a lower average loss and is thus 'better' at producing short cycles.
Construction Heuristics¶
In [21]:
ttest_b_bb = st.ttest_ind(losses['BestInsertion'], losses['BestBestInsertion'])
ttest_b_se = st.ttest_ind(losses['BestInsertion'], losses['ShortestEdge'])
sd_b = st.tstd(losses['BestInsertion'])
sd_bb = st.tstd(losses['BestBestInsertion'])
sd_se =st.tstd(losses['ShortestEdge'])
def acceptance_interval(sd_a, sd_b, size=30):
m = (((sd_a**2/size)+(sd_b**2/size))**2/(((sd_a**2/size)**2/(size+1))+((sd_b**2/size)**2/(size+1))))-2
t_stat = st.t.ppf(.975, m)
return t_stat*np.sqrt((sd_a**2/size)+(sd_b**2/size))
pairwise_construction = pd.DataFrame({'best-bestbest':
{'statistic': ttest_b_bb.statistic,
'pval': ttest_b_bb.pvalue,
'dmean': np.mean(losses['BestInsertion'])-np.mean(losses['BestBestInsertion']),
'ailow': -acceptance_interval(sd_b,sd_bb),
'aihigh': acceptance_interval(sd_b,sd_bb)},
'best-shortest':
{'statistic': ttest_b_se.statistic,
'pval': ttest_b_se.pvalue,
'dmean': np.mean(losses['BestInsertion'])-np.mean(losses['ShortestEdge']),
'ailow': -acceptance_interval(sd_bb,sd_se),
'aihigh': acceptance_interval(sd_bb,sd_se)}}).T
In [22]:
pairwise_construction
Out[22]:
In [23]:
ttest_st = st.ttest_ind(losses['GreedyLocalSearchSwap'], losses['GreedyLocalSearchTranslate'])
ttest_si = st.ttest_ind(losses['GreedyLocalSearchSwap'], losses['GreedyLocalSearchInvert'])
ttest_sm = st.ttest_ind(losses['GreedyLocalSearchSwap'], losses['GreedyLocalSearchMixed'])
ttest_ti = st.ttest_ind(losses['GreedyLocalSearchTranslate'], losses['GreedyLocalSearchInvert'])
ttest_tm = st.ttest_ind(losses['GreedyLocalSearchTranslate'], losses['GreedyLocalSearchMixed'])
ttest_im = st.ttest_ind(losses['GreedyLocalSearchInvert'], losses['GreedyLocalSearchMixed'])
ttest_mh = st.ttest_ind(losses['SimulatedAnnealingMetropolis'], losses['SimulatedAnnealingHeatBath'])
sd_gls_s = st.tstd(losses['GreedyLocalSearchSwap'])
sd_gls_t = st.tstd(losses['GreedyLocalSearchTranslate'])
sd_gls_i = st.tstd(losses['GreedyLocalSearchInvert'])
sd_gls_m = st.tstd(losses['GreedyLocalSearchMixed'])
sd_sa_m = st.tstd(losses['SimulatedAnnealingMetropolis'])
sd_sa_h = st.tstd(losses['SimulatedAnnealingHeatBath'])
pairwise_improvement = pd.DataFrame({'swap-trans':
{'statistic': ttest_st.statistic,
'pval': ttest_st.pvalue,
'dmean': np.mean(losses['GreedyLocalSearchSwap'])-np.mean(losses['GreedyLocalSearchTranslate']),
'ailow': -acceptance_interval(sd_gls_s, sd_gls_t),
'aihigh': acceptance_interval(sd_gls_s, sd_gls_t)},
'swap-invert':
{'statistic': ttest_si.statistic,
'pval': ttest_si.pvalue,
'dmean': np.mean(losses['GreedyLocalSearchSwap'])-np.mean(losses['GreedyLocalSearchInvert']),
'ailow': -acceptance_interval(sd_gls_s, sd_gls_i),
'aihigh': acceptance_interval(sd_gls_s, sd_gls_i)},
'swap-mixed':
{'statistic': ttest_sm.statistic,
'pval': ttest_sm.pvalue,
'dmean': np.mean(losses['GreedyLocalSearchSwap'])-np.mean(losses['GreedyLocalSearchMixed']),
'ailow': -acceptance_interval(sd_gls_s, sd_gls_m),
'aihigh': acceptance_interval(sd_gls_s, sd_gls_m)},
'trans-invert':
{'statistic': ttest_ti.statistic,
'pval': ttest_ti.pvalue,
'dmean': np.mean(losses['GreedyLocalSearchTranslate'])-np.mean(losses['GreedyLocalSearchInvert']),
'ailow': -acceptance_interval(sd_gls_t, sd_gls_i),
'aihigh': acceptance_interval(sd_gls_t, sd_gls_i)},
'trans-mixed':
{'statistic': ttest_tm.statistic,
'pval': ttest_tm.pvalue,
'dmean': np.mean(losses['GreedyLocalSearchTranslate'])-np.mean(losses['GreedyLocalSearchMixed']),
'ailow': -acceptance_interval(sd_gls_t, sd_gls_m),
'aihigh': acceptance_interval(sd_gls_t, sd_gls_m)},
'invert-mixed':
{'statistic': ttest_im.statistic,
'pval': ttest_im.pvalue,
'dmean': np.mean(losses['GreedyLocalSearchInvert'])-np.mean(losses['GreedyLocalSearchMixed']),
'ailow': -acceptance_interval(sd_gls_i, sd_gls_m),
'aihigh': acceptance_interval(sd_gls_i, sd_gls_m)},
'metro-heat':
{'statistic': ttest_mh.statistic,
'pval': ttest_mh.pvalue,
'dmean': np.mean(losses['SimulatedAnnealingMetropolis'])-np.mean(losses['SimulatedAnnealingHeatBath']),
'ailow': -acceptance_interval(sd_sa_m, sd_sa_h),
'aihigh': acceptance_interval(sd_sa_m, sd_sa_h)}
}).T
In [24]:
def build_comparison_plot(results, left_names, right_names, out_name):
fig = plt.figure()
length = len(results.index)
plt.rcParams['figure.figsize'] = (12, 8)
ax1 = fig.add_subplot(111)
ax1.barh(range(length), width = -results.dmean, height = 0.4, color='white',
tick_label = left_names)
ax2 = ax1.twinx() # mirror them
ax2.barh(range(length), width = (results.aihigh-results.ailow), left=results.ailow, height=0.4, tick_label = right_names, alpha=0.7)
ax2.scatter(results.dmean, range(length), marker="|", s=2000, zorder=50,linewidth=3)
ax1.tick_params(axis='both', which='major', labelsize=24,length=0)#
ax2.tick_params(axis='both', which='major', labelsize=24,length=0)
#ax1.axes.get_xaxis().set_ticks([])
#ax2.axes.get_xaxis().set_ticks([])
for spine in ax1.spines.values():
spine.set_visible(False)
for spine in ax2.spines.values():
spine.set_visible(False)
plt.xlabel('Difference in means')
plt.axvline(x=0,linewidth=1, color="black", zorder=50)
save_fig(out_name)
plt.show()
In [25]:
build_comparison_plot(pairwise_construction,
['Best Insertion', 'Best Insertion'],
['Best-best insertion', 'Shortest edge'],
"pairwise_construction")
Improvement Heuristics¶
In [26]:
np.mean(losses['GreedyLocalSearchSwap'])
Out[26]:
In [27]:
np.mean(losses['GreedyLocalSearchInvert'])
Out[27]:
In [28]:
build_comparison_plot(pairwise_improvement,
['Gredy Swap', 'Greedy Swap', 'Greedy Swap', 'Greedy Translate', 'Greedy Translate', 'Greedy Inversion', 'S.A. Metropolis'],
['Greedy Tranlsate', 'Greedy Inversion', 'Greedy Mixed', 'Greedy Inversion', 'Greedy Mixed', 'Greedy Mixed', 'S.A. Heat Bath'],
"pairwise_improvement")
Performance Plots¶
Performance plots are shown for all improvement heuristics. The loss is high at the beginning when the cycle is initialized at random. It can be expected to decrease to a lower value as more and more improvements are made. Convergence is indicated by a constant slope of the curve. The point where this leveling occurs can be interpreted as the minimum number of iterations necessary to arrive at a 'good' solution.
Greedy Local Serach Performance¶
In [29]:
performance_file = Path('./performance.pkl')
if performance_file.is_file():
with open(performance_file, 'rb') as input:
performance = pickle.load(input)
else:
performance = {}
heuristics = [
th.GreedyLocalSearchSwap(tsp),
th.GreedyLocalSearchTranslate(tsp),
th.GreedyLocalSearchInvert(tsp),
th.GreedyLocalSearchMixed(tsp),
th.SimulatedAnnealingMetropolis(tsp),
th.SimulatedAnnealingHeatBath(tsp)
]
log("Starting calculations for performance plots...")
for h in heuristics:
h.calculate_cycle(save_steps=True)
log("Heuristic {} done".format(h.__class__.__name__))
performance[h.__class__.__name__] = h.steps
with open("./performance.pkl", "wb") as output:
pickle.dump(performance, output)
In [30]:
gls_performance = pd.DataFrame({'Swap': performance['GreedyLocalSearchSwap'],
'Translation': performance['GreedyLocalSearchTranslate'],
'Inversion': performance['GreedyLocalSearchInvert'],
'Mixed': performance['GreedyLocalSearchMixed']})
In [31]:
plt.plot(gls_performance)
plt.legend(labels=['Swap','Translation', 'Inversion', 'Mixed'])
plt.xlabel('Iterations')
plt.ylabel('Loss')
save_fig("greedy_local_search_performance")
plt.show()
Simulated Annealing Performance¶
In [32]:
sa_performance_file = Path("sa_performance.pkl")
if sa_performance_file.is_file():
with open(sa_performance_file, 'rb') as input:
sa_performance = pickle.load(input)
else:
sa_performance = {}
log("Starting calculations for simulated annealing performance plots...")
for move in [th.Swap, th.Translate, th.Invert, th.Mixed]:
heuristics = [th.SimulatedAnnealingMetropolis(tsp, move=move),
th.SimulatedAnnealingHeatBath(tsp, move=move)]
log("Doing runs for move {}".format(move.__name__))
for h in heuristics:
h.calculate_cycle(save_steps=True)
log("Heuristic {} done".format(h.__class__.__name__))
sa_performance[h.__class__.__name__+move.__name__] = h.steps
with open(sa_performance_file, "wb") as output:
pickle.dump(sa_performance, output)
In [33]:
dfs = {}
for key, value in sa_performance.items():
dfs[key] = pd.DataFrame(value).T
sa_perf_df = pd.concat(dfs)
sa_perf_df['limits'] = sa_perf_df.apply(lambda row: (row[2], row[0]),axis=1)
sa_perf_df.drop(columns=['min','max'],inplace=True)
sa_perf_df.index.names = ['Heuristic', 'Temperature']
sa_perf_df.reset_index(inplace=True)
sa_perf_df['move'] = sa_perf_df.apply(lambda x: x[0].replace("SimulatedAnnealing","") ,axis=1)
In [34]:
sa_perf_df.apply(lambda x: x[0].replace("SimulatedAnnealing","") ,axis=1)
Out[34]:
In [35]:
def plot_sa_performance(performance_measurements):
for key, value in performance_measurements.items():
p = pd.DataFrame(value).T
plt.plot(p)
plt.legend(labels=['max', 'mean', 'min'])
plt.gca().invert_xaxis()
plt.title("{}".format(key))
plt.xlabel('Temperature')
plt.ylabel('Loss')
save_fig("{}".format(key))
plt.show()
In [36]:
plot_sa_performance(sa_performance)
Further attempts at minimizing loss¶
Combine BestInsertion with greedy local search¶
Combining the best construction heuristic with the best improvement heuristic
In [37]:
runs_improved_file = Path('./runs_improved.pkl')
if runs_improved_file.is_file():
with open(runs_improved_file, 'rb') as input:
runs_improved = pickle.load(input)
else:
n_runs = 30
runs_improved = []
log("Starting calculations for BestInsertion+GLS_Mixed. Doing {} runs...".format(n_runs))
for i in range(n_runs):
sh = th.BestInsertion(tsp)
gls_mixed = th.GreedyLocalSearchMixed(tsp, start_cycle_heuristic=sh)
gls_mixed.calculate_cycle()
runs_improved.append([gls_mixed.loss(), gls_mixed.get_cycle()])
log("Iteration {} done".format(i+1))
with open(runs_improved_file, 'wb') as output:
pickle.dump(runs_improved, output)
In [38]:
stats_runs_improved = pd.DataFrame(runs_improved, columns=['Loss', 'Cycle'])
In [39]:
runs_improved = pd.DataFrame(runs_improved, columns=['Loss', 'Cycle'])
In [40]:
stats = stats_runs_improved.describe().drop(['count','25%','50%','75%']).apply(lambda x: round(x, 1))
In [41]:
stats.T
Out[41]:
In [42]:
min_gls_improved = runs_improved.min()
In [43]:
x, y = np.array(hr.get_cycle_for_plotting(min_gls_improved[1])).T
plt.rcParams['figure.figsize'] = (12, 8)
plt.plot(x, y)
plt.rcParams['axes.titlesize'] = 20
plt.title("$L = {loss:.1f}$".format(loss=min_gls_improved[0]))
plt.axis('off')
save_fig("min_bi_gls_mixed")
plt.show()
Improved Simulated Annealing¶
In [44]:
sa_metropol_improved_file = Path('./sa_metropol_improved.pkl')
if sa_metropol_improved_file.is_file():
with open(sa_metropol_improved_file, 'rb') as input:
sa_metropol_improved = pickle.load(input)
else:
sa_metropol_improved = th.SimulatedAnnealingMetropolis(tsp, move=th.Mixed, max_it=1000, cooling_factor=0.999)
sa_metropol_improved.calculate_cycle(save_steps=True)
with open(sa_metropol_improved_file, 'wb') as output:
pickle.dump(sa_metropol_improved.steps, output)
In [45]:
plt.plot(sa_metropol_improved.T)
plt.legend(labels=['max', 'mean', 'min'])
plt.gca().invert_xaxis()
plt.xlabel('Temperature')
plt.ylabel('Loss')
save_fig("simulated_annealing_metropolis_improved_performance")
plt.show()
