Updates

de2/elements.py

@@ -4,18 +4,18 @@ import pandas as pd
 
 from dataclasses import dataclass
 
+from de2.routing import get_graph, get_closest_vertex
+
 @dataclass
 class Vessel:
 	mmsi: str
-	name: str
 	velocity: float = 0.0
-	distance_left: float = 0.0
-	total_distance: float = 0.0
-	time_until_departure: float = 0.0
 	source: tuple = None
 	target: tuple = None
 	task: str = None
-	# location: str = None
+	distance_left: float = 0.0
+	total_distance: float = 0.0
+	time_until_departure: float = 0.0
 
 	usage_time: float = 0.0
 	idle_time: float = 0.0
@@ -159,10 +159,15 @@ class Tracer(AtomicDEVS):
 	def __init__(self, name, ivef):
 		super(Tracer, self).__init__(name)
 
-		# TODO: compute the distance from the theoretical map?
+		# self.berths = {}
+		# berths = pd.read_csv("berths.csv", index_col=0)
+		# for _, row in berths.iterrows():
+		# 	self.berths[row["lpl_id"]] = row["center_lon"], row["center_lat"]
+
+		# self.ivef = pd.read_csv(ivef, usecols=["mmsi", "start", "source", "target", "task"])
 		self.ivef = pd.read_csv(ivef, usecols=["mmsi", "name", "start",
-		                                     "source_lon", "source_lat", "target_lon", "target_lat",
-		                                     "distance", "task", "velocity"],
+		                                       "source_lon", "source_lat", "target_lon", "target_lat",
+		                                       "distance", "task", "velocity"],
 		                        dtype={"start": float, "distance": float, "velocity": float,
 		                               "source_lon": float, "source_lat": float,
 		                               "target_lon": float, "target_lat": float})
@@ -171,6 +176,24 @@ class Tracer(AtomicDEVS):
 		# self.ivef["start"] -= self.starting_time
 		# self.ivef["end"] -= self.starting_time
 
+		# # obtain the coordinates
+		# self.ivef["source_lon"] = [self.berths[k][0] for k in self.ivef["source"]]
+		# self.ivef["source_lat"] = [self.berths[k][1] for k in self.ivef["source"]]
+		# self.ivef["target_lon"] = [self.berths[k][0] for k in self.ivef["target"]]
+		# self.ivef["target_lat"] = [self.berths[k][1] for k in self.ivef["target"]]
+		#
+		# # TODO: obtain distances using theoretical map
+		# self.graph = get_graph()
+		# self.ivef["distance"] = \
+		# 	self.graph.distances([get_closest_vertex(self.graph, self.berths[k][0], self.berths[k][1]).index for k in self.ivef["source"]],
+		# 	                     [get_closest_vertex(self.graph, self.berths[k][0], self.berths[k][1]).index for k in self.ivef["target"]],
+		# 	                     weights="distance", mode="all")
+		#
+		# # sample velocities
+		# self.ivef["velocity"] = 5.0
+		#
+		# self.ivef["end"] = self.ivef["start"] + self.ivef["distance"] / self.ivef["velocity"]
+
 		self.state = {
 			"index": 0,
 			"time": 0.0
@@ -233,4 +256,18 @@ class Port(CoupledDEVS):
 		for child in imm_children:
 			if isinstance(child, Pool):
 				return child
-		return imm_children[0]
+		return imm_children[0]
+
+	def print_statistics(self):
+		print("TOTAL USAGES:")
+		S = 0
+		for vessel in self.pool.state["waiting"].values():
+			if vessel.usage_time > 0:
+				print("\t%s: %f" % (vessel.mmsi, vessel.usage_time))
+				S += vessel.usage_time
+		for vessel in self.sailer.state["vessels"]:
+			if vessel.usage_time > 0:
+				print("\t%s: %f" % (vessel.mmsi, vessel.usage_time))
+				S += vessel.usage_time
+		print("-------------")
+		print("OVERALL: %f" % (S / 36))

de2/routing.py

@@ -119,6 +119,46 @@ def get_exact_vertex(graph, x, y):
     return None
 
 
+def find_percentage_point_in_path(path, distances, percentage):
+    tot_dist = sum(distances)
+    acpos = tot_dist * percentage
+    perc = 0.
+    s = 0.
+    p1 = path[0]
+    p2 = path[-1]
+    for dx, d in enumerate(distances):
+        s += d
+        if s > acpos:
+            p1 = path[dx]
+            p2 = path[dx + 1]
+            p = sum(distances[dx + 1:]) / tot_dist
+            if 1 - perc - p == 0:
+                percentage = 1
+            else:
+                percentage = (percentage - perc) / (1 - perc - p)
+            break
+        perc = s / tot_dist
+    return p1, p2, percentage
+
+
+def pathfinder(graph, source, target):
+    src, _ = get_closest_vertex(graph, source[0], source[1])
+    tgt, _ = get_closest_vertex(graph, target[0], target[1])
+    path = [source]
+    dists = []
+    if src.index != tgt.index:
+        p = graph.get_k_shortest_paths(src.index, tgt.index, 1, weights="distance", mode="all")
+        if len(p) > 0:
+            for vix in p[0]:
+                vertex = graph.vs[vix]
+                path.append((vertex["x"], vertex["y"]))
+                dists.append(haversine(path[-2][0], path[-2][1], path[-1][0], path[-1][1]))
+    path.append(target)
+    dists.append(haversine(path[-2][0], path[-2][1], path[-1][0], path[-1][1]))
+    return path, dists
+
+
+
 if __name__ == '__main__':
     graph = get_graph()
     graph.write_picklez("paths/cache.pickle.gz")

mapper.py

@@ -1,17 +1,13 @@
-import random
-
 from geoplotlib.layers import BaseLayer
 from geoplotlib.core import BatchPainter
 import geoplotlib
-from geoplotlib.utils import BoundingBox, epoch_to_str, haversine
+from geoplotlib.utils import BoundingBox, epoch_to_str
 import pyglet
 
 import pandas as pd
 import numpy as np
-import math
-import json
 
-from de2.routing import get_graph, euler_distance, get_closest_vertex, ig
+from de2.routing import get_graph, euler_distance, pathfinder, find_percentage_point_in_path as fpp
 
 
 class PoABLayer(BaseLayer):
@@ -46,23 +42,8 @@ class PoABLayer(BaseLayer):
         self.vcache = {}  # vessel -> source
         self.pcache = {}  # vessel -> path, dists
 
-        # with open("paths/endpointsWGS84.json", 'r') as f:
-        #     self.endpoints = json.load(f)["features"]
-        # for feature in self.endpoints:
-        #     loc_BD72 = feature["geometry"]["x"], feature["geometry"]["y"]
-        #     loc_WGS84 = BD72_to_WGS84(loc_BD72[0], loc_BD72[1])
-        #     feature["geometry"]["lat"] = loc_WGS84[0]
-        #     feature["geometry"]["lon"] = loc_WGS84[1]
-
         self.graph = get_graph()
 
-        # with open("paths/assenWGS84.json", 'r') as f:
-        #     axes = json.load(f)["features"]
-        # for feature in axes:
-        #     self.graph.append(feature["geometry"]["paths"][0])
-        # print(self.paths)
-        # self.idx = 0
-
     def draw_shape(self, x, y):
         """
         Draws the shape on the layer, at a specific location.
@@ -73,52 +54,27 @@ class PoABLayer(BaseLayer):
         """
         self.painter.points(x, y, 5)
 
+    def draw_path(self, path, proj):
+        """
+        Draws a path on the layer.
+
+        Args:
+            path:   The path to draw as (x, y) tuples.
+            proj:   The projector.
+        """
+        for ix in range(len(path) - 1):
+            point1 = path[ix]
+            point2 = path[ix + 1]
+            vx1, vy1 = proj.lonlat_to_screen(np.asarray([point1[0]]), np.asarray([point1[1]]))
+            vx2, vy2 = proj.lonlat_to_screen(np.asarray([point2[0]]), np.asarray([point2[1]]))
+            self.painter.lines(vx1, vy1, vx2, vy2, width=1.5)
+
     def draw(self, proj, mouse_x, mouse_y, ui_manager):
         self.painter = BatchPainter()
         tooltips = []
 
         self.painter.set_color(self.__colors__["nodes"])
 
-        # DRAW THE ENDPOINTS
-        # for feature in self.endpoints:
-        #     loc_WGS84 = feature["geometry"]["x"], feature["geometry"]["y"]
-        #     x, y = proj.lonlat_to_screen(np.asarray([loc_WGS84[0]]), np.asarray([loc_WGS84[1]]))
-        #
-        #     if euler_distance(x, y, mouse_x, mouse_y) < 10:
-        #         self.painter.set_color(self.__colors__["highlight"])
-        #         tooltips.append(feature["attributes"]["NUMMER"])
-        #     else:
-        #         self.painter.set_color(self.__colors__["nodes"])
-        #     self.draw_shape(x, y)
-
-        # # DATA ISSUE -- IGNORED
-        # for coord in [(4.24245, 51.27711), (4.24242, 51.27713), (4.24243, 51.2771)]:
-        #     x, y = proj.lonlat_to_screen(np.asarray([coord[0]]), np.asarray([coord[1]]))
-        #     self.draw_shape(x, y)
-        #
-        # poly = [[4.239617449849437, 51.275336155669656], [4.239346838137977, 51.27554571293344],
-        #         [4.242655275582035, 51.277226300904815], [4.242881992657049, 51.27701194846825],
-        #         [4.239617449849437, 51.275336155669656]]
-        # for px, p1 in enumerate(poly[:-1]):
-        #     p2 = poly[px+1]
-        #     x1, y1 = proj.lonlat_to_screen(np.asarray([p1[0]]), np.asarray([p1[1]]))
-        #     x2, y2 = proj.lonlat_to_screen(np.asarray([p2[0]]), np.asarray([p2[1]]))
-        #     self.painter.lines(x1, y1, x2, y2, width=2)
-
-        # DRAW THE PoAB GRAPH
-        # for gx, g in enumerate(self.clusters):
-        #     self.painter.set_color(self.colors[gx])
-        #     for v in g.vs:
-        #         x, y = proj.lonlat_to_screen(np.asarray([v["x"]]), np.asarray([v["y"]]))
-        #         self.draw_shape(x, y)
-        # # self.painter.set_color(self.__colors__["nodes"])
-        #     for e in g.es:
-        #         src = e.source_vertex
-        #         tgt = e.target_vertex
-        #         cx1, cy1 = proj.lonlat_to_screen(np.asarray([src["x"]]), np.asarray([src["y"]]))
-        #         cx2, cy2 = proj.lonlat_to_screen(np.asarray([tgt["x"]]), np.asarray([tgt["y"]]))
-        #         self.painter.lines(cx1, cy1, cx2, cy2, width=2)
-
         if not self.sim_active:
             ui_manager.status("Press Space to Start Simulation")
         else:
@@ -136,32 +92,11 @@ class PoABLayer(BaseLayer):
 
                     if vessel.mmsi not in self.vcache or self.vcache[vessel.mmsi] != vessel.source:
                         self.vcache[vessel.mmsi] = vessel.source
-
-                        src, _ = get_closest_vertex(self.graph, vessel.source[0], vessel.source[1])
-                        tgt, _ = get_closest_vertex(self.graph, vessel.target[0], vessel.target[1])
-                        path = [vessel.source]
-                        dists = []
-                        if src.index != tgt.index:
-                            p = self.graph.get_k_shortest_paths(src.index, tgt.index, 1, weights="distance", mode="all")
-                            if len(p) > 0:
-                                for vix in p[0]:
-                                    vertex = self.graph.vs[vix]
-                                    path.append((vertex["x"], vertex["y"]))
-                                    dists.append(haversine(path[-2][0], path[-2][1], path[-1][0], path[-1][1]))
-                        path.append(vessel.target)
-                        dists.append(haversine(path[-2][0], path[-2][1], path[-1][0], path[-1][1]))
-
-                        self.pcache[vessel.mmsi] = path, dists
+                        self.pcache[vessel.mmsi] = pathfinder(self.graph, vessel.source, vessel.target)
 
                     path, dists = self.pcache[vessel.mmsi]
-                    tot_dist = sum(dists)
 
-                    for ix in range(len(path) - 1):
-                        point1 = path[ix]
-                        point2 = path[ix+1]
-                        vx1, vy1 = proj.lonlat_to_screen(np.asarray([point1[0]]), np.asarray([point1[1]]))
-                        vx2, vy2 = proj.lonlat_to_screen(np.asarray([point2[0]]), np.asarray([point2[1]]))
-                        self.painter.lines(vx1, vy1, vx2, vy2, width=1.5)
+                    self.draw_path(path, proj)
 
                     tot_distance = vessel.total_distance
                     traveled_approx = tot_distance - vessel.distance_left
@@ -169,23 +104,11 @@ class PoABLayer(BaseLayer):
                     traveled_approx += delta * vessel.velocity
                     percentage = traveled_approx / tot_distance
 
-                    acpos = tot_dist * percentage
-                    perc = 0.
-                    s = 0.
-                    p1 = path[0]
-                    p2 = path[-1]
-                    for dx, d in enumerate(dists):
-                        s += d
-                        if s > acpos:
-                            p1 = path[dx]
-                            p2 = path[dx+1]
-                            p = sum(dists[dx+1:]) / tot_dist
-                            if 1 - perc - p == 0:
-                                percentage = 1
-                            else:
-                                percentage = (percentage - perc) / (1 - perc - p)
-                            break
-                        perc = s / tot_dist
+                    if sum(dists) == 0:
+                        p1 = p2 = path[0]
+                        percentage = 1.
+                    else:
+                        p1, p2, percentage = fpp(path, dists, percentage)
                     sx, sy = proj.lonlat_to_screen(np.asarray([p1[0]]), np.asarray([p1[1]]))
                     tx, ty = proj.lonlat_to_screen(np.asarray([p2[0]]), np.asarray([p2[1]]))
 
@@ -220,16 +143,11 @@ if __name__ == '__main__':
     sim = Simulator(port)
     sim.setClassicDEVS()
     sim.setRealTime(scale=0.0002)
-    sim.setTerminationTime(1_209_600)  # 14 days
+    sim.setTerminationTime(14 * 24 * 60 * 60)  # 14 days
 
     layer = PoABLayer(sim)
 
     geoplotlib.set_window_size(640, 760)
-    # geoplotlib.tiles_provider({
-    #     'url': lambda zoom, xtile, ytile: 'http://b.tile.stamen.com/terrain/%d/%d/%d.png' % (zoom, xtile, ytile),
-    #     'tiles_dir': 'mytiles',
-    #     'attribution': 'Map tiles by Stamen Design, under CC BY 3.0. Data @ OpenStreetMap contributors.'
-    # })
     geoplotlib.tiles_provider({
         'url': lambda zoom, xtile, ytile: 'https://tile.openstreetmap.org/%d/%d/%d.png' % (zoom, xtile, ytile),
         'tiles_dir': 'mytiles',
@@ -239,15 +157,4 @@ if __name__ == '__main__':
     geoplotlib.add_layer(layer)
     geoplotlib.show()
 
-    print("TOTAL USAGES:")
-    S = 0
-    for vessel in sim.model.pool.state["waiting"].values():
-        if vessel.usage_time > 0:
-            print("\t%s: %f" % (vessel.mmsi, vessel.usage_time))
-            S += vessel.usage_time
-    for vessel in sim.model.sailer.state["vessels"]:
-        if vessel.usage_time > 0:
-            print("\t%s: %f" % (vessel.mmsi, vessel.usage_time))
-            S += vessel.usage_time
-    print("-------------")
-    print("OVERALL: %f" % (S / 36))
+    sim.model.print_statistics()

parquet/ct2de.py

@@ -6,7 +6,7 @@ import matplotlib.path as mplPath
 from os.path import exists
 import random
 
-from de2.routing import get_graph, get_closest_vertex
+from de2.routing import get_graph, pathfinder, find_percentage_point_in_path as fpp
 
 docks = pd.read_csv("berths.csv")
 polygons = {k["lpl_id"]: eval(k["polygon"]) for _, k in docks.iterrows()}
@@ -66,26 +66,25 @@ def analyze(fname):
 	in_dock = get_dock(first["lon"], first["lat"]) != ""
 	time = first["ts"]
 	for _, row in df.iloc[1:].iterrows():
-		# See if you are at a new location
+		# Check if you are at a new location
 		new_loc = get_dock(row["lon"], row["lat"])
+
 		if new_loc != vessel.location:
 			# fix for teleportation issue
-
 			# add intermediate points on which dock/sailing transition is assumed
 			# it is assumed that these points lie on the theoretical trajectory => APPROXIMATION
-			# FIXME: a better solution would be to divide the theoretical path into equal parts instead
+
+			path, dists = pathfinder(graph, vessel.lonlat_t, (row["lon"], row["lat"]))
+			tot_dist = sum(dists)
+
 			if in_dock and new_loc != "":
-				# check movement to the next dock: divide the trajectory in 3 equal parts
-				center = (row["lon"] + 2 * vessel.lonlat_t[0]) / 3, (row["lat"] + 2 * vessel.lonlat_t[1]) / 3
-				c, _ = get_closest_vertex(graph, center[0], center[1])
-				center = c["x"], c["y"]
-				p2 = (2 * row["lon"] + vessel.lonlat_t[0]) / 3, (2 * row["lat"] + vessel.lonlat_t[1]) / 3
-				p, _ = get_closest_vertex(graph, p2[0], p2[1])
-				p2 = p["x"], p["y"]
-
-				d1 = haversine(vessel.lonlat_t[0], vessel.lonlat_t[1], center[0], center[1])
-				d2 = haversine(p2[0], p2[1], row["lon"], row["lat"])
-				d3 = haversine(center[0], center[1], p2[0], p2[1])
+				# check hidden movement to the next dock: divide the trajectory in 3 equal parts
+				pa1, pb1, percentage = fpp(path, dists, 1/3)
+				center = pa1[0] + (pb1[0] - pa1[0]) * percentage, pa1[1] + (pb1[1] - pa1[1]) * percentage
+				pa2, pb2, percentage = fpp(path, dists, 2/3)
+				p2 = pa2[0] + (pb2[0] - pa2[0]) * percentage, pa2[1] + (pb2[1] - pa2[1]) * percentage
+
+				d1 = d2 = d3 = tot_dist / 3
 				t1 = (row["ts"] + 2 * time) / 3
 				t2 = (2 * row["ts"] + time) / 3
 
@@ -93,13 +92,11 @@ def analyze(fname):
 				res.loc[len(res.index) + 1.5] = [vessel.mmsi, vessel.name, t1, t2, None,
 				                                 center[0], center[1], p2[0], p2[1], d3, vessel.task]
 			else:
-				# divide the trajectory in 2 equal parts
-				center = p2 = (row["lon"] + vessel.lonlat_t[0]) / 2, (row["lat"] + vessel.lonlat_t[1]) / 2
-				c, _ = get_closest_vertex(graph, center[0], center[1])
-				center = c["x"], c["y"]
+				# from dock to free or vice-versa: divide the trajectory in 2 equal parts
+				pa1, pb1, percentage = fpp(path, dists, 1 / 2)
+				center = p2 = pa1[0] + (pb1[0] - pa1[0]) * percentage, pa1[1] + (pb1[1] - pa1[1]) * percentage
 
-				d1 = haversine(vessel.lonlat_t[0], vessel.lonlat_t[1], center[0], center[1])
-				d2 = haversine(center[0], center[1], row["lon"], row["lat"])
+				d1 = d2 = tot_dist / 2
 				t1 = t2 = (row["ts"] + time) / 2
 
 			# store previous trajectory
@@ -127,9 +124,10 @@ def analyze(fname):
 		time = row["ts"]
 		in_dock = new_loc != ""
 
-	res.loc[len(res.index)] = [vessel.mmsi, vessel.name, vessel.start, df.iloc[-1]["ts"], vessel.location,
-	                           vessel.lonlat_s[0], vessel.lonlat_s[1], vessel.lonlat_t[0], vessel.lonlat_t[1],
-	                           vessel.distance, vessel.task]
+	if vessel.end != -1:
+		res.loc[len(res.index)] = [vessel.mmsi, vessel.name, vessel.start, df.iloc[-1]["ts"], vessel.location,
+		                           vessel.lonlat_s[0], vessel.lonlat_s[1], vessel.lonlat_t[0], vessel.lonlat_t[1],
+		                           vessel.distance, vessel.task]
 	# print("AVG dt:", sum(dts) / len(dts), max(dts), min(dts))
 	# res.sort_values(by=["start"], inplace=True)
 	res = res.sort_index().reset_index(drop=True)
@@ -140,11 +138,12 @@ if __name__ == '__main__':
 	# print(get_dock(4.242435, 51.27712))
 	# import sys
 	# analyze(sys.argv[2])
-	# analyze("results/205254890.csv")
+	# analyze("results/205257290.csv")
 
 	TFILE = "tugs.csv"
 	tugs = pd.read_csv(TFILE)
 	mmsi = tugs["MMSI"].dropna().astype(np.int32)
 	for idx in mmsi:
+		print("checking", idx)
 		analyze("results/%s.csv" % str(idx))
 		print(idx, "done")

parquet/parquet-plotter.py

@@ -16,7 +16,7 @@ for idx in mmsi:
 	if not exists(fname):
 		continue
 	tug = pd.read_csv(fname)
-	df = pd.concat([df, tug], ignore_index=True)
+	# df = pd.concat([df, tug], ignore_index=True)
 	d = tug["distance"]  # meters
 	t = (tug["end"] - tug["start"]) / 1000  # seconds
 	v = d / t
@@ -25,10 +25,10 @@ for idx in mmsi:
 	velocities = pd.concat([velocities, v], ignore_index=True)
 
 velocities.dropna(inplace=True)
-df.sort_values(by=["start"], inplace=True)
-df.to_csv("../results-de/IVEF.csv", index=False)
+# df.sort_values(by=["start"], inplace=True)
+# df.to_csv("../results-de/IVEF.csv", index=False)
 
-raise 4
+# raise 4
 
 fig, ax1 = plt.subplots()
 ax1.set_xlabel("velocity (m/s)")