From bdd0ce416cdc21d251262f678b159663b9992baf Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Fri, 29 Apr 2022 23:33:39 +0200
Subject: [PATCH] Rewrite the loop on vertices

Motivation: performance. We eliminate repeated reading of the same
information from the shapefiles. If the length of the segment is
smaller than 14 days then we re-use for the last part of the segment
some of the properties from the first part. Also, we do not compute
properties for the first part of the segment if the segment has no
predecessor and for the last part of the segment if the segment has no
successor.
---
 cost_function.py | 52 +++++++++++++++++++++++++++++-------------------
 1 file changed, 31 insertions(+), 21 deletions(-)

diff --git a/cost_function.py b/cost_function.py
index 2598bfbc..5b5bc26b 100755
--- a/cost_function.py
+++ b/cost_function.py
@@ -144,39 +144,49 @@ n_days_avg = 7 # number of days to average
 print("Iterating on vertices...")
 
 for n in g.vertices():
-    if n.in_degree() != 0 or n.out_degree() != 0:
-        len_seg = len(g.vp.inst_eddies[n])
+    if n.in_degree() != 0:
+        # Define properties for beginning of the segment:
+        properties = node_to_prop(g.vp.inst_eddies[n][:n_days_avg], e_overestim,
+                                  array_d_init, handlers)
+        g.vp.first_av_rad[n], g.vp.first_av_ros[n] \
+            = calculate_radii_rossby(properties)
+        g.vp.pos_first[n] = properties[0]["pos"] # in degrees
 
-        if len_seg > n_days_avg:
-            # The segment is longer than the number of days over which
-            # to average
+    if n.out_degree() != 0:
+        # Define properties for end of the segment:
 
-            # First 7 days calculation
-            properties = node_to_prop(g.vp.inst_eddies[n][:n_days_avg],
-                                      e_overestim, array_d_init, handlers)
-            g.vp.first_av_rad[n], g.vp.first_av_ros[n] \
-                = calculate_radii_rossby(properties)
-            g.vp.pos_first[n] = properties[0]["pos"] # in degrees
+        len_seg = len(g.vp.inst_eddies[n])
+
+        if n.in_degree() == 0 or len_seg > n_days_avg:
+            # We have to read more from the shapefiles and redefine
+            # properties.
+
+            if n.in_degree() == 0 or len_seg > 2 * n_days_avg:
+                # We cannot use part of properties from the beginning
+                # of the segment.
+                properties = node_to_prop(g.vp.inst_eddies[n][- n_days_avg:],
+                                          e_overestim, array_d_init, handlers)
+            else:
+                # assertion: n.in_degree() != 0 and n_days_avg <
+                # len_seg < 2 * n_days_avg
+
+                # We can use part of the properties from the beginning
+                # of the segment.
+                properties = properties[len_seg - n_days_avg:] \
+                    + node_to_prop(g.vp.inst_eddies[n][n_days_avg:],
+                                   e_overestim, array_d_init, handlers)
 
-            # Last 7 days calculation:
-            properties = node_to_prop(g.vp.inst_eddies[n][- n_days_avg:],
-                                      e_overestim, array_d_init, handlers)
             g.vp.last_av_rad[n], g.vp.last_av_ros[n] \
                 = calculate_radii_rossby(properties)
-            g.vp.pos_last[n] = properties[- 1]["pos"] # in degrees
         else:
             # The number of eddies in the segment is lower than or
             # equal to the number of days over which to average. The
             # values for the end of the segment will be the same as
             # for the begining, except for the position.
-            properties = node_to_prop(g.vp.inst_eddies[n], e_overestim,
-                                      array_d_init, handlers)
-            g.vp.first_av_rad[n], g.vp.first_av_ros[n] \
-                = calculate_radii_rossby(properties)
             g.vp.last_av_rad[n] = g.vp.first_av_rad[n]
             g.vp.last_av_ros[n] = g.vp.first_av_ros[n]
-            g.vp.pos_first[n] = properties[0]["pos"] # in degrees
-            g.vp.pos_last[n] = properties[- 1]["pos"] # in degrees
+
+        g.vp.pos_last[n] = properties[- 1]["pos"] # in degrees
 
 t1 = time.perf_counter()
 timings.write(f"iterating on vertices: {t1 - t0:.0f} s\n")
-- 
GitLab