Let's define our initial geometry.
A = -50 + 0j
B = 50 + 0j
M = 0 + 0j
Now we compute the dividing points of the marsh sections.
from math import sqrt
marsh_dx = 50 / sqrt(2) * 2 / 5
marsh_middle = 50 / sqrt(2)
marsh_points = [i * marsh_dx - marsh_middle for i in range(6)]
marsh_points
[-35.35533905932737, -21.213203435596423, -7.0710678118654755, 7.071067811865468, 21.21320343559642, 35.35533905932737]
We can compute the $y$ coordinate of a point $x$ on marsh line $i$ using this formula:
$$ y = x - x_m^i $$Where $x_m^i$ is a seed point from the above list.
This is the logic of the function below:
def compute_coordinate(x, line_index):
"""Returns point on line defined by line_index and with given x coordinate."""
return x + 1j * (x - marsh_points[line_index])
compute_coordinate(-35.35, 0)
(-35.35+0.005339059327369j)
In each section of the marsh, the trajectories are straight lines.
So the total time is computed as a sum of distances / speeds.
speeds = [10, 9, 8, 7, 6, 5, 10]
We can use the abs
function for computing the distances since we use complex numbers.
Let's see if we can replicate the travel time from the example.
points = [A] + [compute_coordinate(point, ind) for ind, point in enumerate(marsh_points)] + [B]
points
[(-50+0j), (-35.35533905932737+0j), (-21.213203435596423+0j), (-7.0710678118654755+0j), (7.071067811865468+0j), (21.21320343559642+0j), (35.35533905932737+0j), (50+0j)]
time = lambda points, speeds: sum(abs(from_ - to_) / speed_ for from_, to_, speed_ in zip(points[:-1], points[1:], speeds))
time(points, speeds)
13.473802361543434
distance = lambda points: sum(abs(from_ - to_) for from_, to_ in zip(points[:-1], points[1:]))
distance(points)
100.0
Ok, we're matching. So now, let's minimize this. It could definitely be interesting to do an autodiff with that.
from random import random, randint
from math import exp
random_sign = lambda : (lambda r: (r - 0.5)/abs(r - 0.5))(random())
random_sign()
1.0
def modify_points(points, index, dx):
"""Modifies points and returns copy."""
new_points = [p for p in points]
p = compute_coordinate(new_points[index].real + dx, index-1)
new_points[index] = p
return new_points
points
[(-50+0j), (-35.35533905932737+0j), (-21.213203435596423+0j), (-7.0710678118654755+0j), (7.071067811865468+0j), (21.21320343559642+0j), (35.35533905932737+0j), (50+0j)]
modify_points(points, 6, 0.01)
[(-50+0j), (-35.35533905932737+0j), (-21.213203435596423+0j), (-7.0710678118654755+0j), (7.071067811865468+0j), (21.21320343559642+0j), (35.36533905932737+0.00999999999999801j), (50+0j)]
Our improvement procedure is simple: choose a point, move it by dx and keep the arrangement if it's better than the previous one (hillclimbing, inspired by Norvig's gesture typing).
def improve(initial_points, max_dx, steps=1000):
"""Hillclimbing strategy to find best trajectory."""
points = initial_points
current_time = time(points, speeds)
for step in range(steps):
dx = random_sign() * max_dx * exp(-step / steps * 20)
new_points = modify_points(points, randint(1, 6), dx)
if time(new_points, speeds) < current_time:
current_time = time(new_points, speeds)
points = new_points
return points
Let's now run our problem:
initial_points = points
new_points = improve(initial_points, 10., steps=100000)
time(new_points, speeds)
13.126510858558497
distance(new_points)
102.13458074876188
Let's look at the optimal solution that was found:
def plot_points(points):
"""plots some points"""
plt.plot([p.real for p in points], [p.imag for p in points], '-o')
import matplotlib.pyplot as plt
%matplotlib inline
plot_points(points)
plot_points(new_points)
plt.axis('equal')
(-55.0, 55.0, -4.5411298101973614, 5.6681232353269273)
Let's round our answer to print it:
round(time(new_points, speeds), 10)
13.1265108586