Skip to content
Snippets Groups Projects
Commit 738166b7 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Polish

parent b79d15fd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
cost_function.py
+ 45
− 58
View file @ 738166b7
#!/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')
segments.py
+ 2
− 2
View file @ 738166b7
......@@ -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")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment