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

@@ -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))