diff --git a/Trajectories/Tests/tests.json b/Trajectories/Tests/tests.json
index ddbe3827b27884d314743a1e949ed7be4c34f0ac..f091ba8dcb9e94878903eaadbaf6efee96866cdc 100644
--- a/Trajectories/Tests/tests.json
+++ b/Trajectories/Tests/tests.json
@@ -43,7 +43,7 @@
"$tests_old_dir/Extraction_eddies_region_2_noise/SHPC",
"Anticyclones",
"$tests_old_dir/Segments/segments.graphml",
- "segments_cost_functions.graphml"
+ "segments_cost_functions.graphml", "--debug"
],
"exclude_cmp": [
"timings.txt"
@@ -64,7 +64,7 @@
"$tests_old_dir/Greece/SHPC",
"Anticyclones",
"$tests_old_dir/Greece_segments/segments.graphml",
- "segments_cost_functions.graphml"
+ "segments_cost_functions.graphml", "--debug"
],
"exclude_cmp": [
"timings.txt"
diff --git a/Trajectories/cost_function.py b/Trajectories/cost_function.py
index 12fe16d7ec0a0bed7eeafe8bb8527b70c2f66921..17a3696480b122b8f0bf7b68e550dadd3ab31336 100755
--- a/Trajectories/cost_function.py
+++ b/Trajectories/cost_function.py
@@ -92,6 +92,8 @@ parser.add_argument("input_segments", help = "input graph of segments without "
"cost, suffix .gt (graph-tool) or .graphml")
parser.add_argument("output_segments", help = "output graph of segments with "
"cost, suffix .gt (graph-tool) or .graphml")
+parser.add_argument("--debug", help = "save properties to output file",
+ action = "store_true")
args = parser.parse_args()
# Set some values needed for the cost function:
@@ -116,12 +118,21 @@ t1 = time.perf_counter()
timings.write(f"loading: {t1 - t0:.0f} s\n")
t0 = t1
-g.vp['pos_first'] = g.new_vp('vector<double>')
-g.vp['pos_last'] = g.new_vp('vector<double>')
-g.vp['first_av_rad'] = g.new_vp('float')
-g.vp['first_av_ros'] = g.new_vp('float')
-g.vp['last_av_rad'] = g.new_vp('float')
-g.vp['last_av_ros'] = g.new_vp('float')
+pos_first = g.new_vp('vector<double>')
+pos_last = g.new_vp('vector<double>')
+first_av_rad = g.new_vp('float')
+first_av_ros = g.new_vp('float')
+last_av_rad = g.new_vp('float')
+last_av_ros = g.new_vp('float')
+
+if args.debug:
+ g.vp['pos_first'] = pos_first
+ g.vp['pos_last'] = pos_last
+ g.vp['first_av_rad'] = first_av_rad
+ g.vp['first_av_ros'] = first_av_ros
+ g.vp['last_av_rad'] = last_av_rad
+ g.vp['last_av_ros'] = last_av_ros
+
g.ep['cost_function'] = g.new_ep('float')
SHPC = util_eddies.SHPC_class(args.SHPC_dir, args.orientation)
n_days_avg = 7 # number of days to average
@@ -132,9 +143,8 @@ for n in g.vertices():
# Define properties for beginning of the segment:
properties = node_to_prop(g.vp.inst_eddies[n][:n_days_avg],
g.gp.e_overestim, SHPC, args.orientation)
- 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
+ first_av_rad[n], first_av_ros[n] = calculate_radii_rossby(properties)
+ pos_first[n] = properties[0]["pos"] # in degrees
if n.out_degree() != 0:
# Define properties for end of the segment:
@@ -161,17 +171,16 @@ for n in g.vertices():
+ node_to_prop(g.vp.inst_eddies[n][n_days_avg:],
g.gp.e_overestim, SHPC, args.orientation)
- g.vp.last_av_rad[n], g.vp.last_av_ros[n] \
- = calculate_radii_rossby(properties)
+ last_av_rad[n], last_av_ros[n] = calculate_radii_rossby(properties)
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.
- g.vp.last_av_rad[n] = g.vp.first_av_rad[n]
- g.vp.last_av_ros[n] = g.vp.first_av_ros[n]
+ last_av_rad[n] = first_av_rad[n]
+ last_av_ros[n] = first_av_ros[n]
- g.vp.pos_last[n] = properties[- 1]["pos"] # in degrees
+ pos_last[n] = properties[- 1]["pos"] # in degrees
t1 = time.perf_counter()
timings.write(f"iterating on vertices: {t1 - t0:.0f} s\n")
@@ -182,31 +191,28 @@ for edge in g.edges():
source_node = edge.source()
target_node = edge.target()
- latitude = math.radians((g.vp.pos_last[source_node][1] +
- g.vp.pos_first[target_node][1]) / 2)
+ latitude = math.radians((pos_last[source_node][1] +
+ pos_first[target_node][1]) / 2)
# (latitude needed for conversion of degrees to kilometers)
# Because of the wrapping issue (360° wrapping incorrectly to 0°),
# we check for that here:
- lon_diff = abs(g.vp.pos_last[source_node][0]
- - g.vp.pos_first[target_node][0])
+ lon_diff = abs(pos_last[source_node][0] - pos_first[target_node][0])
if lon_diff > 300: lon_diff = 360 - lon_diff
Delta_Cent = math.sqrt((lon_diff * 111.32 * math.cos(latitude))**2
- + ((g.vp.pos_last[source_node][1]
- - g.vp.pos_first[target_node][1]) * 110.574)**2)
+ + ((pos_last[source_node][1]
+ - pos_first[target_node][1]) * 110.574)**2)
# Rossbies:
- if (g.vp.first_av_ros[target_node] and g.vp.last_av_ros[source_node]):
- Delta_Ro = g.vp.last_av_ros[source_node] \
- - g.vp.first_av_ros[target_node]
+ if (first_av_ros[target_node] and last_av_ros[source_node]):
+ Delta_Ro = last_av_ros[source_node] - first_av_ros[target_node]
else:
# At least one of the rossbies is invalid.
Delta_Ro = 0
# R_Vmax 1 and 2 already exist, just get the delta
- Delta_R_Vmax = g.vp.last_av_rad[source_node] \
- - g.vp.first_av_rad[target_node]
+ Delta_R_Vmax = last_av_rad[source_node] - first_av_rad[target_node]
# Calculate the cost and assign to the edge:
g.ep.cost_function[edge] \