Initial commit

.gitignore

@@ -0,0 +1,2 @@
+*.parquet
+*.csv

.idea/.gitignore

@@ -0,0 +1,8 @@
+# Default ignored files
+/shelf/
+/workspace.xml
+# Datasource local storage ignored files
+/dataSources/
+/dataSources.local.xml
+# Editor-based HTTP Client requests
+/httpRequests/

.idea/PoAB.iml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<module type="PYTHON_MODULE" version="4">
+  <component name="NewModuleRootManager">
+    <content url="file://$MODULE_DIR$" />
+    <orderEntry type="inheritedJdk" />
+    <orderEntry type="sourceFolder" forTests="false" />
+  </component>
+</module>

.idea/inspectionProfiles/Project_Default.xml

@@ -0,0 +1,34 @@
+<component name="InspectionProjectProfileManager">
+  <profile version="1.0">
+    <option name="myName" value="Project Default" />
+    <inspection_tool class="PyPackageRequirementsInspection" enabled="true" level="WARNING" enabled_by_default="true">
+      <option name="ignoredPackages">
+        <value>
+          <list size="3">
+            <item index="0" class="java.lang.String" itemvalue="tqdm" />
+            <item index="1" class="java.lang.String" itemvalue="matplotlib" />
+            <item index="2" class="java.lang.String" itemvalue="bokeh" />
+          </list>
+        </value>
+      </option>
+    </inspection_tool>
+    <inspection_tool class="PyPep8Inspection" enabled="false" level="WEAK WARNING" enabled_by_default="false">
+      <option name="ignoredErrors">
+        <list>
+          <option value="E128" />
+          <option value="E701" />
+          <option value="E101" />
+        </list>
+      </option>
+    </inspection_tool>
+    <inspection_tool class="PyPep8NamingInspection" enabled="true" level="WEAK WARNING" enabled_by_default="true">
+      <option name="ignoredErrors">
+        <list>
+          <option value="N802" />
+          <option value="N806" />
+          <option value="N803" />
+        </list>
+      </option>
+    </inspection_tool>
+  </profile>
+</component>

.idea/inspectionProfiles/profiles_settings.xml

@@ -0,0 +1,6 @@
+<component name="InspectionProjectProfileManager">
+  <settings>
+    <option name="USE_PROJECT_PROFILE" value="false" />
+    <version value="1.0" />
+  </settings>
+</component>

.idea/misc.xml

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8" project-jdk-type="Python SDK" />
+</project>

.idea/modules.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="ProjectModuleManager">
+    <modules>
+      <module fileurl="file://$PROJECT_DIR$/.idea/PoAB.iml" filepath="$PROJECT_DIR$/.idea/PoAB.iml" />
+    </modules>
+  </component>
+</project>

.idea/other.xml

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project version="4">
+  <component name="PySciProjectComponent">
+    <option name="PY_SCI_VIEW_SUGGESTED" value="true" />
+    <option name="PY_MATPLOTLIB_IN_TOOLWINDOW" value="false" />
+  </component>
+</project>

Port.py

@@ -0,0 +1,65 @@
+from pypdevs.DEVS import CoupledDEVS
+from atomic import *
+from events import Ship
+import pandas as pd
+
+class Port(CoupledDEVS):
+	def __init__(self, name, plan, layer=None):
+		super(Port, self).__init__(name)
+
+		self.planner = self.addSubModel(Planner("planner", plan, layer))
+		self.pool = self.addSubModel(Pool("pool"))
+		self.sail = self.addSubModel(Sailer("sail", 6.235, 4.165)) # TODO: use the values from the data
+		self.tugs = self.addSubModel(Sailer("tugs", 4.455, 2.001))
+
+		self.connectPorts(self.planner.reqs, self.pool.req_in)
+		self.connectPorts(self.pool.ship_out_sail, self.sail.ship_in)
+		self.connectPorts(self.pool.ship_out_tug, self.tugs.ship_in)
+		self.connectPorts(self.sail.ship_out, self.pool.ship_in)
+		self.connectPorts(self.tugs.ship_out, self.pool.ship_in)
+
+	def select(self, imm_children):
+		for child in imm_children:
+			if isinstance(child, Sailer):
+				return child
+		for child in imm_children:
+			if isinstance(child, Pool):
+				return child
+		return imm_children[0]
+
+
+if __name__ == '__main__':
+	from mapper import PoABLayer
+
+	data = pd.read_csv("2022.csv")
+	data = data.sort_values('time')
+
+	layer = PoABLayer()
+
+	port = Port("port", data, layer)
+	set = {}
+	for ix, row in data.iterrows():
+		if row["Sleepboot"] not in set:
+			set[row["Sleepboot"]] = row["Locatie van"]
+	for b, l in set.items():
+		port.pool.state["ships"][b] = Ship(b, l, 0, 0, 0, 0)
+
+	import geoplotlib
+	geoplotlib.set_window_size(640, 760)
+	geoplotlib.tiles_provider({
+		'url': lambda zoom, xtile, ytile: 'http://a.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.add_layer(layer)
+
+	from pypdevs.simulator import Simulator
+	sim = Simulator(port)
+	sim.setClassicDEVS()
+	sim.setRealTime()
+	# sim.setVerbose(None)
+	sim.setTerminationTime(1000)
+	sim.simulate()
+
+	geoplotlib.show()

README.md

@@ -0,0 +1 @@
+# PoAB

atomic.py

@@ -0,0 +1,172 @@
+import copy
+
+from pypdevs.DEVS import AtomicDEVS
+from pypdevs.infinity import INFINITY
+from events import Request
+import numpy as np
+
+
+class Pool(AtomicDEVS):
+	def __init__(self, name):
+		super(Pool, self).__init__(name)
+
+		self.state = {
+			"ships": {},
+			"requests": [],
+			"delta": 0.,
+			"time": 0.
+		}
+
+		self.req_in = self.addInPort("req_in")
+		self.ship_in = self.addInPort("ship_in")
+		self.ship_out_sail = self.addOutPort("ship_out_sail")
+		self.ship_out_tug = self.addOutPort("ship_out_tug")
+
+	def timeAdvance(self):
+		if len(self.state["requests"]) > 0:
+			return self.state["delta"]
+		return INFINITY
+
+	def extTransition(self, inputs):
+		self.state["time"] += self.elapsed
+		if self.req_in in inputs:
+			self.state["requests"].append(inputs[self.req_in])
+		if self.ship_in in inputs:
+			ship = inputs[self.ship_in]
+			if ship.id in self.state["ships"]:
+				raise ValueError("Impossible entry to pool.")
+			self.state["ships"][ship.id] = ship
+		return self.state
+
+	def intTransition(self):
+		self.state["time"] += self.timeAdvance()
+		if len(self.state["requests"]) > 0:
+			if self.state["requests"][0].ship in self.state["ships"]:
+				req = self.state["requests"].pop(0)
+				del self.state["ships"][req.ship]
+				self.state["delta"] = 0
+			else:
+				self.state["delta"] = 0.001
+		return self.state
+
+	def outputFnc(self):
+		if len(self.state["requests"]) > 0:
+			request = self.state["requests"][0]
+			if request.ship not in self.state["ships"]:
+				print("Cannot find free ship", request.ship)
+				return {}
+			ship = self.state["ships"][request.ship]
+			if ship.position != request.source:
+				print("[%s] Ship %s is not at position %s, but at %s" % (self.state["time"], request.ship, request.source, ship.position))
+			ship.x = request.distance
+			ship.v = request.speed
+			ship.position = request.destination
+			if request.tugging:
+				return {
+					self.ship_out_tug: ship
+				}
+			return {
+				self.ship_out_sail: ship
+			}
+		return {}
+
+
+class Sailer(AtomicDEVS):
+	def __init__(self, name, alpha, beta):
+		super(Sailer, self).__init__(name)
+
+		self.alpha = alpha
+		self.beta = beta
+
+		self.state = {
+			"time": 0,
+			"ships": []
+		}
+
+		self.ship_in = self.addInPort("ship_in")
+		self.ship_out = self.addOutPort("ship_out")
+
+	def timeAdvance(self):
+		if len(self.state["ships"]) > 0:
+			return self.state["ships"][0].time_until_departure
+		return INFINITY
+
+	def update_ships(self, elapsed):
+		for ship in self.state["ships"]:
+			ship.x -= ship.v * elapsed
+			ship.total_distance += ship.v * elapsed
+			ship.time_until_departure = round(ship.time_until_departure - elapsed, 6)
+			ship.time_used += elapsed
+
+	def extTransition(self, inputs):
+		self.state["time"] += self.elapsed
+		self.update_ships(self.elapsed)
+		if self.ship_in in inputs:
+			ship = inputs[self.ship_in]
+			# ship.v = np.random.gamma(self.alpha, self.beta)
+			ship.time_start = self.state["time"]
+			ship.time_until_departure = round(ship.x / ship.v, 6)
+			self.state["ships"].append(ship)
+		self.state["ships"] = sorted(self.state["ships"], key=lambda s: s.time_until_departure)
+		return self.state
+
+	def intTransition(self):
+		elapsed = self.timeAdvance()
+		self.state["time"] += elapsed
+		self.state["ships"].pop(0)
+		self.update_ships(elapsed)
+		self.state["ships"] = sorted(self.state["ships"], key=lambda s: s.time_until_departure)
+		return self.state
+
+	def outputFnc(self):
+		if len(self.state["ships"]) > 0:
+			elapsed = self.timeAdvance()
+			ship = copy.deepcopy(self.state["ships"][0])
+			ship.total_distance += ship.v * elapsed
+			ship.time_until_departure -= elapsed
+			ship.time_used += elapsed
+
+			return {
+				self.ship_out: ship
+			}
+		return {}
+
+
+class Planner(AtomicDEVS):
+	def __init__(self, name, plan, layer=None):
+		super(Planner, self).__init__(name)
+
+		self.plan = plan
+		self.state = {
+			"time": 0,
+			"delta": 0,
+			"index": 0,
+			"layer": layer
+		}
+		self.reqs = self.addOutPort("reqs")
+
+	def timeAdvance(self):
+		if self.state["index"] > len(self.plan):
+			return INFINITY
+		return self.state["delta"]
+
+	def intTransition(self):
+		self.state["time"] += self.timeAdvance()
+
+		self.state["layer"].paths = [p for p in self.state["layer"].paths if p["time_end"] < self.state["time"]]
+		record = self.plan.iloc[self.state["index"]]
+		if self.state["layer"] is not None:
+			self.state["layer"].paths.append(record)
+
+		self.state["index"] += 1
+		self.state["delta"] = self.plan.iloc[self.state["index"]]["time"] - self.state["time"]
+		return self.state
+
+	def outputFnc(self):
+		record = self.plan.iloc[self.state["index"]]
+		# TODO: fix location check
+		speed = record["distance"] / (record["time_end"] - record["time"])
+		req = Request(round(self.state["time"] + self.timeAdvance(), 6), record["Sleepboot"], record["distance"], record["tugging"], record["Locatie van"], record["Locatie naar"], speed)
+		return {
+			self.reqs: req
+		}

collect.py

@@ -0,0 +1,30 @@
+import requests
+import json
+import pandas as pd
+
+resp = requests.get(
+	"https://webapps.portofantwerp.com/arcgis/rest/services/extern/detailkaart_fs/FeatureServer/60/query?where=LPL!='kip'&returnGeometry=true&outSR=4326&f=json&outFields=LPL,FILIAAL"
+)
+
+resp_d = json.loads(resp.text)
+
+berths = {
+	"lpl_id": [],
+	"filiaal": [],
+	"polygon": [],
+	"center_lon": [],
+	"center_lat": [],
+}
+
+for x in resp_d["features"]:
+	berths["lpl_id"].append(x["attributes"]["LPL"])
+	berths["filiaal"].append(x["attributes"]["FILIAAL"])
+	ring = x["geometry"]["rings"][0]
+	berths["polygon"].append(ring)
+	berths["center_lon"].append(sum([a for a, _ in ring]) / len(ring))
+	berths["center_lat"].append(sum([a for _, a in ring]) / len(ring))
+
+berths_pdf = pd.DataFrame(berths)
+print(berths_pdf)
+
+berths_pdf.to_csv("berths.csv")

de2/elements.py

@@ -0,0 +1,216 @@
+from pypdevs.DEVS import AtomicDEVS, CoupledDEVS
+from pypdevs.infinity import INFINITY
+import pandas as pd
+
+from dataclasses import dataclass
+
+@dataclass
+class Vessel:
+	mmsi: 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
+
+	usage_time: float = 0.0
+	idle_time: float = 0.0
+
+
+class Pool(AtomicDEVS):
+	def __init__(self, name):
+		super(Pool, self).__init__(name)
+
+		self.state = {
+			"waiting": {},
+			"should_exit": [],
+			"time": 0.0
+		}
+
+		self.req_in = self.addInPort("req_in")
+		self.vessel_in = self.addInPort("vessel_in")
+		self.vessel_out = self.addOutPort("vessel_out")
+
+	def extTransition(self, inputs):
+		self.state["time"] += self.elapsed
+		if self.req_in in inputs:
+			request = inputs[self.req_in]
+			# TODO: prefill with all vessels of the simulation?
+			if request["mmsi"] in self.state["waiting"]:
+				vessel = self.state["waiting"][request["mmsi"]]
+			else:
+				print("VESSEL %s DOES NOT EXIST IN POOL" % str(request["mmsi"]))
+				vessel = Vessel(request["mmsi"])
+
+			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]
+		if self.vessel_in in inputs:
+			vessel = inputs[self.vessel_in]
+			if vessel.mmsi in self.state["waiting"]:
+				raise ValueError("Cannot create duplicate vessels (loc = %s; mmsi = %s; d = %f)" %
+				                 (str(vessel.location), str(vessel.mmsi), vessel.total_distance))
+			else:
+				vessel.source = vessel.target
+				vessel.target = None
+				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"]
+
+				self.state["waiting"][vessel.mmsi] = vessel
+		return self.state
+
+	def timeAdvance(self):
+		if len(self.state["should_exit"]) > 0:
+			return 0.0
+		return INFINITY
+
+	def outputFnc(self):
+		if len(self.state["should_exit"]) > 0:
+			return {
+				self.vessel_out: self.state["should_exit"][0]
+			}
+		return {}
+
+	def intTransition(self):
+		self.state["time"] += self.timeAdvance()
+		if len(self.state["should_exit"]) > 0:
+			self.state["should_exit"].pop(0)
+		return self.state
+
+
+class Sailer(AtomicDEVS):
+	def __init__(self, name):
+		super(Sailer, self).__init__(name)
+
+		self.state = {
+			"vessels": [],
+			"time": 0.0
+		}
+
+		self.vessel_in = self.addInPort("vessel_in")
+		self.vessel_out = self.addOutPort("vessel_out")
+
+	def extTransition(self, inputs):
+		self.state["time"] += self.elapsed
+		for vessel in self.state["vessels"]:
+			x = vessel.velocity * self.elapsed
+			vessel.distance_left = max(0.0, vessel.distance_left - x)
+			vessel.time_until_departure = round(vessel.time_until_departure - self.elapsed, 6)
+		if self.vessel_in in inputs:
+			vessel = inputs[self.vessel_in]
+			self.state["vessels"].append(vessel)
+		self.state["vessels"].sort(key=lambda v: v.time_until_departure)
+		return self.state
+
+	def timeAdvance(self):
+		if len(self.state["vessels"]) > 0:
+			v = self.state["vessels"][0]
+			return v.time_until_departure
+		return INFINITY
+
+	def outputFnc(self):
+		if len(self.state["vessels"]) > 0:
+			return {
+				self.vessel_out: self.state["vessels"][0]
+			}
+		return {}
+
+	def intTransition(self):
+		elapsed = self.timeAdvance()
+		self.state["time"] += elapsed
+		self.state["vessels"].pop(0)
+		for vessel in self.state["vessels"]:
+			x = vessel.velocity * elapsed
+			vessel.distance_left = max(0.0, vessel.distance_left - x)
+			vessel.time_until_departure = round(max(0.0, vessel.time_until_departure - elapsed), 6)
+		if len(self.state["vessels"]) > 0:
+			self.state["vessels"].sort(key=lambda v: v.time_until_departure)
+		return self.state
+
+
+class Tracer(AtomicDEVS):
+	def __init__(self, name, ivef):
+		super(Tracer, self).__init__(name)
+
+		self.ivef = pd.read_csv(ivef)
+		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.state = {
+			"index": 0,
+			"time": 0.0
+		}
+
+		self.reqs = self.addOutPort("reqs")
+
+	def timeAdvance(self):
+		if self.state["index"] < len(self.ivef):
+			return round(self.ivef.iloc[self.state["index"]]["start"] / 1000 - self.state["time"], 6)
+		return INFINITY
+
+	def intTransition(self):
+		self.state["time"] += self.timeAdvance()
+		self.state["index"] += 1
+		return self.state
+
+	def outputFnc(self):
+		if self.state["index"] < len(self.ivef):
+			return {
+				self.reqs: self.ivef.iloc[self.state["index"]]
+			}
+		return {}
+
+
+class Clock(AtomicDEVS):
+	def __init__(self, name, interval=1):
+		super(Clock, self).__init__(name)
+		self.interval = interval
+		self.state = {
+			"time": 0.0
+		}
+
+	def timeAdvance(self):
+		return self.interval
+
+	def intTransition(self):
+		self.state["time"] += self.timeAdvance()
+		return self.state
+
+
+class Port(CoupledDEVS):
+	def __init__(self, name, ivef):
+		super(Port, self).__init__(name)
+
+		self.clock = self.addSubModel(Clock("clock"))
+		self.tracer = self.addSubModel(Tracer("tracer", ivef))
+		self.pool = self.addSubModel(Pool("pool"))
+		self.sailer = self.addSubModel(Sailer("sailer"))
+
+		self.connectPorts(self.tracer.reqs, self.pool.req_in)
+		self.connectPorts(self.pool.vessel_out, self.sailer.vessel_in)
+		self.connectPorts(self.sailer.vessel_out, self.pool.vessel_in)
+
+	def select(self, imm_children):
+		for child in imm_children:
+			if isinstance(child, Sailer):
+				return child
+		for child in imm_children:
+			if isinstance(child, Pool):
+				return child
+		return imm_children[0]

de2/plotter.py

@@ -0,0 +1,32 @@
+import pandas as pd
+import matplotlib.pyplot as plt
+
+
+ivef = pd.read_csv("../results-de/IVEF.csv")
+
+ivef["delta"] = (ivef["end"] - ivef["start"]) / 1000
+ivef["velocity"] = ivef["distance"] / ivef["delta"]
+
+sailing = ivef[ivef["location"].isnull()]
+mooring = ivef[~ivef["location"].isnull()]
+
+buckets1 = {}
+for vel in mooring["velocity"].to_numpy():
+	if vel > 0:
+		buckets1.setdefault(vel, 0)
+		buckets1[vel] += 1
+
+buckets2 = {}
+for vel in sailing["velocity"].to_numpy():
+	if vel > 0:
+		buckets2.setdefault(vel, 0)
+		buckets2[vel] += 1
+
+plt.xlabel("velocity (m/s)")
+plt.ylabel("number of vessels")
+plt.bar(buckets1.keys(), buckets1.values(), 0.01, label="in quay (mooring)")
+plt.bar(buckets2.keys(), buckets2.values(), 0.01, label="on canal/river (sailing)")
+plt.legend()
+plt.show()
+
+# velocity is most likely always 0 and 5 m/s = 18 km/h = 9.71 knots

de2/routing.py

@@ -0,0 +1,124 @@
+from os.path import exists
+import numpy as np
+import igraph as ig
+import json
+from geoplotlib.utils import haversine
+
+
+def euler_distance(x1, y1, x2, y2):
+    """
+    Computes the Euler distance between two points (x1, y1) and (x2, y2).
+    """
+    return np.sqrt((x1 - x2) ** 2 + (y2 - y1) ** 2)
+
+
+def BD72_to_WGS84(E, N):
+    """
+    Converts a coordinate from Belgian Lambert 72 (BD72) to GPS (WGS 84).
+    This basically applies a transformation to the right coordinates, around Belgium.
+
+    :note: BD72 is said to be not 100% accurate. While it is incredibly close, better use the
+    WGS84 datasets!
+
+    :seealso: https://epsg.io/31370, https://epsg.io/3395
+    :seealso: https://gis.stackexchange.com/questions/259121/transformation-functions-for-epsg3395-projection-vs-epsg3857
+    :seealso: https://pubs.usgs.gov/pp/1395/report.pdf, pp. 44-45, 63-64
+
+    This one was used for the computations!
+    :seealso: https://www.linz.govt.nz/guidance/geodetic-system/understanding-coordinate-conversions/projection-conversions/lambert-conformal-conic-geographic-transformation-formulae
+    """
+    a = 6378388
+    f = 1./297
+    lat1 = np.deg2rad(51.1666672333333)
+    lat2 = np.deg2rad(49.8333339)
+    lat0 = np.deg2rad(90)
+    lon0 = np.deg2rad(4.36748666666667)
+    N0 = 5400088.438
+    E0 = 150000.013
+    eps = 1e-10
+
+    # Constants
+    e = np.sqrt(2*f - f**2)
+    mx = lambda latx: np.cos(latx) / np.sqrt(1 - e**2 * np.sin(latx)**2)
+    m1 = mx(lat1)
+    m2 = mx(lat2)
+    tx = lambda latx: np.tan(np.pi/4 - latx/2) / ((1 - e * np.sin(latx))/(1 + e * np.sin(latx)))**(e/2)
+    t0 = tx(lat0)
+    t1 = tx(lat1)
+    t2 = tx(lat2)
+    n = (np.log(m1) - np.log(m2)) / (np.log(t1) - np.log(t2))
+    F = m1 / (n * t1**n)
+    rho0 = a * F * t0**n
+
+    # Computation
+    N -= N0
+    E -= E0
+    rho = np.sqrt(E**2 + (rho0 - N)**2)
+    assert n * rho > 0  # check that rho has the same sign as n
+    t = np.power(rho / (a * F), 1./n)
+    gamma = np.arctan2(E, rho0 - N)
+
+    lat = np.pi/2 - 2 * np.arctan(t)
+    while True:
+        lat2 = np.pi/2 - 2 * np.arctan(t * ((1 - e * np.sin(lat))/(1 + e * np.sin(lat)))**(e/2))
+        if np.abs(lat2 - lat) < eps:
+            lat = lat2
+            break
+        lat = lat2
+
+    lon = gamma / n + lon0
+
+    return np.rad2deg(lat), np.rad2deg(lon)
+
+
+def get_graph():
+    fname = "paths/cache.pickle.gz"
+    if exists(fname):
+        graph = ig.Graph.Read_Picklez(fname)
+        return graph
+    graph = ig.Graph()
+    with open("paths/assenWGS84.json", 'r') as file:
+        features = json.load(file)["features"]
+    for feature in features:
+        for path in feature["geometry"]["paths"]:
+            traj = []
+            for coord in path:
+                vc, d = get_closest_vertex(graph, coord[0], coord[1])
+                if d > 3.:
+                    v = graph.add_vertex()
+                    v["x"] = coord[0]
+                    v["y"] = coord[1]
+                    traj.append(v)
+                else:
+                    traj.append(vc)
+            if len(traj) == 1:
+                print("boo!", feature)
+                continue
+            for i in range(len(traj) - 1):
+                if traj[i].index == traj[i+1].index: continue
+                e = graph.add_edge(traj[i], traj[i+1])
+                e["distance"] = haversine(traj[i]["x"], traj[i]["y"], traj[i+1]["x"], traj[i+1]["y"])
+    return graph
+
+
+def get_closest_vertex(graph, x, y):
+    distance = float('inf')
+    out = None
+    for vertex in graph.vs:
+        d = haversine(vertex["x"], vertex["y"], x, y)
+        if d < distance:
+            distance = d
+            out = vertex
+    return out, distance
+
+
+def get_exact_vertex(graph, x, y):
+    for vertex in graph.vs:
+        if vertex["x"] == x and vertex["y"] == y:
+            return vertex
+    return None
+
+
+if __name__ == '__main__':
+    graph = get_graph()
+    graph.write_picklez("paths/cache.pickle.gz")

events.py

@@ -0,0 +1,22 @@
+from dataclasses import dataclass
+
+@dataclass
+class Request:
+	time: float
+	ship: str
+	distance: float
+	tugging: bool
+	source: str
+	destination: str
+	speed: float = 0.0
+
+@dataclass
+class Ship:
+	id: str
+	position: str
+	v: float
+	x: float
+	time_start: float
+	time_until_departure: float
+	time_used: float = 0
+	total_distance: float = 0

mapper.py

@@ -0,0 +1,228 @@
+import random
+
+from geoplotlib.layers import BaseLayer
+from geoplotlib.core import BatchPainter
+import geoplotlib
+from geoplotlib.utils import BoundingBox, epoch_to_str, haversine
+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
+
+
+class PoABLayer(BaseLayer):
+    """
+    Renders the map.
+    """
+
+    __colors__ = {
+        "highlight": [255, 255, 0, 100],
+        "nodes": [150, 0, 250, 200],
+        "sailing": [0, 150, 250, 200],
+        "tugging": [150, 250, 0, 200]
+    }
+
+    def __init__(self, sim=None):
+        self.sim = sim
+        self.sim_active = False
+
+        self.berths = pd.read_csv("berths.csv")
+        self.berths["used"] = False
+        self.paths = pd.read_csv("2022.csv")
+        self.paths.drop_duplicates(["Locatie van", "Locatie naar"], inplace=True)
+
+        for _, traj in self.paths.iterrows():
+            self.berths.loc[self.berths["lpl_id"] == traj["Locatie van"], "used"] = True
+            self.berths.loc[self.berths["lpl_id"] == traj["Locatie naar"], "used"] = True
+        self.berths = self.berths[self.berths["used"]]
+
+        self.paths = []
+        self.painter = BatchPainter()
+
+        # 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.
+
+        Args:
+            x (numeric):    The x-position to draw this shape at.
+            y (numeric):    The y-position to draw this shape at.
+        """
+        self.painter.points(x, y, 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)
+
+
+        # 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:
+            model = self.sim.model
+            time = model.clock.state["time"]
+            ts = int(time + model.tracer.starting_time / 1000)
+            ui_manager.info("TIME: %s" % epoch_to_str(ts))
+
+            for vessel in model.sailer.state["vessels"]:
+                if vessel.source is not None or vessel.target is not None:
+                    # sx, sy = proj.lonlat_to_screen(np.asarray([vessel.source[0]]), np.asarray([vessel.source[1]]))
+                    # tx, ty = proj.lonlat_to_screen(np.asarray([vessel.target[0]]), np.asarray([vessel.target[1]]))
+                    if not isinstance(vessel.location, str) and math.isnan(vessel.location):
+                        self.painter.set_color(self.__colors__["sailing"])
+                    else:
+                        self.painter.set_color(self.__colors__["tugging"])
+                    # self.painter.lines(sx, sy, tx, ty, width=2)
+
+                    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]))
+                    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)
+
+                    tot_distance = vessel.total_distance
+                    traveled_approx = tot_distance - vessel.distance_left
+                    delta = time - model.sailer.state["time"]
+                    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 - 1]
+                            p2 = path[dx]
+                            p = sum(dists[dx:]) / tot_dist
+                            if 1 - perc - p == 0:
+                                percentage = 1
+                            else:
+                                percentage = (percentage - perc) / (1 - perc - p)
+                            break
+                        perc += s / tot_dist
+                    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]]))
+
+                    dx = (tx - sx) * percentage
+                    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("%.3f m/s" % vessel.velocity)
+                    self.draw_shape(sx + dx, sy + dy)
+
+        ui_manager.tooltip("\n".join(tooltips))
+        self.painter.batch_draw()
+
+    def bbox(self):
+        box = [max(self.berths["center_lat"]), max(self.berths["center_lon"]), min(self.berths["center_lat"]), min(self.berths["center_lon"])]
+        return BoundingBox(north=box[0], east=box[1], south=box[2], west=box[3])
+
+    def on_key_release(self, key, modifiers):
+        if key == pyglet.window.key.SPACE:
+            if self.sim is not None:
+                self.sim.simulate()
+                self.sim_active = True
+
+
+if __name__ == '__main__':
+    from de2.elements import Port
+    from pypdevs.simulator import Simulator
+
+    port = Port("PoAB", "results-de/IVEF.csv")
+
+    sim = Simulator(port)
+    sim.setClassicDEVS()
+    sim.setRealTime(scale=0.0002)
+    sim.setTerminationTime(1_209_600)
+
+    layer = PoABLayer(sim)
+
+    geoplotlib.set_window_size(640, 760)
+    geoplotlib.tiles_provider({
+        'url': lambda zoom, xtile, ytile: 'http://a.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.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))

parquet/check_delays.py

@@ -0,0 +1,12 @@
+import pandas as pd
+
+tug = pd.read_csv("../results/205351190.csv")
+
+prev_time = tug.iloc[0]["ts"]
+
+for rix, row in tug.iterrows():
+	ts = row["ts"]
+	if ts - prev_time > 900_000_000:
+		print(ts)
+		print("Large gap of %f ms at row %d" % (ts - prev_time, rix))
+	prev_time = ts

parquet/ct2de.py

@@ -0,0 +1,85 @@
+import pandas as pd
+import numpy as np
+from dataclasses import dataclass
+from geoplotlib.utils import haversine
+import matplotlib.path as mplPath
+from os.path import exists
+
+docks = pd.read_csv("../berths.csv")
+polygons = {k["lpl_id"]: eval(k["polygon"]) for _, k in docks.iterrows()}
+
+@dataclass
+class Vessel:
+	mmsi: str
+	start: int
+	end: int
+	location: str
+	latlon_s: tuple
+	latlon_t: tuple
+	distance: float
+
+def is_inside(lon, lat, polygon):
+	ppath = mplPath.Path(np.array(polygon))
+	return ppath.contains_point((lon, lat))
+
+
+def get_dock(lat, lon):
+	# loop through all dock indexes and do insideness checks
+	for pid, poly in polygons.items():
+		if is_inside(lon, lat, poly):
+			return pid
+	return ""
+
+def analyze(fname):
+	if not exists(fname):
+		print(fname, "does not exist!")
+		return None
+	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"])
+
+	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)
+
+	if exists("results-de/%s-de.csv" % str(vessel.mmsi)):
+		print("Already done %s" % str(vessel.mmsi))
+		return None
+	# time = first["ts"]
+	# dts = []
+	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"]
+		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"]
+			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]
+	# print("AVG dt:", sum(dts) / len(dts), max(dts), min(dts))
+	res.to_csv("../results-de/%s-de.csv" % str(vessel.mmsi))
+
+if __name__ == '__main__':
+	# import sys
+	# analyze(sys.argv[2])
+	# analyze("results/205254890.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))
+		print(idx, "done")

parquet/parquet-plotter.py

@@ -0,0 +1,70 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+from os.path import exists
+import scipy.stats as stats
+
+TFILE = "../tugs.csv"
+tugs = pd.read_csv(TFILE)
+
+velocities = pd.Series([])
+df = pd.DataFrame()
+
+mmsi = tugs["MMSI"].dropna().astype(np.int32)
+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)
+	df = pd.concat([df, tug], ignore_index=True)
+	d = tug["distance"]  # meters
+	t = (tug["end"] - tug["start"]) / 1000  # seconds
+	v = d / t
+	v.replace([np.inf, -np.inf], np.nan, inplace=True)
+	v.dropna(inplace=True)
+	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)
+
+raise 4
+
+fig, ax1 = plt.subplots()
+ax1.set_xlabel("velocity (m/s)")
+ax2 = ax1.twinx()
+ax1.set_ylabel("amount")
+# ax2.set_ylabel("", color="tab:red")
+ax2.tick_params(axis='y', labelcolor='r')
+
+def plot(series):
+	buckets = {}
+
+	for v in series:
+		vr = round(v, 6)
+		buckets.setdefault(vr, 0)
+		buckets[vr] += 1
+
+	if 0 in buckets:
+		del buckets[0]
+
+	ax1.bar(buckets.keys(), buckets.values(), 0.1)
+	ax1.set_ylim((0, 1.1 * max(buckets.values())))
+
+	mu = series.mean()
+	std = series.std()
+	mx = int(np.ceil(series.max()))
+
+	print(mu, std, mx)
+
+	beta = mu / std**2
+	alpha = mu * beta
+
+	x = np.linspace(0, mx, mx * 2)
+	# y1 = stats.gamma.pdf(x, a=alpha, scale=1./beta)
+	y1 = stats.norm.pdf(x, mu, std)
+	ax2.plot(x, y1, c='r', label=r"$\Gamma(%.3f, %.3f)$" % (mu, std))
+	ax2.set_ylim((0, 1.1 * y1.max()))
+
+plot(velocities)
+plt.show()

parquet/parquet-splitter.py

@@ -0,0 +1,26 @@
+import pandas as pd
+import numpy as np
+
+PFILE = "part-00000-tid-257698335488947142-c2a078b1-adfb-43b6-98b6-947dfaff5eda-5415-1.c000.parquet"
+pq = pd.read_parquet(PFILE)
+
+print("Read Parquet of size", len(pq))
+
+TFILE = "tugs.csv"
+tugs = pd.read_csv(TFILE)
+
+print("Read Tugboat Information")
+
+mmsi = tugs["MMSI"].dropna().astype(np.int32)
+S = 0
+for idx in mmsi:
+	# print(idx)
+	tug = pq[pq["mmsi"] == idx]
+	S += len(tug)
+	if len(tug) > 0:
+		tug.to_csv("results/%s.csv" % str(idx), columns=["mmsi", "ts", "lon", "lat", "speed"], index=False)
+		# tug.to_parquet("results/%s.parquet" % str(idx), compression=None)
+print(S, "events stored")
+
+
+

plotter.py

@@ -0,0 +1,117 @@
+import pandas as pd
+import numpy as np
+import scipy.stats as stats
+import random
+import matplotlib.pyplot as plt
+
+# wind = pd.read_csv("wind.csv", parse_dates=[1, 4], date_format="%Y-%m-%d %H:%M:%S.%f")
+# wind.dropna()
+# start = wind.iloc[0]["Aanvang tijd"]
+# wind["start"] = wind["Aanvang tijd"] - start
+# wind["start"] = wind["start"].apply(lambda x: x.to_numpy() / np.timedelta64(1, 'h'))
+# wind["end"] = wind["Einde tijd"] - start
+# wind["end"] = wind["end"].apply(lambda x: x.to_numpy() / np.timedelta64(1, 'h'))
+#
+# ax1 = plt.subplot(211)
+# plt.tick_params('x', labelbottom=False)
+# ax1.set_ylabel("wind force (B)")
+#
+# ax2 = plt.subplot(212, sharex=ax1)
+# ax2.set_ylabel("wind direction")
+# ax2.set_xlabel("time (h)")
+#
+# for sluis in wind["Sluis"].unique():
+# 	info = wind[wind["Sluis"] == sluis]
+# 	xs, ys, ws = [], [], []
+# 	for _, row in info.iterrows():
+# 		xs.append(row["start"])
+# 		xs.append(row["end"])
+# 		ys.append(row["Aanvang windkracht"])
+# 		ys.append(row["Einde windkracht"])
+# 		ws.append(row["Aanvang windrichting"])
+# 		ws.append(row["Einde windrichting"])
+# 		if len(ws) > 2 and (ws[-1] == np.nan or ws[-2] == np.nan):
+# 			print(row["start"])
+# 	ax1.plot(xs, ys, label=sluis)
+# 	ax2.plot(xs, ws, label=sluis)
+#
+# plt.tight_layout()
+# plt.legend()
+# plt.show()
+#
+# exit()
+
+df = pd.read_csv("2022.csv", parse_dates=[1, 3, 4, 5], date_format="%Y-%m-%d %H:%M:%S")
+
+df["duur"] = df["Einde tijd"] - df["Aanvang tijd"]
+df["duur"] = df["duur"].apply(lambda x: x.to_numpy() / np.timedelta64(1, 'h'))
+df["Afstand"] = df["Afstand"].astype(float)
+df["velocity"] = df["Afstand"] / df["duur"]   # km/h
+
+df.replace([np.inf, -np.inf], np.nan, inplace=True)
+df = df.dropna()
+
+task = df["Taak"] == "Slepen"
+vel = df["velocity"][task]
+
+# remove outliers
+Q1 = vel.quantile(0.25)
+Q3 = vel.quantile(0.75)
+IQR = Q3 - Q1
+trueList = ~((vel < (Q1 - 1.5 * IQR)) | (vel > (Q3 + 1.5 * IQR)))
+print("full:", len(vel))
+inv = vel[~trueList]
+vel = vel[trueList]
+print("IQR:", len(vel), Q1 - 1.5 * IQR, Q3 + 1.5 * IQR)
+print("inv:", len(inv), inv.min(), inv.max())
+
+# unique color for each tug
+colors = { i: "#"+''.join([random.choice('0123456789ABCDEF') for _ in range(6)]) for i in df["Sleepboot"].unique() }
+df["color"] = [colors[b] for b in df["Sleepboot"]]
+
+# plt.scatter(df["time"][task][trueList], vel, s=1, c=df["color"][task][trueList])
+# plt.title("Sailing (Grouped by Start Time)")
+# plt.xlabel("time (h)")
+# plt.ylabel("velocity (km/h)")
+# plt.show()
+
+
+fig, ax1 = plt.subplots()
+ax1.set_xlabel("velocity (km/h)")
+ax2 = ax1.twinx()
+ax1.set_ylabel("amount")
+# ax2.set_ylabel("", color="tab:red")
+ax2.tick_params(axis='y', labelcolor='r')
+
+def plot(series, title):
+	plt.title(title)
+
+	buckets = {}
+
+	for v in series:
+		buckets.setdefault(v, 0)
+		buckets[v] += 1
+
+	if 0 in buckets:
+		del buckets[0]
+
+	ax1.bar(buckets.keys(), buckets.values(), label=title)
+	ax1.set_ylim((0, 1.1 * max(buckets.values())))
+
+	mu = series.mean()
+	std = series.std()
+	mx = int(np.ceil(series.max()))
+
+	beta = mu / std**2
+	alpha = mu * beta
+
+	x = np.linspace(0, mx, mx * 2)
+	# y1 = stats.gamma.pdf(x, a=alpha, scale=1./beta)
+	y1 = stats.norm.pdf(x, mu, std)
+	ax2.plot(x, y1, c='r', label=r"$\Gamma(%.3f, %.3f)$" % (mu, std))
+	ax2.set_ylim((0, 1.1 * y1.max()))
+
+plot(vel, "Sailing")
+
+plt.legend()
+plt.show()

tugs.csv

@@ -0,0 +1,52 @@
+NAME,max pull strength backwards,max pull strength forward,TUG_BASE_STAT,TUG_IN_USE,TUG_OUT_OF_ORDER,SHC_IMO_NUMBER,MMSI
+10,50,50,602,03/04/2009,NULL,NULL,205351190
+11,40,50,602,03/04/2009,NULL,NULL,205351290
+20,40,53,602,03/04/2009,NULL,NULL,205360090
+21 - METHATUG,53,53,NULL,03/04/2009,NULL,NULL,205360190
+22,40,53,1622,03/04/2009,NULL,NULL,205252690
+30,55,55,602,16/01/2012,NULL,NULL,205501790
+31,55,55,602,03/04/2009,NULL,NULL,205504190
+32,55,55,1622,03/04/2009,NULL,NULL,205506390
+40,70,70,602,10/07/2012,NULL,9602095,205625000
+41,70,70,602,19/09/2012,NULL,9602100,205626000
+42,70,70,602,14/01/2013,NULL,9602112,205627000
+43,70,70,NULL,20/02/2013,NULL,9602124,205628000
+51,72,75,602,18/03/2021,NULL,9890123,205094000
+52,72,75,1622,18/06/2021,NULL,9890135,205100000
+70,17,23,NULL,01/05/2012,12/06/2013,NULL,205251490
+73,17,23,NULL,03/05/2012,27/09/2012,NULL,205000000
+82,19,24,NULL,NULL,NULL,7521546,205264290
+83,19,24,NULL,01/04/2009,22/04/2022,7505152,205264890
+84,19,24,NULL,01/04/2009,17/10/2016,NULL,205254890
+85,19,24,NULL,06/04/2009,22/04/2022,7505164,205254790
+86,19,24,NULL,03/04/2009,24/01/2015,NULL,205257290
+90,24,32,1622,03/04/2009,NULL,NULL,205273990
+91,24,32,602,03/04/2009,NULL,NULL,205337590
+92,24,32,602,03/04/2009,NULL,NULL,205203390
+50,72,75,602,31/12/2019,NULL,9842932,205045000
+BUGSIER 21,54,54,NULL,04/03/2021,NULL,9214989,
+FAIRPLAY I,70,70,NULL,01/01/2017,NULL,9365128,305542000
+FAIRPLAY III,70,70,NULL,13/07/2012,NULL,9365116,244616000
+FAIRPLAY TAMADABA,74,74,NULL,01/05/2021,NULL,9817250,305498000
+FAIRPLAY XIV,70,70,NULL,01/05/2012,NULL,9541708,246517000
+FAIRPLAY-23,52,52,NULL,01/05/2012,31/12/2020,9148776,
+MULTRATUG 1,39,39,NULL,01/01/2017,NULL,9035008,244870324
+MULTRATUG 16,36,36,NULL,01/05/2012,NULL,8414154,209126000
+MULTRATUG 19,85,85,NULL,01/07/2012,NULL,9630250,245875000
+MULTRATUG 27,63,63,NULL,01/01/2017,NULL,9667875,244660818
+MULTRATUG 28,72,72,NULL,01/05/2021,NULL,9727962,229835000
+MULTRATUG 29,85,85,NULL,31/05/2021,NULL,9795816,244850093
+MULTRATUG 31,83,83,NULL,01/05/2021,NULL,9695614,244870671
+MULTRATUG 32,70,70,NULL,31/03/2018,04/06/2018,9796042,
+MULTRATUG 33,70,70,NULL,05/06/2018,05/07/2018,9796054,
+MULTRATUG 4,75,75,NULL,31/05/2012,NULL,9360582,256381000
+MULTRATUG 5,65,65,NULL,01/05/2012,31/12/2020,9350161,
+MULTRATUG 7,62,62,NULL,11/05/2021,NULL,9817808,244150102
+VB SEINE,65,65,NULL,24/10/2020,NULL,9476393,205766000
+UNION AMBER,65,65,NULL,01/05/2012,NULL,9365130,205474000
+UNION EAGLE,85,85,NULL,01/05/2012,NULL,9406441,205300000
+UNION GRIZZLY,65,65,NULL,01/01/2014,NULL,9397121,205483000
+UNION HAWK,85,85,NULL,01/05/2012,NULL,9406439,205234000
+UNION JADE,65,65,NULL,01/05/2012,NULL,9365142,205475000
+UNION KODIAK,65,65,NULL,01/05/2012,NULL,9397119,205484000
+UNION PANDA,65,65,NULL,01/05/2012,NULL,9502697,205545000