Use linear programming to find counterexample for a conjecture regarding the Traveling Salesman Problem

tsp_conjecture_counterexample.py

from scipy.optimize import linprog
import numpy as np
import itertools

#Searching for the counterexample for the following conjecture:
# in the shortest hamiltonian path, there is at least one
# node, connected to its 12 nearest neighbors
# https://www.reddit.com/r/math/comments/17rxc28/travelling_salesman_algorithm/

# Use linear programming to find a counterexample
#FIrst define the number of nodes:
N = 5

#Now the variables: they are pairwise distances between nodes.
#THere are N*(N-1)/2 of them

distance_vars = {}
idx = 0
for i in range(N):
    for j in range(i+1, N):
        distance_vars[(i,j)] = idx
        distance_vars[(j,i)] = idx
        idx += 1

M = idx
print(f"Total number of distance variables: {M}")

#Define objective as a sum of all distances, to be sure that the problem is always bounded
objective = np.array([1.0 for _ in range(M)])

#Constraints are of several sets
#First constraints that tell us that the path 0-1-2-...-(N-1)-0 is hamiltonian.
#that is, d_01 + d_12 + d_23 < all other sums

A = []
B = []

#Enumerate all transpositions of the nodes 2..(N-1)
#Construct the duration of the initial cycle
initial_cycle_coeffs = np.zeros(M)
for i in range(N):
    j = (i+1)%N
    initial_cycle_coeffs[distance_vars[i,j]] = 1
print("Constructed the function for the shortest cycle:", initial_cycle_coeffs)

#we would start enumerating all cycles from 0
for idx, perm in enumerate(itertools.permutations(list(range(1, N)))):
    #Skip the first element: it is the zero permutation, base value
    if idx == 0: continue
    #Every cycle would appear twice: in two different directions.
    #Exclude one of them.
    if perm[-1]<perm[0]: continue
    cycle = (0,) + perm
    cycle_coeffs = np.zeros(M)
    for i in range(N):
        ci = cycle[i]
        cj = cycle[(i+1)%N]
        cycle_coeffs[distance_vars[ci,cj]] = 1

    #Condition is that this cycle length is strictly bigger than the shortest cycle
    # base_cycle_len - cycle_len <= -1   
    A.append(initial_cycle_coeffs - cycle_coeffs)
    B.append(-1)

#OK, we now have conditions telling that our cycle is indeed hamiltonian    
#Now construct the counterexample conditions.
#For each node in the shortest cycle 0-1-2-...-(N-1)-0,
#(at least) one of 2 distances in the cycle
#must be bigger than one of the other N-3 distances
#This gives 2*(N-3) choices for one node, and (2*(N-3))^N choices total
#
# to make solution nondegenerate, we require that the difference is bigger than 1.0


#and add the counterexample constrints to them
#  construct choices for each node.
def node_counterexample_choices(n):
    nprev = (n-1)%N
    nnext = (n+1)%N

    #Variables for 2 lengths in the cycle
    len_forward = distance_vars[n, nnext]
    len_backward = distance_vars[n, nprev]

    #Variables for other distances, not used in the cycle
    other_lens = [distance_vars[n, (n+i)%N] for i in range(2,N-1)]
    assert len(other_lens) == N-3
    print(f"For node {n}, counterexample vars are: in_cycle={[len_forward,len_backward]}, out cycle={other_lens}")
    
    return list(itertools.product([len_forward, len_backward], other_lens))

all_counterexample_choices = [node_counterexample_choices(n) for n in range(N)]

#Define bounds for all variables
bounds = [(1.0, np.inf) for _ in range(M)]
cnt = 0
for choice in itertools.product(*all_counterexample_choices):
    #Choice is a list of pairs of 2 variables (variable for edge len in a loop) (variable for some other edge len)
    #Take the hamiltonianity constraints...
    A1 = A[:]
    B1 = B[:]
    print("##################################")
    print("Counterexample choice:", choice)
    
    #And add counterexample constraints to them
    for var_in_cycle, var_out_cycle in choice:
        assert var_in_cycle != var_out_cycle
        #Sonstraint must say that the edge length in cycle is longer than edge outside.
        # var_in_cycle >= var_out_cycle + 1.0
        # var_out_cycle - var_in_cycle <= -1.0
        a = np.zeros(M)
        a[var_in_cycle] = -1
        a[var_out_cycle] = 1
        A1.append(a)
        B1.append(-1.0)
        print(f"   constraint {a}.x <= -1")
    #Now trying to solve the problem!
    opt = linprog(c=objective,
                  A_ub = A1,
                  b_ub = B1,
                  bounds = bounds)
    if opt.success:
        print("COunterexample found!")
        for i in range(N):
            for j in range(i+1,N):
                v = distance_vars[i,j]
                print(f"  d({i},{j})={opt.x[v]}")
        break
    else:
        print("No counterexample")
    cnt += 1
else:
    print(f"Tried {cnt} constructions, no counterexample found!")
    exit(0)

D = np.zeros((N,N))
for i in range(N):
    for j in range(N):
        if j == i: continue
        v = distance_vars[i,j]
        D[i,j] = opt.x[v]
        D[j,i] = opt.x[v]
print("Distance matrix:")
print(repr(D))


print("Checking cycle lengths")

for perm in itertools.permutations(list(range(1,N))):
    cycle = (0,) + perm
    cycle_len = sum(D[cycle[i], cycle[(i+1)%N]] for i in range(N))
    print(cycle,"len=",cycle_len)