From 738166b798620e8e2ffb0ad323a7ebed77fe51c6 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Wed, 20 Apr 2022 13:25:34 +0200
Subject: [PATCH] Polish
---
cost_function.py | 103 +++++++++++++++++++++--------------------------
segments.py | 4 +-
2 files changed, 47 insertions(+), 60 deletions(-)
diff --git a/cost_function.py b/cost_function.py
index 8c27ce7e..d92408f5 100755
--- a/cost_function.py
+++ b/cost_function.py
@@ -1,8 +1,7 @@
#!/usr/bin/env python3
-"""A script that takes a segmented graph in the gt format and performs
-the non-local cost function calculation where each edge will have a
-cost function assigned to it.
+"""A script that takes the graph of segments without cost functions
+and computes the non-local cost functions attributed to edges.
Input:
@@ -33,21 +32,21 @@ import bisect
def calculate_radii_and_rossby(start, end, inc, segment, e_overestim, handlers,
array_d_init):
-
+
radii = 0 #[m]
rossby = 0 #[1/s]
-
+
days_modifier = 0
-
+
for i in range(start, end, inc):
current_eddy = report_graph.node_to_date_eddy(segment[i], e_overestim)
i_SHPC = get_SHPC(array_d_init, current_eddy['date_index'])
-
+
# calculate the location in the shapefile
location = util_eddies.comp_ishape(handlers[i_SHPC],
current_eddy['date_index'],
current_eddy['eddy_index'])
-
+
# now that we have the location in the shapefiles, we need to
# get the radius and the rossby number
shapeRec = handlers[i_SHPC]["readers"]["extremum"].shapeRecord(location)
@@ -55,23 +54,23 @@ def calculate_radii_and_rossby(start, end, inc, segment, e_overestim, handlers,
lat_in_deg = shapeRec.shape.points[0][1]
#[deg]
f = 2*2*math.pi/(24*3600)*math.sin(math.radians(lat_in_deg)) # [1/s]
-
+
V_max = shapeRec.record[4] #[m/s]
R_Vmax = handlers[i_SHPC]["readers"]["max_speed_contour"]\
.record(location)['r_eq_area'] * 1000 #[m]
-
+
if (V_max < 100):
# calculate Ro and Delta_Ro
Ro = V_max / (f * R_Vmax) #[]
else:
Ro = 0
days_modifier += 1
-
+
####### RADII #######
radii += R_Vmax # [m]
####### ROSSBY ######
rossby += Ro # []
-
+
return {"radii": radii, "rossby": rossby, "days_modifier": days_modifier}
def get_SHPC(array_d_ini, date_index):
@@ -97,13 +96,13 @@ delta_cent_std = 8.0388
delta_ro_mean = -0.0025965 # []
delta_ro_std = 5.2168
-delta_r_mean = -0.0094709 * 1000 #[m]
+delta_r_mean = -0.0094709 * 1000 #[m]
delta_r_std = 8.6953 * 1000
# Load the graph_tool file:
t0 = time.perf_counter()
timings = open("timings.txt", "w")
-print('Loading gt file...')
+print('Loading graph...')
g = graph_tool.Graph()
try:
@@ -112,12 +111,12 @@ except FileNotFoundError:
g.load('segments.graphml')
t1 = time.perf_counter()
+print('Loading done...')
print("Input graph:")
print("Number of vertices:", g.num_vertices())
print("Number of edges:", g.num_edges())
print("Internal properties:")
g.list_properties()
-print('Loading done...')
timings.write(f"loading: {t1 - t0} s\n")
timings.close()
@@ -140,19 +139,17 @@ array_d_init = [handler["d_init"] for handler in handlers]
# change if there is a change over the number of days to average
num_of_days_to_avg = 7
-# iterate on the vertices
+print("Iterating on vertices...")
+
for n in g.vertices():
- # Get the segment and the number of days
segment = g.vp.segment[n]
+ num_of_days = len(segment)
+
+ # Calculate the indexes and dates:
- # calculate the indexes and dates
first = report_graph.node_to_date_eddy(segment[0], e_overestim)
first_SHPC = get_SHPC(array_d_init, first['date_index'])
- num_of_days = len(segment)
-
- # start processing
-
last = report_graph.node_to_date_eddy(segment[-1], e_overestim)
last_SHPC = get_SHPC(array_d_init, last['date_index'])
@@ -164,20 +161,19 @@ for n in g.vertices():
last['date_index'],
last['eddy_index'])
- # grab the centers
-
+ # grab the centers:
first_pos = handlers[first_SHPC]["readers"]["extremum"]\
.shape(first_loc).points[0]
last_pos = handlers[last_SHPC]["readers"]["extremum"]\
.shape(last_loc).points[0]
- ##### STORE POSITIONS IN THE VPS ######
+ # Store positions in the vertex properties:
g.vp.pos_first[n] = first_pos # [deg, deg]
g.vp.pos_last[n] = last_pos # [deg, deg]
- # if the segments are longer than the # of days over which to avg
-
if (num_of_days > num_of_days_to_avg):
+ # The segment is longer than the number of days over which to average
+
first_radii = 0 # [m]
last_radii = 0 # [m]
@@ -200,22 +196,17 @@ for n in g.vertices():
# Average and assign the rossbies:
first_rossby = first_res['rossby'] / (num_of_days_to_avg - modifier)
g.vp.first_av_ros[n] = first_rossby
- else:
- # there is division by zero, average rossby is undefinied
- pass
# Last 7 days calculation
- last_res = calculate_radii_and_rossby(len(segment) - 1,
- len(segment) - (num_of_days_to_avg + 1),
- -1,
- segment, e_overestim,
- handlers, array_d_init)
+ last_res = calculate_radii_and_rossby(len(segment) - 1, len(segment) -
+ (num_of_days_to_avg + 1), -1,
+ segment, e_overestim, handlers,
+ array_d_init)
# Average and assign the last radii
last_radii = last_res['radii'] / num_of_days_to_avg
g.vp.last_av_rad[n] = last_radii
-
# grab the days modifier
modifier = last_res['days_modifier']
@@ -223,13 +214,10 @@ for n in g.vertices():
# Average and assign the rossbies:
last_rossby = last_res['rossby'] / (num_of_days_to_avg - modifier)
g.vp.last_av_ros[n] = last_rossby
- else:
- # there is division by zero, average rossby is undefinied
- pass
- # else, the number of eddies in a segment is lower than the # of
- # days over which to average, the values will be the same except
- # for the positions
else:
+ # The number of eddies in a segment is lower than the number
+ # of days over which to average. The values will be the same
+ # except for the positions.
res = calculate_radii_and_rossby(0, num_of_days, 1, segment,
e_overestim, handlers, array_d_init)
@@ -241,19 +229,13 @@ for n in g.vertices():
rossby = res['rossby'] / (num_of_days - modifier)
g.vp.first_av_ros[n] = rossby
g.vp.last_av_ros[n] = rossby
- else:
- # there is division by zero, average rossby is undefinied
- pass
# Average and assign the radii
-
radii = res['radii'] / num_of_days
g.vp.first_av_rad[n] = radii
g.vp.last_av_rad[n] = radii
-###############################
-# Calculate the cost function #
-###############################
+print("Iterating on edges...")
for edge in g.edges():
source_node = edge.source()
@@ -262,25 +244,28 @@ for edge in g.edges():
cf = -10000
lat_for_conv = (g.vp.pos_last[source_node][1] +
- g.vp.pos_first[target_node][1]) / 2 # latitude needed for conversion of degrees to kilometers
+ g.vp.pos_first[target_node][1]) / 2
+ # (latitude needed for conversion of degrees to kilometers)
lat_for_conv = math.radians(lat_for_conv) # need to convert to radians
- # 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])
+ # 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])
if (lon_diff > 300):
lon_diff = 360 - lon_diff
- # calculate Delta_cent: numbers used for conversion obtained from:
- # https://stackoverflow.com/questions/1253499/simple-calculations-for-working-with-lat-lon-and-km-distance
- Delta_Cent = math.sqrt(( lon_diff * 111.32 * math.cos(lat_for_conv) )**2 +
- ( (g.vp.pos_last[source_node][1] - g.vp.pos_first[target_node][1]) * 110.574 )**2)
+ Delta_Cent = math.sqrt((lon_diff * 111.32 * math.cos(lat_for_conv))**2
+ + ((g.vp.pos_last[source_node][1]
+ - g.vp.pos_first[target_node][1]) * 110.574)**2)
# calculate the first term
first_term = ((Delta_Cent - delta_cent_mean)/delta_cent_std) ** 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]
+ Delta_Ro = g.vp.last_av_ros[source_node] \
+ - g.vp.first_av_ros[target_node]
else:
print("At least one of the rossbies is invalid.")
#Delta_Ro = delta_ro_mean
@@ -291,7 +276,8 @@ for edge in g.edges():
# 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 = g.vp.last_av_rad[source_node] \
+ - g.vp.first_av_rad[target_node]
# Calculate the third term
third_term = ((Delta_R_Vmax - delta_r_mean)/delta_r_std) ** 2
@@ -306,5 +292,6 @@ for edge in g.edges():
g.ep.nl_cost_function[edge] = cf
+print("Saving...")
g.save('segments_cost_functions.gt')
print('All done')
diff --git a/segments.py b/segments.py
index 89a696d6..67c6ed1e 100755
--- a/segments.py
+++ b/segments.py
@@ -72,8 +72,8 @@ for v in g.vertices():
t1 = time.perf_counter()
print(f'Done collapsing in {t1 - t0:.0f} s')
t0 = t1
-print('Empty nodes:', len(verts_to_del))
-print('Deleting empty nodes...')
+print('Number of circumvented nodes:', len(verts_to_del))
+print('Deleting circumvented nodes...')
g.remove_vertex(verts_to_del, fast=True)
t1 = time.perf_counter()
print(f"Done deleting in {t1 - t0:.0f} s")
--
GitLab