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

Do not always write all properties

Add option `--debug` to write all properties to output files. The
properties which are not always written are not needed for further
processing of the graph. To implement the option, we define variables
for properties, so we revert in part commit 2e3c6606.
parent c2a77999
No related branches found
No related tags found
No related merge requests found
......@@ -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"
......
......@@ -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] \
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment