see http://hamiltonianplumber.appspot.com/

PlumbersCoins.py

from sys import stdin, stdout

import random, math

def tour_length(G, xs):
  s = sum(G[xs[i-1]][xs[i]] for i in range(len(xs)))
  n = len(xs)
  return s + min(G[0][xs[i]] - max(G[xs[i]][xs[(i+1)%n]], G[xs[i]][xs[(i-1)%n]])
      for i in range(n))

def tsp(G):
  xs = list(set(G.keys()) - set([0]))
  N = len(xs)
  if N <= 1:
    return xs
  random.shuffle(xs)
  # print('init',xs)
  beta = 1.0
  n_accept = 0
  best_energy = None
  energy = tour_length(G, xs)
  for iter in xrange(100000):
      if n_accept == 100:
          beta *= 1.005
          n_accept = 0
      p = random.uniform(0.0, 1.0)
      if p < 0.2:
          i = random.randint(0, N / 2)
          xs = xs[i:] + xs[:i]
          i = random.randint(0, N / 2)
          a = xs[:i]
          a.reverse()
          new_xs = a + xs[i:]
      elif p < 0.6 and N >= 3:
          new_xs = xs[:]
          i = random.randint(1, N - 1)
          a = new_xs.pop(i)
          j = random.randint(1, N - 2)
          new_xs.insert(j, a)
      else:
          new_xs = xs[:]
          i = random.randint(1, N - 1)
          j = random.randint(1, N - 1)
          new_xs[i] = xs[j]
          new_xs[j] = xs[i]
      new_energy = tour_length(G, new_xs)
      if new_energy < energy or random.uniform(0.0, 1.0) < math.exp(- beta * (new_energy - energy)):
          n_accept += 1
          energy = new_energy
          xs = new_xs[:]
          if best_energy is None or energy < best_energy:
             best_energy = energy
             best_tour = xs[:]
             # print('best',best_tour)
      # if iter % 100000 == 0: print energy, iter, 1.0 / beta
  return best_tour[:]

def tsp_iter(G):
  ys = None
  ys_len = float('inf')
  for i in range(3):
    xs = tsp(G)
    xs_len = tour_length(G, xs)
    if xs_len < ys_len:
      ys_len = xs_len
      ys = xs
  return ys

def main():
  ts = [int(t) for t in stdin.readline().split()]
  cs = [int(c) for c in stdin.readline().split()]
  assert ts == sorted(ts)
  assert cs == sorted(cs)
  assert len(set(cs) | set(ts) | set([0])) == len(cs) + len(ts) + 1
  vs = [0] + cs + ts
  G = { u : { v : abs(u - v) for v in vs } for u in vs }
  for i in range(len(ts) / 2):
    u, v = ts[2 * i], ts[2 * i + 1]
    G[u][v] = G[v][u] = 1
  # print(G)
  for k in vs:
    for i in vs:
      for j in vs:
        G[i][j] = min(G[i][j], G[i][k] + G[k][j])
  C = { u : { v : G[u][v] for v in [0] + cs } for u in [0] + cs }
  stdout.write('{}\n'.format(tour_length(C, tsp_iter(C))))

if __name__ == '__main__':
  main()

I should add a comment here: The code is ripped off from the course "Statistical Mechanics Algorithms and Computation" by Werner Krauth et al