ret_cost_test

from __future__ import division

import pickle, csv
import util
import numpy as np
import matplotlib.pyplot as plt
from math import log

mapbox_access_token = 'pk.eyJ1IjoiZ2plZSIsImEiOiJjangzY2F5MDcwMGlpNDhwbWtzbTJ6azBmIn0.P8vS2x_gtfBpWJwWgC3Sbw'

alpha = 48
beta = 78 * 8

def save_color_dict():

	color_dict = {}
	color_dict['age'] = '#0072B2'
	color_dict['oldest'] = '#0072B2'
	color_dict['fragility'] = '#56B4E9'
	color_dict['weakest'] = '#56B4E9'
	color_dict['traffic'] = '#E69F00'
	color_dict['busiest'] = '#E69F00'
	color_dict['composite'] = '#009E73'
	color_dict['OAT'] = '#D55E00'
	color_dict['OAT total'] = '#D55E00'
	color_dict['Sobol, exp. cost'] = '#CC79A7'
	color_dict['Sobol exp. total'] = '#CC79A7'
	color_dict['Sobol, perc'] = '#CC79A7'
	color_dict['Sobol'] = '#CC79A7'
	color_dict['p0.2'] = '#DDACE3'
	color_dict['p0.8'] = '#824E49'

	with open('color_dict.pkl','wb') as f:
		pickle.dump(color_dict,f)

def get_color_dict():

	with open('color_dict.pkl','rb') as f:
		color_dict = pickle.load(f)

	return color_dict

def get_color(series):

	color_dict = get_color_dict()

	return color_dict[series]

def get_master_dict():
	with open('input/20140114_master_bridge_dict.pkl',
			  'rb') as f:  # this version includes area for computation of repair cost and bridge-specific retrofit factor omega
		master_dict = pickle.load(f)  # has 1743 keys. One per highway bridge. (NOT BART)

	return master_dict

def get_sf_fullr_dict():
	with open('input/sf_fullr_dict.pkl', 'rb') as f:
		sf_dict = pickle.load(f)

	with open('input/sf_fullr_bridge_ids.pkl', 'rb') as f:
		bridge_ids = pickle.load(f)

	return sf_dict, bridge_ids

def load_undamaged_stats():
	# G = mahmodel.get_graph()
	#
	# demand = bd.build_demand('input/BATS2000_34SuperD_TripTableData.csv',
	# 						 'input/superdistricts_centroids_dummies.csv')
	#
	# no_damage_travel_time, no_damage_vmt, no_damage_trips_made, _ = cbs.precompute_network_performance()
	#
	# undamaged_stats = [no_damage_travel_time, no_damage_vmt, no_damage_trips_made]
	#
	# with open('undamaged_stats_sf_full.pkl', 'wb') as f:
	# 	pickle.dump(undamaged_stats, f)

	with open('undamaged_stats_sf_full.pkl', 'rb') as f:
		undamaged_stats = pickle.load(f)

	return undamaged_stats

def load_individual_undamaged_stats():
	undamaged_stats = load_undamaged_stats()

	return undamaged_stats[0], undamaged_stats[1], undamaged_stats[2]  # tt, vmt, trips made

def compute_weighted_average_performance(lnsas, map_weights, num_damage_maps, travel_times, vmts, trips_made,
										 no_damage_travel_time, no_damage_vmt, no_damage_trips_made, direct_costs):
	# Compute weighted average of performance metrics for a single sample of a fragility function vector.

	scenarios = len(lnsas) # number of scenarios under consideration

	# GB ADDITION -- computed weighted average (expectation) of travel time and other metrics of interest
	average_travel_time = 0
	average_trips_made = 0
	average_vmt = 0
	average_direct_costs = 0

	for j in range(0, len(lnsas)):  # for every scenario considered
		w = map_weights[j]
		temp_times = np.asarray(travel_times[np.arange(start=j, stop=scenarios*num_damage_maps, step=scenarios)])
		temp_trips = np.asarray(trips_made[np.arange(start=j, stop=scenarios*num_damage_maps, step=scenarios)])
		temp_vmts = np.asarray(vmts[np.arange(start=j, stop=scenarios*num_damage_maps, step=scenarios)])
		temp_direct_costs = np.asarray(direct_costs[np.arange(start=j, stop=scenarios*num_damage_maps, step=scenarios)])
		assert temp_trips.shape[0] == num_damage_maps, 'Error -- wrong number of trips.'
		assert temp_times.shape[0] == num_damage_maps, 'Error -- wrong number of times.'
		assert temp_vmts.shape[0] == num_damage_maps, 'Error -- wrong number of vmts.'
		average_travel_time += w *np.average(temp_times)
		average_trips_made += w *np.average(temp_trips)
		average_vmt += w*np.average(temp_vmts)
		average_direct_costs += w*np.average(temp_direct_costs)

	# add the scenario of no earthquake
	average_travel_time += (1 - sum(map_weights)) * no_damage_travel_time
	average_trips_made += (1 - sum(map_weights)) * no_damage_trips_made
	average_vmt += (1 - sum(map_weights)) * no_damage_vmt

	average_delay_cost = alpha*max(0,((average_travel_time - no_damage_travel_time) / 3600))  # travel times are in seconds, so convert to units of monetary units/hour*hours --> monetary units per day; assume travel times increase with damage

	average_connectivity_cost = beta*max(0, (no_damage_trips_made - average_trips_made))  # units of monetary units/hour*lost trips/day*hours/(trips*days)--> monetary units per day; assume total flows decrease with damage

	assert average_delay_cost >= 0, 'ERROR in compute_indirect_costs(): delay cost is negative.'
	assert average_connectivity_cost >= 0, 'ERROR in compute_indirect_costs(): connectivity cost is negative.'

	average_indirect_cost = average_delay_cost + average_connectivity_cost

	return average_travel_time, average_vmt, average_trips_made, average_direct_costs, average_delay_cost, average_connectivity_cost, average_indirect_cost

def compute_percent_reduction(baseline, new_value):

	return (new_value-baseline)/baseline*100


def get_baseline_retrofit_results(n_retrofits, training, detailed=False): # assuming S = 45 #TODO -- this is for retrofits with revised average

	no_damage_travel_time, no_damage_vmt, no_damage_trips_made = load_individual_undamaged_stats()

	filename = '_sf_fullr'

	if not training:
		output_folder = 'sobol_output/retrofits/ret_revised_avg/r' + str(n_retrofits) + '/'
		scenarios = 45
		map_indices_input = 'sobol_input/sf_fullr_testing_map_indices.pkl'  # S = 45 for training sf_fullr
		map_weights_input = 'sobol_input/sf_fullr_testing_map_weights.pkl'  # S = 45 for training sf_fullr
	else:
		output_folder = 'sobol_output/retrofits/ret_revised_avg_s30/r' + str(n_retrofits) + '/'
		scenarios = 30
		map_indices_input = 'sobol_input/sf_fullr_training_map_indices.pkl'  # S = 30 for training sf_fullr
		map_weights_input = 'sobol_input/sf_fullr_training_map_weights.pkl'  # S = 30 for training sf_fullr

	# store the results
	fX_times_output = output_folder + 'fX_times' + filename  # travel times for f_X
	fX_trips_output = output_folder + 'fX_trips' + filename  # trips made for f_X
	fX_vmts_output = output_folder + 'fX_vmts' + filename  # VMTs for f_X
	fX_avg_times_output = output_folder + 'fX_avg_time' + filename  # average TT
	fX_avg_trips_output = output_folder + 'fX_avg_trips' + filename  # average trips made
	fX_avg_vmts_output = output_folder + 'fX_avg_vmts' + filename  # average VMT
	fX_delay_costs_output = output_folder + 'fX_delay_costs' + filename
	fX_conn_costs_output = output_folder + 'fX_conn_costs' + filename
	fX_indirect_costs_output = output_folder + 'fX_indirect_costs' + filename
	fX_direct_costs_output = output_folder + 'fX_direct_costs' + filename
	fX_exp_indirect_cost_output = output_folder + 'fX_exp_indirect_costs' + filename
	fX_exp_direct_cost_output = output_folder + 'fX_exp_direct_costs' + filename
	fX_expected_cost_output = output_folder + 'fX_exp_costs' + filename

	damage_x_output = output_folder + 'damage_x' + filename

	# save data for f_X
	with open(damage_x_output, 'rb') as f:
		damage_tracker = pickle.load(f)

	with open(fX_times_output, 'rb') as f:  # save raw performance data
		f_X_times = pickle.load(f)
	with open(fX_trips_output, 'rb') as f:
		f_X_trips = pickle.load(f)
	with open(fX_vmts_output, 'rb') as f:
		f_X_vmts = pickle.load(f)

	with open(fX_avg_times_output, 'rb') as f:  # save average (expected) performance data
		f_X_avg_time = pickle.load(f)
	with open(fX_avg_trips_output, 'rb') as f:
		f_X_avg_trip = pickle.load(f)
	with open(fX_avg_vmts_output, 'rb') as f:
		f_X_avg_vmt = pickle.load(f)

	with open(fX_delay_costs_output, 'rb') as f:
		f_X_delay_costs = pickle.load(f)
	with open(fX_conn_costs_output, 'rb') as f:
		f_X_conn_costs = pickle.load(f)
	with open(fX_direct_costs_output, 'rb') as f:
		f_X_direct_costs = pickle.load(f)
	with open(fX_indirect_costs_output, 'rb') as f:
		f_X_indirect_costs = pickle.load(f)

	with open(fX_exp_direct_cost_output, 'rb') as f:
		f_X_exp_direct_cost = pickle.load(f)
	with open(fX_exp_indirect_cost_output, 'rb') as f:
		f_X_exp_indirect_cost = pickle.load(f)
	with open(fX_expected_cost_output, 'rb') as f:
		f_X_exp_cost = pickle.load(f)

	with open(map_indices_input, 'rb') as f:
		map_indices = pickle.load(f)

	with open(map_weights_input, 'rb') as f:
		map_weights = pickle.load(f)

	if len(map_indices) != scenarios:
		map_indices = map_indices[0]
		map_weights = map_weights[0]

	## GB: this gets hazard-consistent maps that we created from Miller's subsetting procedure
	sa_matrix_full = util.read_2dlist('input/sample_ground_motion_intensity_maps_road_only_filtered.txt',
									  delimiter='\t')
	sa_matrix = [sa_matrix_full[i] for i in
				 map_indices]  # GB: get the ground_motions for just the scenarios we are interested in

	lnsas = []
	magnitudes = []
	for row in sa_matrix:
		lnsas.append([log(float(sa)) for sa in row[4:]])
		magnitudes.append(float(row[2]))

	# no_damage_travel_time, no_damage_vmt, no_damage_trips_made, _ = precompute_network_performance()

	# get just the information relating to the strategy of interest
	if n_retrofits == 0 or n_retrofits == 71:
		# print strategy, strategy_index, n_retrofits, f_X_times.shape, f_X_times[strategy_index].shape

		average_travel_time, average_vmt, average_trips_made, average_direct_cost, average_delay_cost, average_connectivity_cost, \
		average_indirect_cost = compute_weighted_average_performance(lnsas, map_weights, num_damage_maps=10,
																	 travel_times=f_X_times[0],
																	 vmts=f_X_vmts[0],
																	 trips_made=f_X_trips[0],
																	 no_damage_travel_time=no_damage_travel_time,
																	 no_damage_vmt=no_damage_vmt,
																	 no_damage_trips_made=no_damage_trips_made,
																	 direct_costs=f_X_direct_costs[0])

		f_X_avg_time = average_travel_time
		f_X_avg_vmt = average_vmt
		f_X_avg_trip = average_trips_made
		f_X_exp_direct_cost = average_direct_cost
		f_X_exp_indirect_cost = average_indirect_cost  # hourly
		f_X_exp_cost = 24 * 125 * average_indirect_cost + average_direct_cost

	if not detailed:
		return f_X_exp_cost
	else:
		return f_X_delay_costs, f_X_conn_costs, f_X_indirect_costs, f_X_direct_costs, f_X_exp_cost

def import_retrofit_list(csv_filepath):

	strategy = []

	with open(csv_filepath, 'r') as csvfile:
		reader = csv.reader(csvfile, delimiter=',')
		count = 0
		for row in reader:
			# print row[0:]
			s = [x for x in row[0:] if x]  # filter out empty cells in the csv
			#print s
			strategy.append(s[0])
			count += 1
		# print 'count = ', count

	return strategy

def make_incremental_retrofit_lists(retrofit_list):

	ret_lists = []
	for i in range(1,len(retrofit_list)+1):
		ret_lists.append(retrofit_list[0:i])

	return ret_lists

def get_retrofit_cost(strategy, n_retrofits, unit_cost=293, rep_cost_ratio=0.25,
							  ret_cost_ratio=0.25): # given a strategy and the number of retrofits, compute the cost of all completed retrofits

	# originally, ret_cost_ratio = 0.30

	# Modify the name of the strategy to get the correct filepath.
	if strategy == 'oldest':
		strategy = 'age'
	elif strategy == 'busiest':
		strategy = 'traffic'
	elif strategy == 'weakest':
		strategy = 'fragility'
	elif strategy == 'OAT total':
		strategy = 'oat'
	elif strategy == 'Sobol exp. total':
		strategy = 'sobol'
		retrofit_list_filepath = 'sobol_input/sf_fullr_2020 revised averages rankings - ranking_370.csv'
	else:
		pass

	if strategy != 'sobol' and strategy != 'Sobol exp. total':
		retrofit_list_filepath = 'sobol_input/sf_fullr_2020 - ' + strategy + '_rets.csv'

	# print strategy, retrofit_list_filepath

	retrofit_list = import_retrofit_list(retrofit_list_filepath)
	all_ret_lists = make_incremental_retrofit_lists(retrofit_list)

	# get relevant list of retrofits for n_retrofits
	ret_list = all_ret_lists[n_retrofits-1]

	# print 'n_retrofits, ret_list', n_retrofits, ret_list

	bridge_dict, bridge_ids = get_sf_fullr_dict()

	retrofit_cost = 0
	for r in ret_list: # COPIED IN PART FROM compute_retrofit_cost() in mahmodel_road_only.py
		temp_cost = bridge_dict[r]['area'] * unit_cost * rep_cost_ratio * ret_cost_ratio
		retrofit_cost += temp_cost

	return retrofit_cost


def get_retrofit_results_from_file(n_retrofits, strategy, training=True):

	no_damage_travel_time, no_damage_vmt, no_damage_trips_made = load_individual_undamaged_stats()

	if not training:
		output_folder = 'sobol_output/retrofits/ret_revised_avg/r' + str(n_retrofits) + '/'
		scenarios = 45
		map_indices_input = 'sobol_input/sf_fullr_testing_map_indices.pkl'  # S = 45 for training sf_fullr
		map_weights_input = 'sobol_input/sf_fullr_testing_map_weights.pkl'  # S = 45 for training sf_fullr
	else:
		output_folder = 'sobol_output/retrofits/ret_revised_avg_s30/r' + str(n_retrofits) + '/'
		scenarios = 30
		map_indices_input = 'sobol_input/sf_fullr_training_map_indices.pkl'  # S = 30 for training sf_fullr
		map_weights_input = 'sobol_input/sf_fullr_training_map_weights.pkl'  # S = 30 for training sf_fullr

	strategies = ['oldest', 'busiest', 'weakest', 'composite', 'OAT total', 'OAT indirect', 'OAT direct', 'Sobol exp. total', 'Sobol exp. indirect', 'Sobol exp. direct']

	strategy_index = strategies.index(strategy)


	# n_scenarios = 45
	# dam_maps_per_scenario = 10

	# output_folder = 'sobol_output/retrofits/max_cost/max_cost_' + str(n_retrofits) + '/'
	if training: # S = 30
		output_folder = 'sobol_output/retrofits/s30/r' + str(n_retrofits) + '/'
	else: # S = 45
		output_folder = 'sobol_output/retrofits/rets_to_local/r' + str(n_retrofits) + '/'


	filename = '_sf_fullr'

	# store the results
	fX_times_output = output_folder + 'fX_times' + filename  # travel times for f_X
	fX_trips_output = output_folder + 'fX_trips' + filename  # trips made for f_X
	fX_vmts_output = output_folder + 'fX_vmts' + filename  # VMTs for f_X
	fX_avg_times_output = output_folder + 'fX_avg_time' + filename  # average TT
	fX_avg_trips_output = output_folder + 'fX_avg_trips' + filename  # average trips made
	fX_avg_vmts_output = output_folder + 'fX_avg_vmts' + filename  # average VMT
	fX_delay_costs_output = output_folder + 'fX_delay_costs' + filename
	fX_conn_costs_output = output_folder + 'fX_conn_costs' + filename
	fX_indirect_costs_output = output_folder + 'fX_indirect_costs' + filename
	fX_direct_costs_output = output_folder + 'fX_direct_costs' + filename
	fX_exp_indirect_cost_output = output_folder + 'fX_exp_indirect_costs' + filename
	fX_exp_direct_cost_output = output_folder + 'fX_exp_direct_costs' + filename
	fX_expected_cost_output = output_folder + 'fX_exp_costs' + filename

	damage_x_output = output_folder + 'damage_x' + filename

	# save data for f_X
	with open(damage_x_output, 'rb') as f:
		damage_tracker = pickle.load(f)

	with open(fX_times_output, 'rb') as f:  # save raw performance data
		f_X_times = pickle.load(f)
	with open(fX_trips_output, 'rb') as f:
		f_X_trips = pickle.load(f)
	with open(fX_vmts_output, 'rb') as f:
		f_X_vmts = pickle.load(f)

	with open(fX_avg_times_output, 'rb') as f:  # save average (expected) performance data
		f_X_avg_time = pickle.load(f)
	with open(fX_avg_trips_output, 'rb') as f:
		f_X_avg_trip = pickle.load(f)
	with open(fX_avg_vmts_output, 'rb') as f:
		f_X_avg_vmt = pickle.load(f)

	with open(fX_delay_costs_output, 'rb') as f:
		f_X_delay_costs = pickle.load(f)
	with open(fX_conn_costs_output, 'rb') as f:
		f_X_conn_costs = pickle.load(f)
	with open(fX_direct_costs_output, 'rb') as f:
		f_X_direct_costs = pickle.load(f)
	with open(fX_indirect_costs_output, 'rb') as f:
		f_X_indirect_costs = pickle.load(f)

	with open(fX_exp_direct_cost_output, 'rb') as f:
		f_X_exp_direct_cost = pickle.load(f)
	with open(fX_exp_indirect_cost_output, 'rb') as f:
		f_X_exp_indirect_cost = pickle.load(f)
	with open(fX_expected_cost_output, 'rb') as f:
		f_X_exp_cost = pickle.load(f)

	# n_samples = f_X_delay_costs.shape[0]
	# print 'shapes = ', f_X_indirect_costs.shape, f_X_times.shape

	with open(map_indices_input, 'rb') as f:
		map_indices = pickle.load(f)

	with open(map_weights_input, 'rb') as f:
		map_weights = pickle.load(f)

	if len(map_indices) != scenarios:
		map_indices = map_indices[0]
		map_weights = map_weights[0]

	## GB: this gets hazard-consistent maps that we created from Miller's subsetting procedure
	sa_matrix_full = util.read_2dlist('input/sample_ground_motion_intensity_maps_road_only_filtered.txt',
									  delimiter='\t')
	sa_matrix = [sa_matrix_full[i] for i in
				 map_indices]  # GB: get the ground_motions for just the scenarios we are interested in

	lnsas = []
	magnitudes = []
	for row in sa_matrix:
		lnsas.append([log(float(sa)) for sa in row[4:]])
		magnitudes.append(float(row[2]))

	# no_damage_travel_time, no_damage_vmt, no_damage_trips_made, _ = precompute_network_performance()

	# get just the information relating to the strategy of interest
	if n_retrofits == 0 or n_retrofits == 71:
		# print strategy, strategy_index, n_retrofits, f_X_times.shape, f_X_times[strategy_index].shape

		average_travel_time, average_vmt, average_trips_made, average_direct_cost, average_delay_cost, average_connectivity_cost, \
		average_indirect_cost = compute_weighted_average_performance(lnsas, map_weights, num_damage_maps=10,
																	 travel_times=f_X_times[strategy_index],
																	 vmts=f_X_vmts[strategy_index],
																	 trips_made=f_X_trips[strategy_index],
																	 no_damage_travel_time=no_damage_travel_time,
																	 no_damage_vmt=no_damage_vmt,
																	 no_damage_trips_made=no_damage_trips_made,
																	 direct_costs=f_X_direct_costs[strategy_index])

		f_X_avg_time = average_travel_time
		f_X_avg_vmt = average_vmt
		f_X_avg_trip = average_trips_made
		f_X_exp_direct_cost = average_direct_cost
		f_X_exp_indirect_cost = average_indirect_cost  # hourly
		f_X_exp_cost = 24 * 125 * average_indirect_cost + average_direct_cost

	# get just the information relating to the strategy of interest
	if n_retrofits > 0 and n_retrofits < 71:
		average_travel_time, average_vmt, average_trips_made, average_direct_cost, average_delay_cost, average_connectivity_cost, \
		average_indirect_cost = compute_weighted_average_performance(lnsas, map_weights, num_damage_maps=10,
																	 travel_times=f_X_times[strategy_index],
																	 vmts=f_X_vmts[strategy_index],
																	 trips_made=f_X_trips[strategy_index],
																	 no_damage_travel_time=no_damage_travel_time,
																	 no_damage_vmt=no_damage_vmt,
																	 no_damage_trips_made=no_damage_trips_made,
																	 direct_costs=f_X_direct_costs[strategy_index])

		f_X_avg_time = average_travel_time
		f_X_avg_vmt = average_vmt
		f_X_avg_trip = average_trips_made
		f_X_exp_direct_cost = average_direct_cost
		f_X_exp_indirect_cost = average_indirect_cost  # hourly
		f_X_exp_cost = 24 * 125 * average_indirect_cost + average_direct_cost

		f_X_delay_costs = f_X_delay_costs[strategy_index, :]
		f_X_conn_costs = f_X_conn_costs[strategy_index, :]
		f_X_direct_costs = f_X_direct_costs[strategy_index, :]
		f_X_indirect_costs = f_X_indirect_costs[strategy_index, :]
		# f_X_exp_cost = f_X_exp_cost[strategy_index]
		# print strategy, strategy_index, f_X_exp_cost

	return f_X_delay_costs, f_X_conn_costs, f_X_indirect_costs, f_X_direct_costs, f_X_exp_cost

def get_sobol_ret_cost_retrofit_results(n_retrofits, tag='ret_ratio_reverse'):


	output_folder = 'sobol_output/retrofits/' + tag + '/r' + str(n_retrofits) + '/'

	filename = '_sf_fullr'

	# store the results
	fX_times_output = output_folder + 'fX_times' + filename  # travel times for f_X
	fX_trips_output = output_folder + 'fX_trips' + filename  # trips made for f_X
	fX_vmts_output = output_folder + 'fX_vmts' + filename  # VMTs for f_X
	fX_avg_times_output = output_folder + 'fX_avg_time' + filename  # average TT
	fX_avg_trips_output = output_folder + 'fX_avg_trips' + filename  # average trips made
	fX_avg_vmts_output = output_folder + 'fX_avg_vmts' + filename  # average VMT
	fX_delay_costs_output = output_folder + 'fX_delay_costs' + filename
	fX_conn_costs_output = output_folder + 'fX_conn_costs' + filename
	fX_indirect_costs_output = output_folder + 'fX_indirect_costs' + filename
	fX_direct_costs_output = output_folder + 'fX_direct_costs' + filename
	fX_exp_indirect_cost_output = output_folder + 'fX_exp_indirect_costs' + filename
	fX_exp_direct_cost_output = output_folder + 'fX_exp_direct_costs' + filename
	fX_expected_cost_output = output_folder + 'fX_exp_costs' + filename

	damage_x_output = output_folder + 'damage_x' + filename

	# save data for f_X
	with open(damage_x_output, 'rb') as f:
		damage_tracker = pickle.load(f)

	with open(fX_times_output, 'rb') as f:  # save raw performance data
		f_X_times = pickle.load(f)
	with open(fX_trips_output, 'rb') as f:
		f_X_trips = pickle.load(f)
	with open(fX_vmts_output, 'rb') as f:
		f_X_vmts = pickle.load(f)

	with open(fX_avg_times_output, 'rb') as f:  # save average (expected) performance data
		f_X_avg_time = pickle.load(f)
	with open(fX_avg_trips_output, 'rb') as f:
		f_X_avg_trip = pickle.load(f)
	with open(fX_avg_vmts_output, 'rb') as f:
		f_X_avg_vmt = pickle.load(f)

	with open(fX_delay_costs_output, 'rb') as f:
		f_X_delay_costs = pickle.load(f)
	with open(fX_conn_costs_output, 'rb') as f:
		f_X_conn_costs = pickle.load(f)
	with open(fX_direct_costs_output, 'rb') as f:
		f_X_direct_costs = pickle.load(f)
	with open(fX_indirect_costs_output, 'rb') as f:
		f_X_indirect_costs = pickle.load(f)

	with open(fX_exp_direct_cost_output, 'rb') as f:
		f_X_exp_direct_cost = pickle.load(f)
	with open(fX_exp_indirect_cost_output, 'rb') as f:
		f_X_exp_indirect_cost = pickle.load(f)
	with open(fX_expected_cost_output, 'rb') as f:
		f_X_exp_cost = pickle.load(f)

	# with open(map_indices_input, 'rb') as f:
	# 	map_indices = pickle.load(f)
	#
	# with open(map_weights_input, 'rb') as f:
	# 	map_weights = pickle.load(f)
	#
	# if len(map_indices) != scenarios:
	# 	map_indices = map_indices[0]
	# 	map_weights = map_weights[0]
	#
	# ## GB: this gets hazard-consistent maps that we created from Miller's subsetting procedure
	# sa_matrix_full = util.read_2dlist('input/sample_ground_motion_intensity_maps_road_only_filtered.txt',
	# 								  delimiter='\t')
	# sa_matrix = [sa_matrix_full[i] for i in
	# 			 map_indices]  # GB: get the ground_motions for just the scenarios we are interested in
	#
	# lnsas = []
	# magnitudes = []
	# for row in sa_matrix:
	# 	lnsas.append([log(float(sa)) for sa in row[4:]])
	# 	magnitudes.append(float(row[2]))
	#
	# # no_damage_travel_time, no_damage_vmt, no_damage_trips_made, _ = precompute_network_performance()
	#
	# # get just the information relating to the strategy of interest
	# if n_retrofits == 0 or n_retrofits == 71:
	# 	# print strategy, strategy_index, n_retrofits, f_X_times.shape, f_X_times[strategy_index].shape
	#
	# 	average_travel_time, average_vmt, average_trips_made, average_direct_cost, average_delay_cost, average_connectivity_cost, \
	# 	average_indirect_cost = compute_weighted_average_performance(lnsas, map_weights, num_damage_maps=10,
	# 																 travel_times=f_X_times[0],
	# 																 vmts=f_X_vmts[0],
	# 																 trips_made=f_X_trips[0],
	# 																 no_damage_travel_time=no_damage_travel_time,
	# 																 no_damage_vmt=no_damage_vmt,
	# 																 no_damage_trips_made=no_damage_trips_made,
	# 																 direct_costs=f_X_direct_costs[0])
	#
	# 	f_X_avg_time = average_travel_time
	# 	f_X_avg_vmt = average_vmt
	# 	f_X_avg_trip = average_trips_made
	# 	f_X_exp_direct_cost = average_direct_cost
	# 	f_X_exp_indirect_cost = average_indirect_cost  # hourly
	# 	f_X_exp_cost = 24 * 125 * average_indirect_cost + average_direct_cost


	return f_X_delay_costs, f_X_conn_costs, f_X_indirect_costs, f_X_direct_costs, f_X_exp_cost

def plot_retrofit_results_vs_cumulative_retrofit_cost(tag):
	# tag can be 'ret_cost', 'ret_ratio', 'ret_ratio_reverse'

	training = False

	# folder = 'figs_diff_p/'
	# folder = 'exp_cost_vs_cum_cost/'
	folder = 'figs/'

	fig_folder = folder

	# Get baseline, i.e. result when R = 0.
	baseline = get_baseline_retrofit_results(n_retrofits=0, training=training)  # expected cost when R = 0 (no retrofit costs)
	# Get best case, i.e. result when R = 71. (all bridges retrofitted)
	best_exp_cost = get_baseline_retrofit_results(n_retrofits=71, training=training)
	best_ret_cost = get_retrofit_cost('age', n_retrofits=71)

	print 'baseline, best_exp_cost, best_ret_cost = ', baseline, best_exp_cost, best_ret_cost

	# first, get all the data for the other strategies -- these are for S = 45
	strategies = ['oldest', 'busiest', 'weakest', 'composite', 'OAT total', 'Sobol exp. total']
	n_retrofits = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 71]
	R = len(n_retrofits)
	exp_cost_results = {} # y axis
	ret_cost_results = {} # x axis
	for s in strategies:
		temp_exp_results = np.zeros((R,))
		temp_ret_results = np.zeros((R,))
		i = 0
		for r in n_retrofits:
			if r == 0:
				temp_exp_results[i] = baseline
				temp_ret_results[i] = 0
			elif r == 71:
				temp_exp_results[i] = best_exp_cost
				temp_ret_results[i] = best_ret_cost
			else:
				_, _, _, _ , expected_total_cost = get_retrofit_results_from_file(r, strategy=s,training=training)
				retrofit_cost = get_retrofit_cost(s, r)
				temp_exp_results[i] = expected_total_cost
				temp_ret_results[i] = retrofit_cost

			i += 1

		# exp_cost_results[s] = [(t-baseline)/baseline*-100 for t in temp_exp_results] # -100 is to convert to a positive percentage
		ret_cost_results[s] = [t for t in temp_ret_results]
		exp_cost_results[s] = [-1*compute_percent_reduction(baseline, t) for t in temp_exp_results]

	# Now get the results of the Sobol' index-based strategy with respect to the total cost, INCLUDING retrofit costs.
	sobol_exp_results = np.zeros((len(n_retrofits),))
	sobol_ret_results = np.zeros((len(n_retrofits),))
	i = 0
	for r in n_retrofits:
		if r == 0:
			# sobol_exp_results[i] = baseline
			sobol_exp_results[i] = 0
			sobol_ret_results[i] = 0
		elif r == 71:
			sobol_exp_results[i] = (best_exp_cost-baseline)/baseline*-100
			sobol_ret_results[i] = best_ret_cost
		else:
			f_X_delay_costs, f_X_conn_costs, f_X_indirect_costs, f_X_direct_costs, expected_total_cost = get_sobol_ret_cost_retrofit_results(
				n_retrofits=r, tag=tag)
			retrofit_cost = get_retrofit_cost(tag, n_retrofits=r)
			# sobol_exp_results[i] = (expected_total_cost-baseline)/baseline*-100
			sobol_exp_results[i] = -1*compute_percent_reduction(baseline, expected_total_cost)
			sobol_ret_results[i] = retrofit_cost
			# print i, r, retrofit_cost
		i += 1

	print 'oldest results = ', exp_cost_results['oldest']
	print 'sobol results = ', sobol_exp_results
	print 'sobol ret costs = ', sobol_ret_results

	title = 'sf_fullr_exp_cost_reduction_vs_ret_cost_' + tag
	# title = 'sf_fullr_exp_cost_vs_ret_cost'

	strategy_labels = ['age', 'traffic volume', 'fragility', 'composite', 'OAT', 'Sobol, $\\mathbb{E}[C]$']

	marker_style = 'o'

	fig = plt.figure()
	ax = fig.add_subplot(111)
	i = 0
	for s in strategies:
		ax.plot(ret_cost_results[s], exp_cost_results[s], color=get_color(s), marker=marker_style, label=strategy_labels[i])

		i += 1
	ax.plot(sobol_ret_results, sobol_exp_results, marker='s', ls='--', color=get_color('Sobol'),
			label='Sobol, cost ratio')
	ax.plot(ret_cost_results['oldest'], [exp_cost_results['oldest'][-1] for i in ret_cost_results['oldest']], color='black', label='all retrofitted')
	ax.set_xlim([0, best_ret_cost+1e7])
	# plt.legend(loc='best', prop={'size': 8})
	# plt.xlabel('Number of retrofitted bridges, $R$')
	plt.xlabel('Cumulative retrofit cost, $C_R$ [USD]')
	plt.ylabel('$\\%$ reduction in $\\mathbb{E}[C]$  of network performance')
	ax.spines['top'].set_visible(False)
	ax.spines['right'].set_visible(False)
	ax.yaxis.set_ticks_position('left')
	ax.xaxis.set_ticks_position('bottom')

	handles, labels = ax.get_legend_handles_labels()

	handles = [handles[7], handles[6], handles[5], handles[0], handles[4], handles[1], handles[3], handles[2]]
	labels = [labels[7], labels[6], labels[5], labels[0], labels[4], labels[1], labels[3], labels[2]]

	ax.legend(handles, labels, prop={'size':11}, loc='upper left', bbox_to_anchor = (0.2, 0.4, 0.5, 0.5), frameon=False)
	plt.savefig(fig_folder + title + '.png', bbox_inches='tight')


	# ZOOMED in version
	fig = plt.figure()
	ax = fig.add_subplot(111)
	i = 0
	for s in strategies:
		ax.plot(ret_cost_results[s], exp_cost_results[s], color=get_color(s), marker=marker_style,
				label=strategy_labels[i])

		i += 1
	ax.plot(sobol_ret_results, sobol_exp_results, marker='s', ls='--', color=get_color('Sobol'),
			label='Sobol, cost ratio')
	ax.plot(ret_cost_results['oldest'], [exp_cost_results['oldest'][-1] for i in ret_cost_results['oldest']],
			color='black', label='all retrofitted')
	ax.set_xlim([0, sobol_ret_results[15]])
	ax.set_ylim([0, sobol_exp_results[15]])
	# plt.legend(loc='best', prop={'size': 8})
	# plt.xlabel('Number of retrofitted bridges, $R$')
	plt.xlabel('Cumulative retrofit cost, $C_R$ [USD]')
	plt.ylabel('$\\%$ reduction in $\\mathbb{E}[C]$  of network performance')
	ax.spines['top'].set_visible(False)
	ax.spines['right'].set_visible(False)
	ax.yaxis.set_ticks_position('left')
	ax.xaxis.set_ticks_position('bottom')

	handles, labels = ax.get_legend_handles_labels()

	handles = [handles[7], handles[6], handles[5], handles[0], handles[4], handles[1], handles[3], handles[2]]
	labels = [labels[7], labels[6], labels[5], labels[0], labels[4], labels[1], labels[3], labels[2]]

	ax.legend(handles, labels, prop={'size': 12}, loc='upper left', bbox_to_anchor = (0.01, 0.4, 0.5, 0.5), frameon=False)
	plt.savefig(fig_folder + title + '_zoom.png', bbox_inches='tight')

	plt.show()


def main():
	plot_retrofit_results_vs_cumulative_retrofit_cost('ret_ratio_reverse')

if __name__ == "__main__":
	# execute only if run as a script
	main()