From 480800c9a468acad7e52eb9d92cbb9facdfd5cf2 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 14 Sep 2023 20:56:27 +0200
Subject: [PATCH] Rewrite using new algorithm
The old algorithm did not give correct results on tests
`Greece_trajectories`, `Test_order_edges` and `Traj_component`. We
abandon the idea of processing edges in date order. Instead, we now
process vertices in topological order.
---
Trajectories/trajectories.py | 173 ++++++++++++-----------------------
1 file changed, 61 insertions(+), 112 deletions(-)
diff --git a/Trajectories/trajectories.py b/Trajectories/trajectories.py
index 7742cf5c..b734f6c7 100755
--- a/Trajectories/trajectories.py
+++ b/Trajectories/trajectories.py
@@ -13,81 +13,48 @@ instantaneous eddies in each trajectory are in ascending order.
"""
import graph_tool
+from graph_tool import topology
import json
import sys
-import numpy as np
-import report_graph
-
-def assign_traj_edge(source, target, traj_prop, ind_traj, g, traj_segm,
- expanded_traj):
- """Assign trajectories to the source and target of a given edge. This
- means updating ind_traj, traj_prop, traj_segm and expanded_traj.
-
- """
-
- if traj_prop[source] == - 1 and traj_prop[target] == - 1:
- ind_traj += 1
- traj_prop[source] = ind_traj
- traj_prop[target] = ind_traj
- segments = [int(g.vertex_properties.name[source]),
- int(g.vertex_properties.name[target])]
- traj_segm.append(segments)
- inst_eddies = list(g.vertex_properties.inst_eddies[source]) \
- + list(g.vertex_properties.inst_eddies[target])
- expanded_traj.append(inst_eddies)
- elif traj_prop[source] == - 1:
- # assert traj_prop[target] != - 1
- # target is already in a trajectory. Can we extend the
- # trajectory of target backward?
- in_neighbors = g.get_in_neighbors(target, [traj_prop])
-
- if np.any(in_neighbors[:, 1] == traj_prop[target]):
- # target already has a direct predecessor in its
- # trajectory so we create a new trajectory for source.
- ind_traj += 1
- traj_prop[source] = ind_traj
- segments = [int(g.vertex_properties.name[source])]
- traj_segm.append(segments)
- inst_eddies = list(g.vertex_properties.inst_eddies[source])
- expanded_traj.append(inst_eddies)
- else:
- traj_prop[source] = traj_prop[target]
- segment = int(g.vertex_properties.name[source])
- traj_segm[traj_prop[source]].insert(0, segment)
- inst_eddies = list(g.vertex_properties.inst_eddies[source])
- expanded_traj[traj_prop[source]] = inst_eddies \
- + expanded_traj[traj_prop[source]]
- elif traj_prop[target] == - 1:
- # assert traj_prop[source] != - 1
- # source is already in a trajectory. Can we extend the
- # trajectory of source forward?
- out_neighbors = g.get_out_neighbors(source, [traj_prop])
-
- if np.any(out_neighbors[:, 1] == traj_prop[source]):
- # source already has a direct successor in its trajectory
- # so we create a new trajectory for target.
- ind_traj += 1
- traj_prop[target] = ind_traj
- segments = [int(g.vertex_properties.name[target])]
- traj_segm.append(segments)
- inst_eddies = list(g.vertex_properties.inst_eddies[target])
- expanded_traj.append(inst_eddies)
- else:
- traj_prop[target] = traj_prop[source]
- segment = int(g.vertex_properties.name[target])
- traj_segm[traj_prop[target]].append(segment)
- inst_eddies = list(g.vertex_properties.inst_eddies[target])
- expanded_traj[traj_prop[target]].extend(inst_eddies)
+def new_traj(ind_traj, traj_prop, n, g, traj_segm, expanded_traj):
+ """Assign new trajectory to node index or vertex descriptor n."""
+
+ ind_traj += 1
+ traj_prop[n] = ind_traj
+ segments = [int(g.vertex_properties.name[n])]
+ traj_segm.append(segments)
+ inst_eddies = list(g.vertex_properties.inst_eddies[n])
+ expanded_traj.append(inst_eddies)
return ind_traj
if len(sys.argv) != 2: sys.exit("Required argument: input-graph")
g = graph_tool.load_graph(sys.argv[1])
-traj_prop = g.new_vertex_property('int', val = - 1)
-# Trajectory index to which a segment belongs. - 1 means the segment
-# does not belong yet to any trajectory. The trajectory index is also
-# the index in lists traj_segm and expanded_traj.
+if "cost_function" in g.edge_properties:
+ closest_succ = g.new_vertex_property('int')
+ # Index of the closest direct successor of a node: tail of the
+ # out-edge with smallest cost. - 1 means there is no successor.
+
+ # Set closest_succ for all nodes:
+ for n in g.iter_vertices():
+ all_cost = g.get_out_edges(n,
+ eprops = [g.edge_properties.cost_function])
+ # numpy array with dtype float64
+
+ if all_cost.size == 0:
+ # n has no successor
+ closest_succ[n] = - 1
+ else:
+ i_min = all_cost[:, 2].argmin()
+ closest_succ[n] = all_cost[i_min, 1]
+else:
+ # This happens if there are only isolated segments.
+ print("No cost values in the input file.")
+
+traj_prop = g.new_vertex_property('int')
+# Trajectory index to which a segment belongs. The trajectory index is
+# also the index in lists traj_segm and expanded_traj.
ind_traj = - 1 # initialize index of trajectory
traj_segm = [] # an element of traj_segm is a trajectory, represented
@@ -95,53 +62,35 @@ traj_segm = [] # an element of traj_segm is a trajectory, represented
expanded_traj = [] # an element of expanded_traj is a trajectory,
# represented by a list of instantaneous eddies
-if "cost_function" in g.edge_properties:
- # Build my_subgraphs, a dictionary of nodes at each date:
-
- my_subgraphs = {}
-
- for v in g.vertices():
- date_index = report_graph.node_to_date_eddy(int(g.vp.name[v]),
- g.gp.e_overestim,
- only_date = True)
-
- if date_index in my_subgraphs:
- my_subgraphs[date_index].append(int(v))
+for n in topology.topological_sort(g):
+ if g.vertex(n).in_degree() >= 1:
+ # Find the index of the closest direct predecessor of a node:
+ # head of the in-edge with smallest cost.
+
+ all_cost = g.get_in_edges(n, eprops = [g.edge_properties.cost_function])
+ # numpy array with dtype float64
+
+ i_min = all_cost[:, 2].argmin()
+ closest_pred = all_cost[i_min, 0]
+
+ if closest_succ[closest_pred] == n:
+ # Assign to n the trajectory of closest_pred. This means
+ # updating ind_traj, traj_prop, traj_segm and
+ # expanded_traj. closest_pred is already in a trajectory,
+ # traj_prop[closest_pred]] is defined. We extend the
+ # trajectory of closest_pred forward.
+ traj_prop[n] = traj_prop[closest_pred]
+ segment = int(g.vertex_properties.name[n])
+ traj_segm[traj_prop[n]].append(segment)
+ inst_eddies = list(g.vertex_properties.inst_eddies[n])
+ expanded_traj[traj_prop[n]].extend(inst_eddies)
else:
- my_subgraphs[date_index] = [int(v)]
-
- for date_index in sorted(my_subgraphs):
- # Consider all outgoing edges at a given date:
-
- edges = []
- all_cost = []
-
- for n in my_subgraphs[date_index]:
- for source, target, cost \
- in g.iter_out_edges(n, [g.edge_properties.cost_function]):
- edges.append([n, target])
- all_cost.append(cost)
-
- # Array of indices that sort on cost:
- ind_cost = np.argsort(all_cost)
-
- # Loop over edges in ascending order of cost:
- for i in ind_cost:
- ind_traj = assign_traj_edge(*edges[i], traj_prop, ind_traj, g,
- traj_segm, expanded_traj)
-else:
- # This happens if there are only isolated segments.
- print("No cost values in the input file.")
+ ind_traj = new_traj(ind_traj, traj_prop, n, g, traj_segm,
+ expanded_traj)
+ else:
+ ind_traj = new_traj(ind_traj, traj_prop, n, g, traj_segm, expanded_traj)
-# Create trajectories for isolated segments:
-for v in g.vertices():
- if traj_prop[v] == - 1:
- ind_traj += 1
- traj_prop[v] = ind_traj
- segments = [int(g.vertex_properties.name[v])]
- traj_segm.append(segments)
- inst_eddies = list(g.vertex_properties.inst_eddies[v])
- expanded_traj.append(inst_eddies)
+ # traj_prop[n] is defined
print("Number of trajectories: ", len(traj_segm))
--
GitLab