Average over a maximum number of dates

The old behaviour: averaging over a maximum number of instantaneous
eddies, becomes an option.

We do not use bisect because `max_delta` should be small, around 6,
and most often there should be only one time step between
instantaneous eddies, so the first guess at distance `max_delta`
should be close to the result. Moreover, we would need to define a
date function if we used bisect.

Trajectories/cost_function.py

@@ -103,12 +103,78 @@ def node_to_prop(node_list, e_overestim, SHPC, orientation):
     return properties
 
 
-def search_beg(inst_eddies, max_delta):
-    return min(max_delta + 1, len(inst_eddies))
+def search_beg(inst_eddies, max_delta, avg_fix, e_overestim):
+    """Return an index ip in inst_eddies. inst_eddies is a list of
+    inter-date indices of instantaneous eddies. len(inst_eddies) must
+    be >= 1. max_delta must be >= 0. If avg_fix then inst_eddies[:ip]
+    is a fixed number of instantaneous eddies, if the list is long
+    enought. If not avg_fix then ip = bisect.bisect_right(inst_eddies,
+    date(inst_eddies[0]) + max_delta, key = date), although the
+    computation is not implemented that way.
 
+    """
+    ip = min(max_delta + 1, len(inst_eddies))
+
+    if not avg_fix and ip >= 2:
+        d_max = (
+            util_eddies.node_to_date_eddy(
+                inst_eddies[0], e_overestim, only_date=True
+            )
+            + max_delta
+        )
+        # {date(elem) > d_max for elem in inst_eddies[ip:]}
+
+        while (
+            ip >= 2
+            and util_eddies.node_to_date_eddy(
+                inst_eddies[ip - 1], e_overestim, only_date=True
+            )
+            > d_max
+        ):
+            ip -= 1
+
+        # {date(elem) <= d_max for elem in inst_eddies[:ip] and
+        # date(elem) > d_max for elem in inst_eddies[ip:]}
+
+    # {1 <= ip <= min(max_delta + 1, len(inst_eddies))}
+    return ip
+
+
+def search_end(inst_eddies, max_delta, avg_fix, e_overestim):
+    """Return an index ip in inst_eddies. inst_eddies is a list of
+    inter-date indices of instantaneous eddies. len(inst_eddies) must
+    be >= 1. max_delta must be >= 0. If avg_fix then inst_eddies[ip:]
+    is a fixed number of instantaneous eddies, if the list is long
+    enought. If not avg_fix then ip = bisect.bisect_left(inst_eddies,
+    date(inst_eddies[-1]) - max_delta, key = date), although the
+    computation is not implemented that way.
+
+    """
+    ip = max(len(inst_eddies) - max_delta - 1, 0)
+
+    if not avg_fix and ip <= len(inst_eddies) - 2:
+        d_min = (
+            util_eddies.node_to_date_eddy(
+                inst_eddies[-1], e_overestim, only_date=True
+            )
+            - max_delta
+        )
+        # {date(elem) < d_min for elem in inst_eddies[:ip]}
+
+        while (
+            ip <= len(inst_eddies) - 2
+            and util_eddies.node_to_date_eddy(
+                inst_eddies[ip], e_overestim, only_date=True
+            )
+            < d_min
+        ):
+            ip += 1
 
-def search_end(inst_eddies, max_delta):
-    return max(len(inst_eddies) - max_delta - 1, 0)
+        # {date(elem) < d_min for elem in inst_eddies[:ip] and
+        # date(elem) >= d_min for elem in inst_eddies[ip:]}
+
+    # {max(len(inst_eddies) - max_delta - 1, 0) <= ip <= len(inst_eddies) - 1}
+    return ip
 
 
 t0 = time.perf_counter()
@@ -129,6 +195,11 @@ parser.add_argument(
 parser.add_argument(
     "--debug", help="save properties to output file", action="store_true"
 )
+parser.add_argument(
+    "--avg_fix",
+    help="average over a fixed number of instantaneous eddies",
+    action="store_true",
+)
 args = parser.parse_args()
 
 # Set some values needed for the cost function:
@@ -185,7 +256,9 @@ print("Iterating on vertices...")
 for n in g.vertices():
     if n.in_degree() != 0:
         # Define properties for beginning of the segment:
-        ip_beg = search_beg(g.vp.inst_eddies[n], max_delta)
+        ip_beg = search_beg(
+            g.vp.inst_eddies[n], max_delta, args.avg_fix, g.gp.e_overestim
+        )
         properties = node_to_prop(
             g.vp.inst_eddies[n][:ip_beg],
             g.gp.e_overestim,
@@ -203,7 +276,9 @@ for n in g.vertices():
         if ip_beg < len(g.vp.inst_eddies[n]):
             # We have to read more from the shapefiles and redefine
             # properties.
-            ip_end = search_end(g.vp.inst_eddies[n], max_delta)
+            ip_end = search_end(
+                g.vp.inst_eddies[n], max_delta, args.avg_fix, g.gp.e_overestim
+            )
 
             if ip_beg <= ip_end:
                 # We cannot use part of properties from the beginning