diff --git a/gt_segmented_cf.py b/gt_segmented_cf.py
new file mode 100644
index 0000000000000000000000000000000000000000..fd3ddc85fea1a75e9926b8ec34562146db093adc
--- /dev/null
+++ b/gt_segmented_cf.py
@@ -0,0 +1,383 @@
+#!/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 each will have a cost
+# function assigned to it.
+# Inputs (to be changed within the script):
+ # Orientation
+ # node_id_param_file location and name
+ # gt graph name and location
+ # shapefiles names and location
+# Output: a gt graph with cost functions and an xml graph with cost functions
+
+import graph_tool
+import time
+import json
+import math
+from os import path
+import shapefile
+import datetime
+from numpy import loadtxt
+
+######################
+# Utility functions: #
+######################
+
+def node_to_date_eddy(node_index, e_overestim, d_init):
+ """Convert from node identification in graph to eddy identification in
+ shapefiles.
+
+ """
+ k = (node_index - 1) // e_overestim
+ d = d_init + k
+ return {"eddy_index": node_index - k * e_overestim, "date_index": k,
+ "days_1950": d}
+
+def comp_ishape(d_init, ishape_last, days_1950, eddy_index):
+ k = days_1950 - d_init
+
+ if k == 0:
+ ishape = eddy_index - 1
+ else:
+ ishape = ishape_last[k - 1] + eddy_index
+
+ return ishape
+
+#######################
+# Necessary functions #
+#######################
+
+def calculate_radii_and_rossby(start, end, inc, segment, e_overestim, d_init, ext_shp_rcd, max_rcd, ishape_last):
+
+ radii = 0 #[m]
+ rossby = 0 #[1/s]
+
+ days_modifier = 0
+
+ for i in range(start, end, inc):
+ current_eddy = node_to_date_eddy(segment[i], e_overestim, d_init)
+ date = datetime.date(1950,1,1) + datetime.timedelta(current_eddy['days_1950'])
+ year = date.year
+
+ # calculate the location in the shapefile
+ d_init_1 = ext_shp_rcd[year][0].record[1]
+ location = comp_ishape(d_init_1, ishape_last[year],
+ current_eddy['days_1950'], current_eddy['eddy_index'])
+
+ # now that we have the location in the shapefiles, we need to get the radius and the rossby number
+
+ # rossby:
+ lat_in_deg = ext_shp_rcd[year][location].shape.points[0][1] #[deg]
+ f = 2*2*math.pi/(24*3600)*math.sin(math.radians(lat_in_deg)) # [1/s]
+
+ V_max = ext_shp_rcd[year][location].record[4] #[m/s]
+ R_Vmax = max_rcd[year][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}
+
+###############################
+# grab e_overestim and d_init #
+###############################
+
+node_id_param_file = "/data/misic/Convert_Global/node_id_param.json"
+with open(node_id_param_file) as f: node_id_param = json.load(f)
+# assign attributes to the whole graphs
+e_overestim = node_id_param["e_overestim"]
+d_init = node_id_param["d_init"]
+
+################################################
+# set some values needed for the cost function #
+################################################
+
+delta_cent_mean = 3.8481 # [km]
+delta_cent_std = 8.0388
+
+delta_ro_mean = -0.0025965 # []
+delta_ro_std = 5.2168
+
+delta_r_mean = -0.0094709 * 1000 #[m]
+delta_r_std = 8.6953 * 1000
+
+orientation = 'anti'
+
+##########################
+# load the edgelist file #
+##########################
+
+t0 = time.perf_counter()
+
+print('Loading gt file...')
+g = graph_tool.Graph()
+g.load(f'global_segmented_{orientation}.gt')
+t1 = time.perf_counter()
+print(f'Loading done: {g}')
+print(t1 - t0)
+
+g.set_fast_edge_removal()
+
+##################################
+# Calculating the Cost Function: #
+##################################
+
+#######################
+# Load the shapefiles #
+#######################
+
+# load the shapefiles and the ishape_last files
+
+print('Loading extrema and ishape_last files.')
+
+# setup the dicts for the shapefiles and the ishape_last files
+extremum = {}
+max_speed = {}
+ishape_last = {}
+year_bool = {}
+
+shp_tr_dir = '/data/guez/Oceanic_eddies/Global_Matlab'
+
+for year in range(1993,2019):
+ shp_dir = path.join(shp_tr_dir, f'SHPC_{orientation}_{year}')
+ extremum[year] = shapefile.Reader(path.join(shp_dir, 'extremum.shp'))
+ max_speed[year] = shapefile.Reader(path.join(shp_dir, 'max_speed_contour.shp'))
+ ishape_last[year] = loadtxt(path.join(shp_dir, 'ishape_last.txt'),
+ dtype = int, delimiter = "\n", unpack = False)
+ year_bool[year] = False
+
+print('Loading of the shapefiles and ishape_last files done.')
+
+count = 1
+
+t0 = time.perf_counter()
+# change if there is a change over the number of days to average
+num_of_days_to_avg = 7
+
+# we iterate on the years because the segments are non sequential in time
+for yr in range(1993, 2019):
+
+ print(f'Current year: {yr}')
+
+ # grab the shapefiles
+
+ ext_shp_rcd = {}
+ max_rcd = {}
+
+ ext_shp_rcd[yr] = extremum[yr].shapeRecords()
+ max_rcd[yr] = max_speed[yr].records()
+
+ # check if we are reaching the last two years or the last year
+ if (yr <= 2017):
+ ext_shp_rcd[yr+1] = extremum[yr + 1].shapeRecords()
+ max_rcd[yr+1] = max_speed[yr + 1].records()
+
+ if (yr <= 2016):
+ ext_shp_rcd[yr+2] = extremum[yr + 2].shapeRecords()
+ max_rcd[yr+2] = max_speed[yr + 2].records()
+
+ # iterate on the vertices
+ for n in g.vertices():
+ # Get the segment and the number of days
+ segment = g.vp.segment[n]
+
+ # calculate the indexes and dates
+ first = node_to_date_eddy(segment[0], e_overestim, d_init)
+
+ # get days since 1950
+ first_days = first['days_1950']
+
+ # get year
+ first_date = datetime.date(1950,1,1) + datetime.timedelta(first['days_1950'])
+ first_year = first_date.year
+
+ #print(f'Segment size: {len(segment)}, first: {first_year}, last: {last_year}, yr: {yr}')
+
+ # chech if we are in the current year
+
+ if (first_year == yr):
+ num_of_days = len(segment)
+
+ # start processing
+
+ first_pos = []
+ last_pos = []
+
+ last = node_to_date_eddy(segment[-1], e_overestim, d_init)
+ last_days = last['days_1950']
+
+ last_date = datetime.date(1950,1,1) + datetime.timedelta(last['days_1950'])
+ last_year = last_date.year
+
+ print(f'Segment size: {len(segment)}, first: {first_year}, last: {last_year}, yr: {yr}')
+
+ # calculate the location in the shapefile
+
+ d_init_1 = ext_shp_rcd[first_year][0].record[1]
+ d_init_2 = ext_shp_rcd[last_year][0].record[1]
+
+ first_loc = comp_ishape(d_init_1, ishape_last[first_year], first['days_1950'], first['eddy_index'])
+ last_loc = comp_ishape(d_init_2, ishape_last[last_year], last['days_1950'], last['eddy_index'])
+
+ # grab the centers
+
+ first_pos = ext_shp_rcd[first_year][first_loc].shape.points[0]
+ last_pos = ext_shp_rcd[last_year][last_loc].shape.points[0]
+
+ ##### STORE POSITIONS IN THE VPS ######
+ 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):
+ first_radii = 0 # [m]
+ last_radii = 0 # [m]
+
+ first_rossby = 0 # []
+ last_rossby = 0 # []
+
+ # First 7 days calculation
+ first_res = calculate_radii_and_rossby(0, num_of_days_to_avg, 1, segment, e_overestim,
+ d_init, ext_shp_rcd, max_rcd, ishape_last)
+
+ # average and assign radii
+ first_radii = first_res['radii'] / num_of_days_to_avg
+ g.vp.first_av_rad[n] = first_radii
+
+ # grab the days modifier
+ modifier = first_res['days_modifier']
+
+ if (num_of_days_to_avg - modifier > 0):
+ # 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, d_init, ext_shp_rcd,
+ max_rcd, ishape_last)
+
+ # 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']
+
+ if (num_of_days_to_avg - modifier > 0):
+ # 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:
+ res = calculate_radii_and_rossby(0, num_of_days, 1, segment, e_overestim,
+ d_init, ext_shp_rcd, max_rcd, ishape_last)
+
+ # grab the days modifier
+ modifier = res['days_modifier']
+
+ if (num_of_days - modifier > 0):
+ # Average and assign the rossbies:
+ 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 #
+###############################
+
+for edge in g.edges():
+ source_node = edge.source()
+ target_node = edge.target()
+
+ 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
+ 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])
+ 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)
+
+ # 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]
+ else:
+ print("At least one of the rossbies is invalid.")
+ #Delta_Ro = delta_ro_mean
+ Delta_Ro = 0
+
+ # Calculate the second term
+ second_term = ((Delta_Ro - delta_ro_mean)/delta_ro_std ) ** 2
+
+ # 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]
+
+ # Calculate the third term
+ third_term = ((Delta_R_Vmax - delta_r_mean)/delta_r_std) ** 2
+
+ #############################
+ # calculate the cost function
+ #############################
+
+ cf = math.sqrt(first_term + second_term + third_term)
+
+ # assign as weight to the edge
+ g.ep.nl_cost_function[edge] = cf
+
+
+################################
+# Writing to an GT Graph File #
+################################
+
+#g.save("segments_cyclo_gl_cf.xml")
+#print("Done writing xml.")
+g.save(f'segmented_{orientation}_cf.gt')
+print("Done writing gt.")
+g.save(f'segmented_{orientation}_cf.xml')
+print("Done writing xml.")
+#g.save(f'{orientation}_segmented_cf.gv', fmt = "dot")
+#print("Done writing dot file.")
+
+print('All done')