From 77ec1d5fe1659911cd492b18a9c67888c6675e7b Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Wed, 29 May 2024 09:03:40 +0200
Subject: [PATCH] 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 | 87 ++++++++++++++++++++++++++++++++---
1 file changed, 81 insertions(+), 6 deletions(-)
diff --git a/Trajectories/cost_function.py b/Trajectories/cost_function.py
index b9bfa1ab..b795e0dd 100755
--- a/Trajectories/cost_function.py
+++ b/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
--
GitLab