tsp_olympic_park.py

"""
Find an optimal way to explore the Olympic Park west of Seattle

Ref: http://www.nps.gov/olym/planyourvisit/gettingaround.htm
"""
import z3

# import matplotlib.pyplot as plt
# import seaborn as sns
# sns.set_context("paper") # also: talk, poster, etc.

places = [
    "Elwha",
    "Fairholme",
    "Heart O'the Hills",
    "Hoh Rain Forest",
    "Hurricane Ridge",
    "Kalaloch",
    "Lake Crescent (Barnes Point)",
    "Mora",
    "Port Angeles",
    "Seattle (via ferry)",
    "Sol Duc",
    "Staircase",
    "Quinault"]

num_places = len(places)
rng = range(num_places)
map_name2idx = dict(zip(places, rng))
map_idx2name = dict(zip(rng, places))

D = [[0, 23, 18, 79, 30, 81, 16, 60, 11, 85, 39, 111, 123],
     [23, 0, 31, 57, 43, 59, 7, 36, 26, 100, 16, 126, 102],
     [18, 31, 0, 96, 12, 100, 26, 75, 5, 77, 45, 105, 133],
     [79, 57, 96, 0, 108, 40, 65, 47, 91, 156, 72, 191, 75],
     [30, 43, 12, 108, 0, 112, 38, 87, 17, 89, 57, 117, 146],
     [112, 59, 100, 40, 112, 0, 67, 49, 95, 160, 72, 195, 33],
     [16, 7, 26, 65, 38, 67, 0, 44, 21, 93, 16, 121, 108],
     [60, 36, 75, 47, 87, 49, 44, 0, 70, 137, 49, 170, 84],
     [11, 26, 5, 91, 17, 95, 21, 70, 0, 72, 40, 100, 128],
     [85, 100, 77, 156, 89, 160, 93, 137, 72, 0, 114, 110, 362],
     [39, 16, 45, 72, 57, 72, 16, 49, 40, 114, 0, 140, 132],
     [111, 126, 105, 191, 117, 195, 121, 170, 100, 110, 140, 0, 126],
     [123, 102, 133, 75, 146, 33, 108, 84, 128, 362, 132, 126, 0]]

# Construct a list of nodes
bvar = [z3.Int('bv%s' % i) for i in rng] # choice of places (0 or 1)
cvar = [z3.Int('cv%s' % i) for i in rng] # the path assignment
dvar = [z3.Int('dv%s' % i) for i in rng] # distance

solver = z3.Solver()

# Basics 
solver.add(cvar[ map_name2idx['Seattle (via ferry)'] ] == 0) # from Seattle
solver.add(bvar[ map_name2idx['Seattle (via ferry)'] ] == 1)

tot_places = bvar[0]
for i in rng[1:]:
    tot_places = tot_places + bvar[i]

# The main path constraints
for i in rng:
    solver.add( z3.Or( bvar[i] == 0, bvar[i] == 1 ) )
    solver.add( dvar[i] >= 0 )
    solver.add( z3.Or( bvar[i] == 0,
                       z3.And( bvar[i] == 1,
                               z3.Or( z3.And(cvar[i] + 1 == tot_places,                                             
                                             dvar[i] == D[i][ map_name2idx['Seattle (via ferry)'] ] ),
                                      z3.Or([ z3.And( bvar[i] == 1,
                                                      bvar[j] == 1,
                                                      cvar[j] == (cvar[i] + 1) % num_places,
                                                      dvar[i] == D[i][j] )
                                              for j in rng if j != i ]) ))))
                
# Places are assigned to meaninful ordering numbers
solver.add( z3.And([ z3.Or( bvar[i] == 0,
                            z3.And( bvar[i] == 1,
                                    z3.And( bvar[i] == 1,
                                            cvar[i] >= 0,
                                            cvar[i] < tot_places ) ))
                     for i in rng ] ))

# Places are assigned unique ordering numbers
for i in rng:
    solver.add( z3.And([ cvar[j] != cvar[i]
                         for j in rng if j != i ]) )
    
# The total distance constraint
dist = dvar[0]
for i in rng[1:]:
    dist = dist + dvar[i] * bvar[i]

##---------------------------------------------------
# Desired constraints
solver.add(bvar[ map_name2idx['Port Angeles'] ] == 1)
solver.add(bvar[ map_name2idx['Hurricane Ridge'] ] == 1)
#solver.add(bvar[ map_name2idx['Hoh Rain Forest'] ] == 1)
solver.add( tot_places > 7 )    
solver.add( dist < 300 )

# Compute the SAT problem
print(solver.check())
model = solver.model()

# Print the solution
travel_places = 0
travel_dist = 0
places_chosen = {}
for i in rng:
    if model[bvar[i]].as_long() == 0:
        continue
    travel_places += 1
    travel_dist += model[dvar[i]].as_long()
    idx = model[cvar[i]].as_long()
    #print map_idx2name[i], idx
    places_chosen[ idx ] = map_idx2name[i]

places_chosen[ len(places_chosen) ] = 'Seattle (via ferry)'

# Validate that the distances computed is correct
actual_dist = 0
for i in xrange(travel_places):
    u = map_name2idx[ places_chosen[i] ]
    v = map_name2idx[ places_chosen[i+1] ]
    actual_dist += D[u][v]

print "Total places %d, distance %d (%d)" % (travel_places, travel_dist, actual_dist)
print places_chosen