Simulation updates + datagen v2

de2/elements.py

@@ -7,13 +7,15 @@ from dataclasses import dataclass
 @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
-	location: str = None
+	task: str = None
+	# location: str = None
 
 	usage_time: float = 0.0
 	idle_time: float = 0.0
@@ -42,21 +44,31 @@ class Pool(AtomicDEVS):
 				vessel = self.state["waiting"][request["mmsi"]]
 			else:
 				print("VESSEL %s DOES NOT EXIST IN POOL" % str(request["mmsi"]))
-				vessel = Vessel(request["mmsi"])
+				vessel = Vessel(request["mmsi"], request["name"])
 
 			vessel.idle_time += self.state["time"] - vessel.time_until_departure
 
-			vessel.time_until_departure = (request["end"] - request["start"]) / 1000
-			vessel.velocity = request["distance"] / vessel.time_until_departure
-			vessel.distance_left = request["distance"]
-			vessel.total_distance = request["distance"]
-			vessel.source = (request["source_lon"], request["source_lat"])
-			vessel.target = (request["target_lon"], request["target_lat"])
-			vessel.location = request["location"]
-
-			self.state["should_exit"].append(vessel)
-			if vessel.mmsi in self.state["waiting"]:
-				del self.state["waiting"][vessel.mmsi]
+			# vessel.time_until_departure = (request["end"] - request["start"]) / 1000
+			# vessel.velocity = request["distance"] / vessel.time_until_departure
+			# vessel.distance_left = request["distance"]
+			# vessel.total_distance = request["distance"]
+			# vessel.source = (request["source_lon"], request["source_lat"])
+			# vessel.target = (request["target_lon"], request["target_lat"])
+			# vessel.location = request["location"]
+
+			vessel.velocity = request["velocity"]
+			if vessel.velocity > 0:
+				# requests for idle times should be ignored!
+				vessel.distance_left = request["distance"]
+				vessel.total_distance = vessel.distance_left
+				vessel.time_until_departure = vessel.distance_left / vessel.velocity
+				vessel.source = (request["source_lon"], request["source_lat"])
+				vessel.target = (request["target_lon"], request["target_lat"])
+				vessel.task = request["task"]
+
+				self.state["should_exit"].append(vessel)
+				if vessel.mmsi in self.state["waiting"]:
+					del self.state["waiting"][vessel.mmsi]
 		if self.vessel_in in inputs:
 			vessel = inputs[self.vessel_in]
 			if vessel.mmsi in self.state["waiting"]:
@@ -65,7 +77,8 @@ class Pool(AtomicDEVS):
 			else:
 				vessel.source = vessel.target
 				vessel.target = None
-				vessel.usage_time += vessel.total_distance / vessel.velocity
+				if vessel.velocity > 0:
+					vessel.usage_time += vessel.total_distance / vessel.velocity
 
 				# temporarily use this variable to store the "idle from" time
 				vessel.time_until_departure = self.state["time"]
@@ -146,11 +159,17 @@ class Tracer(AtomicDEVS):
 	def __init__(self, name, ivef):
 		super(Tracer, self).__init__(name)
 
-		self.ivef = pd.read_csv(ivef)
+		# TODO: compute the distance from the theoretical map?
+		self.ivef = pd.read_csv(ivef, usecols=["mmsi", "name", "start",
+		                                     "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})
 		self.ivef.sort_values(by=["start"])
 		self.starting_time = self.ivef.iloc[0]["start"]
-		self.ivef["start"] -= self.starting_time
-		self.ivef["end"] -= self.starting_time
+		# self.ivef["start"] -= self.starting_time
+		# self.ivef["end"] -= self.starting_time
 
 		self.state = {
 			"index": 0,
@@ -161,7 +180,8 @@ class Tracer(AtomicDEVS):
 
 	def timeAdvance(self):
 		if self.state["index"] < len(self.ivef):
-			return round(self.ivef.iloc[self.state["index"]]["start"] / 1000 - self.state["time"], 6)
+			start = self.ivef.iloc[self.state["index"]]["start"] - self.starting_time
+			return round(start / 1000 - self.state["time"], 6)
 		return INFINITY
 
 	def intTransition(self):

mapper.py

@@ -91,6 +91,19 @@ class PoABLayer(BaseLayer):
         #         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):
@@ -116,7 +129,7 @@ class PoABLayer(BaseLayer):
 
             for vessel in model.sailer.state["vessels"]:
                 if vessel.source is not None or vessel.target is not None:
-                    if not isinstance(vessel.location, str) and math.isnan(vessel.location):
+                    if vessel.task == "sailing":
                         self.painter.set_color(self.__colors__["sailing"])
                     else:
                         self.painter.set_color(self.__colors__["tugging"])
@@ -180,7 +193,7 @@ class PoABLayer(BaseLayer):
                     dy = (ty - sy) * percentage
                     if euler_distance(sx + dx, sy + dy, mouse_x, mouse_y) < 10:
                         self.painter.set_color(self.__colors__["highlight"])
-                        tooltips.append(str(vessel.mmsi))
+                        tooltips.append("%s (MMSI: %s)" % (str(vessel.name), str(vessel.mmsi)))
                         tooltips.append("%.3f m/s" % vessel.velocity)
                     self.draw_shape(sx + dx, sy + dy)
 
@@ -212,8 +225,13 @@ if __name__ == '__main__':
     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: 'http://a.tile.stamen.com/terrain/%d/%d/%d.png' % (zoom, xtile, ytile),
+        'url': lambda zoom, xtile, ytile: 'https://tile.openstreetmap.org/%d/%d/%d.png' % (zoom, xtile, ytile),
         'tiles_dir': 'mytiles',
         'attribution': 'Map tiles by Stamen Design, under CC BY 3.0. Data @ OpenStreetMap contributors.'
     })

parquet/ct2de.py

@@ -4,26 +4,36 @@ from dataclasses import dataclass
 from geoplotlib.utils import haversine
 import matplotlib.path as mplPath
 from os.path import exists
+import random
 
-docks = pd.read_csv("../berths.csv")
+from de2.routing import get_graph, get_closest_vertex
+
+docks = pd.read_csv("berths.csv")
 polygons = {k["lpl_id"]: eval(k["polygon"]) for _, k in docks.iterrows()}
 
+tugs = pd.read_excel("20230405_Tugs.xlsx", dtype={"MMSI": str, "NAME": str})
+tmap = {}
+for _, row in tugs.iterrows():
+	tmap[row["MMSI"]] = row["NAME"]
+
 @dataclass
 class Vessel:
 	mmsi: str
+	name: str
 	start: int
 	end: int
 	location: str
-	latlon_s: tuple
-	latlon_t: tuple
+	lonlat_s: tuple
+	lonlat_t: tuple
 	distance: float
+	task: str
 
 def is_inside(lon, lat, polygon):
 	ppath = mplPath.Path(np.array(polygon))
 	return ppath.contains_point((lon, lat))
 
 
-def get_dock(lat, lon):
+def get_dock(lon, lat):
 	# loop through all dock indexes and do insideness checks
 	for pid, poly in polygons.items():
 		if is_inside(lon, lat, poly):
@@ -37,49 +47,104 @@ def analyze(fname):
 	df = pd.read_csv(fname, dtype={"mmsi": str})
 	df.sort_values(by=["ts"], inplace=True)
 
-	res = pd.DataFrame(columns=["mmsi", "start", "end", "location", "source_lon", "source_lat", "target_lon", "target_lat", "distance"])
+	res = pd.DataFrame(columns=["mmsi", "name", "start", "end", "location", "source_lon", "source_lat",
+	                            "target_lon", "target_lat", "distance", "task"])
+
+	graph = get_graph()
+
+	task = random.choice([0, 1])
+	tasks = ["tugging", "sailing"]
 
 	first = df.iloc[0]
-	vessel = Vessel(first["mmsi"], first["ts"], -1, get_dock(first["lat"], first["lon"]),
-	                (first["lat"], first["lon"]), (first["lat"], first["lon"]), 0.0)
+	vessel = Vessel(first["mmsi"], tmap[str(first["mmsi"])], first["ts"], -1, get_dock(first["lon"], first["lat"]),
+	                (first["lon"], first["lat"]), (first["lon"], first["lat"]), 0.0, tasks[task])
 
-	if exists("results-de/%s-de.csv" % str(vessel.mmsi)):
+	if exists("../results-de/%s-de.csv" % str(vessel.mmsi)):
 		print("Already done %s" % str(vessel.mmsi))
 		return None
-	# time = first["ts"]
-	# dts = []
+
+	in_dock = get_dock(first["lon"], first["lat"]) != ""
+	time = first["ts"]
 	for _, row in df.iloc[1:].iterrows():
-		# find the distance between the previous point and the current point
-		vessel.distance += haversine(vessel.latlon_t[1], vessel.latlon_t[0], row["lon"], row["lat"])
-		vessel.latlon_t = row["lat"], row["lon"]
-		new_loc = get_dock(row["lat"], row["lon"])
-		# dts.append(row["ts"] - time)
-		# time = row["ts"]
+		# See if you are at a new location
+		new_loc = get_dock(row["lon"], row["lat"])
 		if new_loc != vessel.location:
-			# new location!
-			vessel.end = row["ts"]
-			res.loc[len(res.index)] = [vessel.mmsi, vessel.start, vessel.end, vessel.location,
-			                           vessel.latlon_s[1], vessel.latlon_s[0], vessel.latlon_t[1], vessel.latlon_t[0],
-			                           vessel.distance]
-			vessel.start = row["ts"]
+			# 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
+			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])
+				t1 = (row["ts"] + 2 * time) / 3
+				t2 = (2 * row["ts"] + time) / 3
+
+				# add a movement to the next dock
+				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"]
+
+				d1 = haversine(vessel.lonlat_t[0], vessel.lonlat_t[1], center[0], center[1])
+				d2 = haversine(center[0], center[1], row["lon"], row["lat"])
+				t1 = t2 = (row["ts"] + time) / 2
+
+			# store previous trajectory
+			res.loc[len(res.index)] = [vessel.mmsi, vessel.name, vessel.start, t1, vessel.location,
+			                           vessel.lonlat_s[0], vessel.lonlat_s[1],
+			                           center[0], center[1], vessel.distance + d1, vessel.task]
+
+			if new_loc != "":
+				# Randomly pick tugging or sailing
+				task = random.choice([0, 1])
+				vessel.task = tasks[task]
+
+			# store beginning of new trajectory
+			vessel.start = t2
 			vessel.end = -1
 			vessel.location = new_loc
-			vessel.distance = 0.0
-			vessel.latlon_s = vessel.latlon_t
-	res.loc[len(res.index)] = [vessel.mmsi, vessel.start, df.iloc[-1]["ts"], vessel.location,
-	                           vessel.latlon_s[1], vessel.latlon_s[0], vessel.latlon_t[1], vessel.latlon_t[0],
-	                           vessel.distance]
+			vessel.distance = d2
+			vessel.lonlat_s = p2
+			vessel.lonlat_t = row["lon"], row["lat"]
+		else:
+			# find the distance between the previous point and the current point
+			vessel.distance += haversine(vessel.lonlat_t[0], vessel.lonlat_t[1], row["lon"], row["lat"])
+			vessel.lonlat_t = row["lon"], row["lat"]
+
+		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]
 	# print("AVG dt:", sum(dts) / len(dts), max(dts), min(dts))
-	res.to_csv("../results-de/%s-de.csv" % str(vessel.mmsi))
+	# res.sort_values(by=["start"], inplace=True)
+	res = res.sort_index().reset_index(drop=True)
+	res["velocity"] = res["distance"] / ((res["end"] - res["start"]) / 1000)
+	res.to_csv("results-de/%s-de.csv" % str(vessel.mmsi), index=False)
 
 if __name__ == '__main__':
+	# print(get_dock(4.242435, 51.27712))
 	# import sys
 	# analyze(sys.argv[2])
 	# analyze("results/205254890.csv")
 
-	TFILE = "../tugs.csv"
+	TFILE = "tugs.csv"
 	tugs = pd.read_csv(TFILE)
 	mmsi = tugs["MMSI"].dropna().astype(np.int32)
 	for idx in mmsi:
-		analyze("../results/%s.csv" % str(idx))
+		analyze("results/%s.csv" % str(idx))
 		print(idx, "done")

parquet/parquet-plotter.py

@@ -15,7 +15,7 @@ for idx in mmsi:
 	fname = "../results-de/%s-de.csv" % str(idx)
 	if not exists(fname):
 		continue
-	tug = pd.read_csv(fname, index_col=0)
+	tug = pd.read_csv(fname)
 	df = pd.concat([df, tug], ignore_index=True)
 	d = tug["distance"]  # meters
 	t = (tug["end"] - tug["start"]) / 1000  # seconds