Inspire a community of professionals skilled in spatiotemporal analysis for community health and safety.
Draft algorithm to designate which poles should be fitted with a street lamp.
Roy Hyunjin Han
We want to make sure that each customer is near a street lamp.
20171111-0330 -
+ Draft algorithm
Press SHIFT-ENTER to run each cell.
20171111-0330 - 20171111-0400: 30 minutes
Since we can only put a street lamp on a pole and the number of poles is finite, this could be solved combinatorially.
from itertools import product
items = 0, 1
list(product(items, repeat=4))
For each customer, get the distance to the nearest street lamp. Sum these distances over all customers. This is a simple application of kdtree. We'll use this as our metric function and then simply choose the combination above that minimizes our metric function.
+ Minimize the summed distance of a customer to a street lamp over all customers
+ Place a street lamp on selected poles to satisfy the above constraint
Given a set of customer locations and street lamp locations, compute the metric function.
20171117-1700 - 20171117-1800: 60 minutes
import numpy as np
mask = [True, True, False, False, True]
np.array([1, 2, 3, 4, 5])[mask]
# Represent each pole by a boolean
# Generate boolean masks
# For each boolean mask, compute metric
# Choose boolean mask with best metric
import numpy as np
from itertools import product
def compute_street_lamp_metric(customers, poles, mask):
lamp_poles = np.array(poles)[list(mask)]
if not len(lamp_poles):
return np.inf
metric = 0
for pole in lamp_poles:
pole_point = pole.geometry
metric += sum(pole_point.distance(
x.geometry) for x in customers)
return metric
"""
def place_street_lamps(customers, poles):
best_metric = 0
best_mask = None
for mask in product((0, 1), repeat=len(poles)):
metric = compute_street_lamp_metric(customers, poles, mask)
if metric < best_metric:
best_metric = metric
best_mask = mask
for pole in np.array(poles)[list(mask)]:
pole.has_street_lamp = True
return poles
place_street_lamps(customers, poles)
""";
20171117-1800 - 20171117-1830: 30 minutes
Unfortunately, there is a problem with the above approach. The number of combinations is too large. Perhaps we can just use the same pole finding algorithm for the lamps and then assign each lamp to the nearest pole?
Another strategy we could use is hill climbing, where we place one and then place the next. It is not optimal, but it will provide an acceptable solution.
(10, 1) > (0, 10)
from scipy.spatial import KDTree
tree = KDTree([(0, 0)])
tree.query((1, 1), distance_upper_bound=0.5)
import attr
from shapely.geometry import Point
@attr.s
class Customer(object):
geometry = attr.ib()
@attr.s
class Pole(object):
geometry = attr.ib()
def get_geometries(xs):
return [x.geometry for x in xs]
customers = [
Customer(Point(0, 0)),
Customer(Point(1, 0)),
Customer(Point(1, 1)),
Customer(Point(9, 9)),
]
poles = [
Pole(Point(0.5, 0.5)),
Pole(Point(0.75, 0.75)),
Pole(Point(8, 8)),
]
maximum_street_lamp_distance_in_meters = 5
from geotable import ColorfulGeometryCollection
from shapely.geometry import GeometryCollection
ColorfulGeometryCollection([
GeometryCollection(get_geometries(customers)),
GeometryCollection(get_geometries(poles)),
], colors=['yellow', 'red'])
def sum_distances(point):
return sum(point.distance(x.geometry) for x in customers)
sum_distances(poles[0].geometry)
from scipy.spatial import KDTree
pole = poles[0]
pole_point = pole.geometry
pole_xy = pole_point.x, pole_point.y
customer_tree = KDTree([(_.x, _.y) for _ in get_geometries(customers)])
indices = customer_tree.query_ball_point(
pole_xy, maximum_street_lamp_distance_in_meters)
indices
customer_count = len(customers)
distances, indices = customer_tree.query(
pole_xy,
k=customer_count,
distance_upper_bound=maximum_street_lamp_distance_in_meters)
print(distances)
print(indices)
len(customers), len(customer_tree.data), type(distances)
import numpy as np
[x for x in distances if x < np.inf]
[x for x in indices if x < customer_count]
20171129-1300 - 20171129-1315: 15 minutes
We assembled a solution in the Assemble Algorithm (1st Implementation) notebook, which we will duplicate here. The approach is not optimal, but it works.
import numpy as np
def get_best_pole(
customer_tree, poles, maximum_lamp_pole_distance_in_meters):
best_metric = np.inf, np.inf
best_pole = None
for pole in poles:
metric = compute_lamp_pole_metric(
customer_tree, pole,
maximum_lamp_pole_distance_in_meters)
if metric < best_metric:
best_metric = metric
best_pole = pole
return best_pole
def compute_lamp_pole_metric(
customer_tree, pole, maximum_lamp_pole_distance_in_meters):
customer_count = len(customer_tree.data)
pole_point = pole.geometry
pole_xy = pole_point.x, pole_point.y
distances, indices = customer_tree.query(
pole_xy,
k=customer_count,
distance_upper_bound=maximum_lamp_pole_distance_in_meters)
if not hasattr(distances, '__iter__'):
distances = [distances]
indices = [indices]
distances = [x for x in distances if x < np.inf]
indices = [x for x in indices if x != customer_count]
return -len(indices), sum(distances)
from copy import copy
from scipy.spatial import KDTree
def get_lamp_poles(
customers, poles, maximum_lamp_pole_distance_in_meters):
remaining_customers = copy(customers)
remaining_poles = copy(poles)
lamp_poles = []
while remaining_customers:
customer_tree = KDTree([(_.x, _.y) for _ in get_geometries(
remaining_customers)])
pole = get_best_pole(
customer_tree,
remaining_poles,
maximum_lamp_pole_distance_in_meters)
lamp_poles.append(pole)
# Remove lamp pole from remaining_poles
remaining_poles.remove(pole)
# Remove satisfied customers from remaining_customers
pole_point = pole.geometry
pole_xy = pole_point.x, pole_point.y
indices = customer_tree.query_ball_point(
pole_xy,
maximum_lamp_pole_distance_in_meters)
for customer in np.array(remaining_customers)[indices]:
remaining_customers.remove(customer)
for pole in lamp_poles:
pole.has_street_lamp = True
return lamp_poles
maximum_lamp_pole_distance_in_meters = 5
lamp_poles = get_lamp_poles(
customers, poles, maximum_lamp_pole_distance_in_meters)
len(lamp_poles)
from geotable import ColorfulGeometryCollection
from shapely.geometry import GeometryCollection
ColorfulGeometryCollection([
GeometryCollection(get_geometries(customers)),
GeometryCollection(get_geometries(lamp_poles)),
], colors=['yellow', 'magenta'])