# HG changeset patch
# User Teemu Piippo <teemu@hecknology.net>
# Date 1537902949 -10800
# Node ID c3b022f517048feadac7727e076c19085150467d
# Parent  759b0302ee7a73aebdd845623a9e397574dd28a0
optimized region computation with a blockmap

diff -r 759b0302ee7a -r c3b022f51704 compute-regions.py
--- a/compute-regions.py	Tue Sep 25 22:15:31 2018 +0300
+++ b/compute-regions.py	Tue Sep 25 22:15:49 2018 +0300
@@ -10,27 +10,60 @@
 
 regions = parse_regions(sys.argv[2])
 bus_stops = {}
+block_factor = 0.01
 
-with ZipFile(sys.argv[1]) as archive:
-	with archive.open('stops.txt') as file:
-		for row in read_csv(map(bytes.decode, file)):
-			location = Location(float(row['stop_lat']), float(row['stop_lon']))
-			reference = row['stop_id']
-			bus_stops[reference] = location
+class Blockmap:
+	def __init__(self, default = None):
+		from collections import defaultdict
+		self.blocks = default or defaultdict(set)
+	def __getitem__(self, blockid):
+		return self.blocks[blockid]
+	def blockpoint(self, point):
+		return int(point.x / block_factor), int(point.y / block_factor)
+	def block_min_corner(self, blockpoint):
+		return Location(blockpoint[0] * block_factor, blockpoint[1] * block_factor)
+	def block_max_corner(self, blockpoint):
+		return Location((blockpoint[0] + 1) * block_factor, (blockpoint[1] + 1) * block_factor)
 
-region_shapes = []
-districts = {}
-bus_stop_regions = {}
+def blocks_in_shape(blockmap, shape):
+	min_x = min(point.x for point in shape.points)
+	min_y = min(point.y for point in shape.points)
+	max_x = max(point.x for point in shape.points)
+	max_y = max(point.y for point in shape.points)
+	min_blockpoint = blockmap.blockpoint(Location(min_x, min_y))
+	max_blockpoint = blockmap.blockpoint(Location(max_x, max_y))
+	for x in range(min_blockpoint[0], max_blockpoint[0] + 1):
+		for y in range(min_blockpoint[1], max_blockpoint[1] + 1):
+			yield blockmap[x, y]
 
-for stop_id, stop_position in bus_stops.items():
+def main():
+	with ZipFile(sys.argv[1]) as archive:
+		with archive.open('stops.txt') as file:
+			for row in read_csv(map(bytes.decode, file)):
+				location = Location(float(row['stop_lat']), float(row['stop_lon']))
+				reference = row['stop_id']
+				bus_stops[reference] = location
+	region_shapes = list()
+	districts = dict()
+	bus_stop_regions = dict()
+	blocks_per_region = dict()
+	blockmap = Blockmap()
 	for region in regions.values():
-		if region['shape'].contains_point(stop_position):
-			bus_stop_regions[stop_id] = region['name']
-			break
-	else:
-		bus_stop_regions[stop_id] = None
+		blocks_per_region[region['name']] = list(blocks_in_shape(blockmap, region['shape']))
+		for block in blocks_per_region[region['name']]:
+			set.add(block, region['name'])
+	for stop_id, stop_position in bus_stops.items():
+		for region_name in blockmap[blockmap.blockpoint(stop_position)]:
+			region = regions[region_name]
+			if region['shape'].contains_point(stop_position):
+				bus_stop_regions[stop_id] = region['name']
+				break
+		else:
+			bus_stop_regions[stop_id] = None
+	covered = sum(1 if value else 0 for value in bus_stop_regions.values())
+	total = len(bus_stops)
+	print('%.1f%% bus stops covered.' % (covered * 100 / total), file = sys.stderr)
+	json.dump(bus_stop_regions, sys.stdout, indent = 2)
 
-covered = sum(1 if value else 0 for value in bus_stop_regions.values())
-total = len(bus_stops)
-print('%.1f%% bus stops covered.' % (covered * 100 / total), file = sys.stderr)
-json.dump(bus_stop_regions, sys.stdout, indent = 2)
+if __name__ == '__main__':
+	main()