//	Subnetwork Origin-Destination Trip Table Estimation Tool
//
//	Chi Xie, 08-11-09.
//	Copyright (c) 2009 The University of Texas at Austin
//	$Revision: 1.3$  $Date: 2011/01/06       :  :   $
//	This program is an implementation of the Frank-Wolfe method to compute the 
//	origin-destination trip table of an abstract network from the link flows in 
//	the format generated by Bar-Gera's OBA code

//	step 1 >> read the original network file
//	step 2 >> read the abstract network file
//	step 3 >> identify nodes included in the abstract network
//	step 4 >> read link flows from the original network flow file
//	step 5 >> compute the O-D tables of the abstract network
//	step 6 >> write the O-D table into the output files

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <new>
using namespace std;

void matrix_inversion(double **, int, double **);
double matrix_determinant(double **, int);
bool bellman_ford_shortest_path(int, int, long int *[][2], double *, int, long int *, long int **, double *);
void floyd_warshall_shortest_path(int, long int *, long int ***, double **);

const int max_num_node = 1000;
const int max_num_link = 5000;
const int max_num_route = 100000;
const double negative_infinity = -999999999;
const double positive_infinity = 999999999999999;
const double positive_small = 0.01;

int main(int argc, char *argv[]) {

	cout << "\n";
	cout << "   Origin-Destination Travel Table Estimation Tool for Sketch Planning\n";
	cout << "   Designed and Coded by Chi Xie, Ph.D.\n";
	cout << "   Version 1.4, Updated on March 31, 2011\n";
	cout << "   (C) Center for Transportation Research, The University of Texas at Austin\n";
	cout << "\n";
	int starttime, endtime, cputime;
	starttime = clock() / CLOCKS_PER_SEC;
	cout << "   Clock" << "\t\t" << starttime << " : Starting" << "\n";

	// step 1 >> read the original network file

	if (!argv[1])
	{
		cout << "argument(s) missing...";
		return 1;
	}
	ifstream network_in(argv[1]);

	char metadata[50];
	char header[200];
	char c_link_row[500];
	long int link_id[max_num_link][2];
	double link_attrbt[max_num_link][8];
	
	for (int i = 0; i < 7; i++)
	{
		network_in.getline(metadata, 50);
	}
	network_in.getline(header, 200);
	int l = 0;
	while(!network_in.getline(c_link_row, 500).eof())
	{
		char *pEnd;
		link_id[l][0] = strtol(c_link_row, &pEnd, 10);
		link_id[l][1] = strtol(pEnd, &pEnd, 10);
		link_attrbt[l][0] = strtod(pEnd, &pEnd);
		link_attrbt[l][1] = strtod(pEnd, &pEnd);
		link_attrbt[l][2] = strtod(pEnd, &pEnd);
		link_attrbt[l][3] = strtod(pEnd, &pEnd);
		link_attrbt[l][4] = strtod(pEnd, &pEnd);
		link_attrbt[l][5] = strtod(pEnd, &pEnd);
		link_attrbt[l][6] = strtod(pEnd, &pEnd);
		link_attrbt[l][7] = strtod(pEnd, &pEnd);
		/*
		cout << link_id[l][0] << "\t" << link_id[l][1] << "\t" 
			<< link_attrbt[l][0] << "\t" << link_attrbt[l][1] << "\t" 
			<< link_attrbt[l][2] << "\t" << link_attrbt[l][3] << "\t" 
			<< link_attrbt[l][4] << "\t" << link_attrbt[l][5] << "\t" 
			<< link_attrbt[l][6] << "\t" << link_attrbt[l][7] << "\n";
		*/
		l++;
	}
	int num_link = l;
	network_in.close();

	long int node_id[max_num_node];
	int p = 0;
	bool node_existence = false;
	for (int m = 0; m < num_link; m++)
	{
		for (int n = 0; n < p; n++)
		{
			if (link_id[m][0] == node_id[n])
			{
				node_existence = true;
				break;
			}
		}
		if (node_existence == false)
		{
			node_id[p] = link_id[m][0];
			p++;
		}
		node_existence = false;
	}
	for (int m = 0; m < num_link; m++)
	{
		for (int n = 0; n < p; n++)
		{
			if (link_id[m][1] == node_id[n])
			{
				node_existence = true;
				break;
			}
		}
		if (node_existence == false)
		{
			node_id[p] = link_id[m][1];
			p++;
		}
		node_existence = false;
	}
	int num_node = p;

	
	// step 2 >> read the abstract network file
	
	if (!argv[2])
	{
		cout << "argument(s) missing...";
		return 1;
	}
	ifstream abs_network_in(argv[2]);

	long int link_id_abs[max_num_link][2];
	long int *link_address_abs[max_num_link][2];
	double link_attrbt_abs[max_num_link][10];
	
	/*
	while(1)
	{
		abs_network_in.getline(metadata, 50);
		abs_network_in.get(header, 2);
		if (!strcmp(header, "~")) break;
	}
	*/
	for (int i = 0; i < 7; i++)
	{
		abs_network_in.getline(metadata, 50);
	}

	abs_network_in.getline(header, 200);
	l = 0;
	while(!abs_network_in.getline(c_link_row, 500).eof())
	{
		char *pEnd;
		link_id_abs[l][0] = strtol(c_link_row, &pEnd, 10); // initial id
		link_id_abs[l][1] = strtol(pEnd, &pEnd, 10); // terminal id
		link_attrbt_abs[l][0] = strtod(pEnd, &pEnd); // capacity
		link_attrbt_abs[l][1] = strtod(pEnd, &pEnd); // length
		link_attrbt_abs[l][2] = strtod(pEnd, &pEnd); // free flow time
		link_attrbt_abs[l][3] = strtod(pEnd, &pEnd); // base
		link_attrbt_abs[l][4] = strtod(pEnd, &pEnd); // power
		link_attrbt_abs[l][5] = strtod(pEnd, &pEnd); // speed limit
		link_attrbt_abs[l][6] = strtod(pEnd, &pEnd); // toll
		link_attrbt_abs[l][7] = strtod(pEnd, &pEnd); // type
		link_attrbt_abs[l][8] = 0; // volume
		link_attrbt_abs[l][9] = 0; // cost
		/*
		cout << link_id_abs[l][0] << "\t" << link_id_abs[l][1] << "\t" 
			<< link_attrbt_abs[l][0] << "\t" << link_attrbt_abs[l][1] << "\t" 
			<< link_attrbt_abs[l][2] << "\t" << link_attrbt_abs[l][3] << "\t" 
			<< link_attrbt_abs[l][4] << "\t" << link_attrbt_abs[l][5] << "\t" 
			<< link_attrbt_abs[l][6] << "\t" << link_attrbt_abs[l][7] << "\n";
		*/
		l++;
	}
	int num_link_abs = l;
	abs_network_in.close();

	long int node_id_abs[max_num_node];
	p = 0;
	node_existence = false;
	for (int m = 0; m < num_link_abs; m++)
	{
		for (int n = 0; n <= p; n++)
		{
			if (link_id_abs[m][0] == node_id_abs[n])
			{
				node_existence = true;
				break;
			}
		}
		if (node_existence == false)
		{
			node_id_abs[p] = link_id_abs[m][0];
			p++;
		}
		node_existence = false;
	}
	for (int m = 0; m < num_link_abs; m++)
	{
		for (int n = 0; n <= p; n++)
		{
			if (link_id_abs[m][1] == node_id_abs[n])
			{
				node_existence = true;
				break;
			}
		}
		if (node_existence == false)
		{
			node_id_abs[p] = link_id_abs[m][1];
			p++;
		}
		node_existence = false;
	}
	int num_node_abs = p;
	for (int m = 0; m < num_link_abs; m++)
	{
		for (int n = 0; n <= num_node_abs; n++)
		{
			if (link_id_abs[m][0] == node_id_abs[n])
			{
				link_address_abs[m][0] = &node_id_abs[n];
			}
			if (link_id_abs[m][1] == node_id_abs[n])
			{
				link_address_abs[m][1] = &node_id_abs[n];
			}
		}
	}


	// step 3 >> identify nodes included in the abstract network

	bool *node_id_abs_in;
	node_id_abs_in = new bool [num_node_abs];
	for (int n1 = 0; n1 < num_node; n1++)
	{
		node_id_abs_in[n1] = false;
		for (int n2 = 0; n2 < num_node_abs; n2++)
		{
			if (node_id[n1] == node_id_abs[n2])
			{
				node_id_abs_in[n1] = true;
				break;
			}
		}
	}


	// step 4 >> read link flows from the original network flow file
	
	if (!argv[3])
	{
		cout << "argument(s) missing...";
		return 1;
	}
	ifstream linkflow_in(argv[3]);
	
	long int fr_node_id, to_node_id;
	for (int i = 0; i < 5; i++)
	{
		linkflow_in.getline(metadata, 50);
	}
	linkflow_in.getline(header, 200);
	while(!linkflow_in.getline(c_link_row, 500).eof())
	{
		char *pEnd;
		fr_node_id = strtol(c_link_row, &pEnd, 10);
		to_node_id = strtol(pEnd, &pEnd, 10);
		for (int m = 0; m < num_link; m++)
		{
			if (fr_node_id == link_id_abs[m][0] && to_node_id == link_id_abs[m][1])
			{
				pEnd = pEnd + 3; // move the pointer pEnd forward by 3 addresses
				link_attrbt_abs[m][8] = strtod(pEnd, &pEnd); // volume
				link_attrbt_abs[m][9] = strtod(pEnd, &pEnd); // cost
				/*
				cout << link_id_abs[m][0] << "\t" << link_id_abs[m][1] << "\t" 
					<< link_attrbt_abs[m][8] << "\t" << link_attrbt_abs[m][9] << "\n";
				*/
			}
		}
	}
	linkflow_in.close();


	// step 5 >> compute the O-D tables of the abstract network

	double **od_flow_abs, **pre_od_flow_abs, **left_od_flow_abs, **right_od_flow_abs, **mid_od_flow_abs, **od_imp_abs, **aux_od_imp_abs, **od_cost_abs;
	od_flow_abs = new double *[num_node_abs];
	pre_od_flow_abs = new double *[num_node_abs];
	left_od_flow_abs = new double *[num_node_abs];
	right_od_flow_abs = new double *[num_node_abs];
	mid_od_flow_abs = new double *[num_node_abs];
	od_imp_abs = new double *[num_node_abs];
	aux_od_imp_abs = new double *[num_node_abs];
	od_cost_abs = new double *[num_node_abs];
	for (int n = 0; n < num_node_abs; n++)
	{
		od_flow_abs[n] = new double [num_node_abs];
		pre_od_flow_abs[n] = new double [num_node_abs];
		left_od_flow_abs[n] = new double [num_node_abs];
		right_od_flow_abs[n] = new double [num_node_abs];
		mid_od_flow_abs[n] = new double [num_node_abs];
		od_imp_abs[n] = new double [num_node_abs];
		aux_od_imp_abs[n] = new double [num_node_abs];
		od_cost_abs[n] = new double [num_node_abs];
	}
	//
	int num_basic_path_abs = num_link_abs;
	long int ***basic_path_address_abs;
	basic_path_address_abs = new long int **[num_basic_path_abs];
	for (int k = 0; k < num_basic_path_abs; k++)
	{
		basic_path_address_abs[k] = new long int *[2];
	}
	double *basic_path_flow_abs, *basic_path_imp_abs, *link_cost_abs, *link_multiplier_abs;
	basic_path_flow_abs = new double [num_basic_path_abs];
	basic_path_imp_abs = new double [num_basic_path_abs];
	link_cost_abs = new double [num_link_abs];
	link_multiplier_abs = new double [num_link_abs];
	//
	double *entering_col_abs;
	entering_col_abs = new double [num_link_abs];
	double entering_flow_abs, temp_entering_flow_abs, min_flow_abs;
	double product_inv_matrix_link_flow, product_inv_matrix_entering_col;
	//
	double **basic_matrix, **inv_basic_matrix;
	basic_matrix = new double *[num_link_abs];
	inv_basic_matrix = new double *[num_link_abs];
	for (int m = 0; m < num_link_abs; m++)
	{
		basic_matrix[m] = new double [num_link_abs];
		inv_basic_matrix[m] = new double [num_link_abs];
	}



	double **identity_matrix;
	identity_matrix = new double *[num_link_abs];
	for (int m = 0; m < num_link_abs; m++)
	{
		identity_matrix[m] = new double [num_link_abs];
	}




	//
	long int **predecessor_node_id_abs;
	predecessor_node_id_abs = new long int *[num_node_abs];
	double *node_cost_abs;
	node_cost_abs = new double [num_node_abs];
	//
	long int ***od_predecessor_node_id_abs;
	od_predecessor_node_id_abs = new long int **[num_node_abs];
	for (int n = 0; n < num_node_abs; n++)
	{
		od_predecessor_node_id_abs[n] = new long int *[num_node_abs];
	}
	//
	int fr_node_index, to_node_index, ori_node_index, des_node_index, leaving_col_index;
	bool loop_existence, basic_path_existence;
	double min_reduced_cost;

	cout << "   Estimate the travel table in the abstract network\n";
	
	// step 5.0 >> initialization
	for (int fr = 0; fr < num_node_abs; fr++)
	{
		for (int to = 0; to < num_node_abs; to++)
		{
			od_flow_abs[fr][to] = 0; // zero O-D flow rates
			od_imp_abs[fr][to] = negative_infinity; // infinite O-D entropy impedances
			if (fr == to)
			{
				od_imp_abs[fr][to] = 0; // zero O-D entropy impedances
			}
		}
	}
	for (int fr = 0; fr < num_node_abs; fr++)
	{
		for (int to = 0; to < num_node_abs; to++)
		{
			for (int m = 0; m < num_link_abs; m++)
			{
				if (node_id_abs[fr] == link_id_abs[m][0] && node_id_abs[to] == link_id_abs[m][1])
				{
					if  (link_attrbt_abs[m][8] == 0)
					{
						link_attrbt_abs[m][8] = 0.01;
					}
					od_flow_abs[fr][to] = link_attrbt_abs[m][8]; // initial O-D flow rates
					od_imp_abs[fr][to] = log(od_flow_abs[fr][to]); // initial O-D entropy impedances
				}
			}
		}
	}
	int num_intervals = 20;
	double derivative, accml_error, avg_error;
	double epsilon = 0.001;
	cout << "\n";
	// steps 5.1-.4 >> direction finding, line search, solution update and convergence test
	int iteration = 0;
	int pivot = 0;
	do
	{
		iteration++;
		cout << "   ==========================================\n";
		cout << "   Iteration: " << iteration << "\n";
		cout << "   Pivot move over the linear space ";
		pivot = 0;
		for (int fr = 0; fr < num_node_abs; fr++)
		{
			for (int to = 0; to < num_node_abs; to++)
			{
				pre_od_flow_abs[fr][to] = od_flow_abs[fr][to];
			}
		}
		// step 5.1 >> direction finding (the linearized problem)
		// initialize the basic path variables
		for (int k = 0; k < num_basic_path_abs; k++)
		{
			basic_path_address_abs[k][0] = link_address_abs[k][0];
			basic_path_address_abs[k][1] = link_address_abs[k][1];
			fr_node_index = basic_path_address_abs[k][0] - node_id_abs;
			to_node_index = basic_path_address_abs[k][1] - node_id_abs;
			basic_path_flow_abs[k] = link_attrbt_abs[k][8];
			basic_path_imp_abs[k] = od_imp_abs[fr_node_index][to_node_index];
		}
		// initialize the basic link-path incidence matrix
		for (int m = 0; m < num_link_abs; m++)
		{
			for (int k = 0; k < num_basic_path_abs; k++)
			{
				if (m == k)
				{
					basic_matrix[m][k] = 1;
				}
				else
				{
					basic_matrix[m][k] = 0;
				}
			}
		}
		do
		{
			pivot++;
			cout << ".";
			matrix_inversion(basic_matrix, num_link_abs, inv_basic_matrix);
			int b = 0;//********************
			for (int m = 0; m < num_link_abs; m++)
			{
				for (int n = 0; n < num_link_abs; n++)
				{
					int a = 0;
					for (int k = 0; k < num_basic_path_abs; k++)
					{
						a += basic_matrix[m][k] * inv_basic_matrix[k][n];
					}
					identity_matrix[m][n] = a;
					if (m == n && a != 1)
					{
						b = 1;
					}
					if (m != n && a != 0)
					{
						b = 1;
					}
				}
			}
			if (b == 1)
			{
				for (int m = 0; m < num_link_abs; m++)
				{
					for (int n = 0; n < num_link_abs; n++)
					{
						cout << identity_matrix[m][n] << " ";
					}
					cout << "\n";
				}
				cout << "\n";
			}//********************
			for (int m = 0; m < num_link_abs; m++)
			{
				link_multiplier_abs[m] = 0;
				for (int k = 0; k < num_basic_path_abs; k++)
				{
					link_multiplier_abs[m] += basic_path_imp_abs[k] * inv_basic_matrix[k][m];
				}
			}
			// implementation of Bellman-Ford shortest path algorithm
			for (int n = 0; n < num_node_abs; n++)
			{
				node_cost_abs[n] = positive_infinity;
			}
			int fr = 0;
			node_cost_abs[fr] = 0;
			for (int m = 0; m < num_link_abs; m++)
			{
				link_cost_abs[m] = -1 * link_multiplier_abs[m];
			}
			loop_existence = bellman_ford_shortest_path(fr, num_link_abs, link_address_abs, link_cost_abs, 
				num_node_abs, node_id_abs, predecessor_node_id_abs, node_cost_abs);
			// implementation of Floyd-Warshall shortest path algorithm
			if (loop_existence == true)
			{
				for (int m = 0; m < num_link_abs; m++)
				{
					if (link_cost_abs[m] < 0)
					{
						link_cost_abs[m] = positive_small;
					}
				}
			}
			for (int n1 = 0; n1 < num_node_abs; n1++)
			{
				for (int n2 = 0; n2 < num_node_abs; n2++)
				{
					if (n1 == n2)
					{
						od_cost_abs[n1][n2] = 0;
						od_predecessor_node_id_abs[n1][n2] = &node_id_abs[n1];
					}
					else
					{
						od_cost_abs[n1][n2] = positive_infinity;
					}
				}
			}
			for (int m = 0; m < num_link_abs; m++)
			{
				fr_node_index = link_address_abs[m][0] - node_id_abs;
				to_node_index = link_address_abs[m][1] - node_id_abs;
				od_cost_abs[fr_node_index][to_node_index] = link_cost_abs[m];
				od_predecessor_node_id_abs[fr_node_index][to_node_index] = link_address_abs[m][0];
			}
			floyd_warshall_shortest_path(num_node_abs, node_id_abs, od_predecessor_node_id_abs, od_cost_abs);
			// search for the path with the minimum reduced cost (the entering path)
			min_reduced_cost = positive_infinity;
			for (int n1 = 0; n1 < num_node_abs; n1++)
			{
				for (int n2 = 0; n2 < num_node_abs; n2++)
				{
					for (int m = 0; m < num_link_abs; m++)
					{
					//	if (node_id_abs[n1] == link_id_abs[m][0] && node_id_abs[n2] == link_id_abs[m][1] && link_attrbt_abs[m][8] == 0)
					//	{
					//	}
					//	else
					//	{
							if (od_imp_abs[n1][n2] + od_cost_abs[n1][n2] < min_reduced_cost)
							{
								// retrieve the candidate entering path
								ori_node_index = n1;
								des_node_index = n2;
								for (int m = 0; m < num_link_abs; m++)
								{
									entering_col_abs[m] = 0;
								}
								fr_node_index = des_node_index;
								do
								{
									to_node_index = fr_node_index;
									fr_node_index = od_predecessor_node_id_abs[ori_node_index][to_node_index] - node_id_abs;
									for (int m = 0; m < num_link_abs; m++)
									{
										if (link_id_abs[m][0] == node_id_abs[fr_node_index] && link_id_abs[m][1] == node_id_abs[to_node_index])
										{
											entering_col_abs[m] = 1;
											break;
										}
									}
								}
								while (ori_node_index != fr_node_index);
								// check the overlap between the entering path and the basic paths
								for (int k = 0; k < num_basic_path_abs; k++)
								{
									basic_path_existence = false;
									for (int m = 0; m < num_link_abs; m++)
									{
										if (basic_matrix[m][k] == entering_col_abs[m])
										{
											basic_path_existence = true;
										}
										else
										{
											basic_path_existence = false;
											break;
										}
									}
									if (basic_path_existence == true)
									{
										break;
									}
								}
								if (basic_path_existence == false)
								{
									min_reduced_cost = od_imp_abs[n1][n2] + od_cost_abs[n1][n2];
									ori_node_index = n1;
									des_node_index = n2;
								}
							}
					//	}
					}
				}
			}
			if (min_reduced_cost < 0 && pivot < 500)
			{
				// retrieve the link-path incidence vector of the entering path
				for (int m = 0; m < num_link_abs; m++)
				{
					entering_col_abs[m] = 0;
				}
				fr_node_index = des_node_index;
				do
				{
					to_node_index = fr_node_index;
					fr_node_index = od_predecessor_node_id_abs[ori_node_index][to_node_index] - node_id_abs;
					for (int m = 0; m < num_link_abs; m++)
					{
						if (link_id_abs[m][0] == node_id_abs[fr_node_index] && link_id_abs[m][1] == node_id_abs[to_node_index])
						{
							entering_col_abs[m] = 1;
							break;
						}
					}
				}
				while (ori_node_index != fr_node_index);
				// determine the entering path flow rate
				entering_flow_abs = positive_infinity;
				for (int k = 0; k < num_basic_path_abs; k++)
				{
					product_inv_matrix_link_flow = 0;
					product_inv_matrix_entering_col = 0;
					for (int m = 0; m < num_link_abs; m++)
					{
						product_inv_matrix_link_flow += inv_basic_matrix[k][m] * link_attrbt_abs[m][8];
						product_inv_matrix_entering_col += inv_basic_matrix[k][m] * entering_col_abs[m];
					}
					if (product_inv_matrix_entering_col > 0)
					{
						temp_entering_flow_abs = product_inv_matrix_link_flow / product_inv_matrix_entering_col;
						if (temp_entering_flow_abs >= 0 && temp_entering_flow_abs < entering_flow_abs)
						{
							entering_flow_abs = temp_entering_flow_abs;
						}
					}
				}
				// update the basic path variables
				for (int k = 0; k < num_basic_path_abs; k++)
				{
					basic_path_flow_abs[k] = 0;
					for (int m = 0; m < num_link_abs; m++)
					{
						basic_path_flow_abs[k] += inv_basic_matrix[k][m] * 
							(link_attrbt_abs[m][8] - entering_col_abs[m] * entering_flow_abs);
					}
				}
				min_flow_abs = positive_infinity;
				for (int k = 0; k < num_basic_path_abs; k++)
				{
					if (min_flow_abs >= basic_path_flow_abs[k])
					{
						min_flow_abs = basic_path_flow_abs[k];
						leaving_col_index = k;
					}
				}
				basic_path_flow_abs[leaving_col_index] = entering_flow_abs;
				basic_path_address_abs[leaving_col_index][0] = &node_id_abs[ori_node_index];
				basic_path_address_abs[leaving_col_index][1] = &node_id_abs[des_node_index];
				for (int k = 0; k < num_basic_path_abs; k++)
				{
					fr_node_index = basic_path_address_abs[k][0] - node_id_abs;
					to_node_index = basic_path_address_abs[k][1] - node_id_abs;
					basic_path_imp_abs[k] = od_imp_abs[fr_node_index][to_node_index];
				}
				// update the O-D flow matrix
				for (int fr = 0; fr < num_node_abs; fr++)
				{
					for (int to = 0; to < num_node_abs; to++)
					{
						od_flow_abs[fr][to] = 0;
					}
				}
				for (int k = 0; k < num_basic_path_abs; k++)
				{
					fr_node_index = basic_path_address_abs[k][0] - node_id_abs;
					to_node_index = basic_path_address_abs[k][1] - node_id_abs;
					od_flow_abs[fr_node_index][to_node_index] += basic_path_flow_abs[k];
				}
				// update the basic link-path incidence matrix
				for (int m = 0; m < num_link_abs; m++)
				{
					basic_matrix[m][leaving_col_index] = entering_col_abs[m];
				}
				// check the invertibility of the basic link-path incident matrix
				if (matrix_determinant(basic_matrix, num_link_abs) == 0)
				{
					break;
				}
			}
			else
			{
				break;
			}
		}
		while (1);
		// step 5.2 >> line search
		for (int fr = 0; fr < num_node_abs; fr++)
		{
			for (int to = 0; to < num_node_abs; to++)
			{
				left_od_flow_abs[fr][to] = pre_od_flow_abs[fr][to];
				right_od_flow_abs[fr][to] = od_flow_abs[fr][to];
			}
		}
		for (int l = 0; l < num_intervals; l++)
		{
			derivative = 0;
			for (int fr = 0; fr < num_node_abs; fr++)
			{
				for (int to = 0; to < num_node_abs; to++)
				{
					if (right_od_flow_abs[fr][to] + left_od_flow_abs[fr][to] > 0)
					{
						derivative += (right_od_flow_abs[fr][to] - left_od_flow_abs[fr][to]) * 
						log(0.5 * (right_od_flow_abs[fr][to] + left_od_flow_abs[fr][to]));
					}
					mid_od_flow_abs[fr][to] = 0.5 * (left_od_flow_abs[fr][to] + right_od_flow_abs[fr][to]);
				}
			}
			if (derivative < 0)
			{
				for (int fr = 0; fr < num_node_abs; fr++)
				{
					for (int to = 0; to < num_node_abs; to++)
					{
						left_od_flow_abs[fr][to] = mid_od_flow_abs[fr][to];
					}
				}
			}
			else if (derivative > 0)
			{
				for (int fr = 0; fr < num_node_abs; fr++)
				{
					for (int to = 0; to < num_node_abs; to++)
					{
						right_od_flow_abs[fr][to] = mid_od_flow_abs[fr][to];
					}
				}
			}
			else
			{
				break;
			}
		}
		// step 5.3 >> solution update
		for (int fr = 0; fr < num_node_abs; fr++)
		{
			for (int to = 0; to < num_node_abs; to++)
			{
				od_imp_abs[fr][to] = negative_infinity; // infinite O-D entropy impedances
				if (fr == to)
				{
					od_imp_abs[fr][to] = 0; // zero O-D entropy impedances
				}
			}
		}
		for (int fr = 0; fr < num_node_abs; fr++)
		{
			for (int to = 0; to < num_node_abs; to++)
			{
				od_flow_abs[fr][to] = mid_od_flow_abs[fr][to]; // update O-D flow rates
				if (od_flow_abs[fr][to] > 0)
				{
					od_imp_abs[fr][to] = log(od_flow_abs[fr][to]); // update O-D entropy impedances
				}
			}
		}
		cout << "\n";
		// step 5.4 >> convergence test
		accml_error = 0;
		for (int fr = 0; fr < num_node_abs; fr++)
		{
			for (int to = 0; to < num_node_abs; to++)
			{
				if (od_flow_abs[fr][to] > 0)
				{
					accml_error += abs(od_flow_abs[fr][to] - pre_od_flow_abs[fr][to]) / od_flow_abs[fr][to];
				}
			}
		}
		avg_error = accml_error / (pow((double) num_node_abs, 2) - (double) num_node_abs); // accumulated errors divided by number of O-D pairs
		cout << "   Average O-D flow difference: " << avg_error << "\n";
		if (avg_error < epsilon || iteration >= 25)
		{
			cout << "\n";
			break;
		}
	}
	while (1);

	// calculate the total trip rates over the network
	double total_trip_abs = 0;
	for (int fr = 0; fr < num_node_abs; fr++)
	{
		for (int to = 0; to < num_node_abs; to++)
		{
			total_trip_abs += abs(od_flow_abs[fr][to]);
		}
	}
	

	// step 6 >> write the O-D table into the output files
	
	if (!argv[4])
	{
		argv[4] = "default_od_abs.txt";
		endtime = clock() / CLOCKS_PER_SEC;
		cout << "\n";
		cout << "   Clock" << "\t\t" << endtime << " : Argument(s) missing...Default output file name used" << "\n";
	}
	ofstream odflow_out(argv[4]);

	odflow_out << "<NUMBER OF ZONES> " << node_id_abs[num_node_abs - 1] << "\n";
	odflow_out << "<TOTAL OD FLOW> " << total_trip_abs << "\n";
	odflow_out << "<END OF METADATA>" << "\n\n\n";
	for (int fr = 0; fr < num_node; fr++)
	{
		if (node_id_abs_in[fr] == true)
		{
			odflow_out << "Origin" << " " << "\t" << node_id[fr] << " \n";
			int count = 0;
			for (int to = 0; to < num_node; to++)
			{
				if (node_id_abs_in[to] == true)
				{
					odflow_out.width(5); odflow_out << right << node_id[to] << " :";
					odflow_out.width(9); odflow_out << right << setiosflags(ios::fixed) << setprecision(1) << abs(od_flow_abs[fr][to]) << ";";
					if ((count + 1) / 5 * 5 == (count + 1))
					{
						odflow_out << " \n";
					}
					count++;
				}
			}
			if (count / 5 * 5 != count + 1)
			{
				odflow_out << " \n";
			}
			odflow_out << "\n";
		}
	}

	odflow_out.close();

	endtime = clock() / CLOCKS_PER_SEC;
	cputime = endtime - starttime;
	cout << "   Clock" << "\t\t" << endtime << " : Finishing" << "\n";
	cout << "   CPU time elapsed: " << cputime << " sec" << "\n";

	return 0;

}

// calculate the inverse matrix by the Gauss-Jordan elimination
void matrix_inversion(double **basic, int size, double **inverse)
{
	double *row, coefficient, **original, **temp_original, **temp_inverse;
	row = new double [size];
	for (int i = 0; i < size; i++)
	{
		row[i] = 0;
	}
	// initialize the auxiliary matrix
	original = new double *[size];
	temp_original = new double *[size];
	temp_inverse = new double *[size];
	for (int i = 0; i < size; i++)
	{
		original[i] = new double [size];
		temp_original[i] = new double [size];
		temp_inverse[i] = new double [size];
		for (int j = 0; j < size; j++)
		{
			original[i][j] = basic[i][j];
			inverse[i][j] = 0;
		}
		inverse[i][i] = 1;
	}
	// perform the "scale and add" row operation
	for (int j = 0; j < size; j++)
	{
		for (int i = 0; i < size; i++)
		{
			if (original[i][j] != 0 && row[i] == 0)
			{
				for (int m = 0; m < size; m++)
				{
					if (m != i)
					{
						coefficient = original[m][j] / original[i][j];
						for (int n = 0; n < size; n++)
						{
							original[m][n] = original[m][n] - coefficient * original[i][n];
							inverse[m][n] = inverse[m][n] - coefficient * inverse[i][n];
						}
					}
				}
				row[i] = 1;
				break;
			}
		}
	}
	// perform the "swap" row operation
	for (int j = 0; j < size; j++)
	{
		for (int i = 0; i < size; i++)
		{
			if (original[i][j] != 0)
			{
				for (int n = 0; n < size; n++)
				{
					temp_original[j][n] = original[i][n];
					temp_inverse[j][n] = inverse[i][n];
				}
				break;
			}
		}
	}
	// perform the "scale" row operation
	for (int i = 0; i < size; i++)
	{
		coefficient = temp_original[i][i];
		temp_original[i][i] = 1;
		for (int j = 0; j < size; j++)
		{
			temp_inverse[i][j] = temp_inverse[i][j] / coefficient;
		}
	}
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			original[i][j] = temp_original[i][j];
			inverse[i][j] = temp_inverse[i][j];
		}
	}
	// release memory
	for (int i = 0; i < size; i++)
	{
		delete [] original[i];
		delete [] temp_original[i];
		delete [] temp_inverse[i];
	}
	delete [] original;
	delete [] temp_original;
	delete [] temp_inverse;
}

// calculate the matrix determinant by the Doolittle decomposition
double matrix_determinant(double **square, int size)
{
	double **lower, **upper, determinant;
	double positive_tiny = 0.000000000000001;
	// initialize the triangular matrices
	lower = new double *[size];
	upper = new double *[size];
	for (int i = 0; i < size; i++)
	{
		lower[i] = new double [size];
		upper[i] = new double [size];
		for (int j = 0; j < size; j++)
		{
			if (i == j)
			{
				lower[i][j] = 1;
			}
			else if (i < j)
			{
				lower[i][j] = 0;
			}
			else if (i > j)
			{
				upper[i][j] = 0;
			}

		}
	}
	// perform the row and column operations
	for (int i = 0; i < size; i++)
	{
		for (int j = i; j < size; j++)
		{
			upper[i][j] = square[i][j];
			for (int k = 0; k < i; k++)
			{
				upper[i][j] = upper[i][j] - lower[i][k] * upper[k][j];
			}
		}
		for (int j = i + 1; j < size; j++)
		{
			lower[j][i] = square[j][i];
			for (int k = 0; k < i; k++)
			{
				lower[j][i] = lower[j][i] - lower[j][k] * upper[k][i];
			}
			if (upper[i][i] != 0)
			{
				lower[j][i] = lower[j][i] / upper[i][i];
			}
			else
			{
				lower[j][i] = lower[j][i] / positive_tiny;
			}
		}
	}
	// calculate the determinant of the square matrix
	determinant = 1;
	for (int i = 0; i < size; i++)
	{
		determinant *= upper[i][i];
	}
	// release memory
	for (int i = 0; i < size; i++)
	{
		delete [] lower[i];
		delete [] lower[i];
	}
	delete [] lower;
	delete [] upper;
	// return the determinant value
	return determinant;
}

// detect the existence of negative loops by the Bellman-Ford algorithm
bool bellman_ford_shortest_path(int fr, int num_link, long int *link_address[][2], double *link_multiplier, 
								int num_node, long int *node_id, long int **predecessor_node_id, double *node_cost)
{
	int origin = fr;
	int fr_node_index, to_node_index;
	for (int n = 0; n < num_node; n++)
	{
		for (int m = 0; m < num_link; m++)
		{
			fr_node_index = link_address[m][0] - node_id;
			to_node_index = link_address[m][1] - node_id;
			if (node_cost[to_node_index] > node_cost[fr_node_index] + link_multiplier[m])
			{
				node_cost[to_node_index] = node_cost[fr_node_index] + link_multiplier[m];
				predecessor_node_id[to_node_index] = link_address[m][0];
			}
		}
	}
	bool loop_existence = false;
	for (int m = 0; m < num_link; m++)
	{
		fr_node_index = link_address[m][0] - node_id;
		to_node_index = link_address[m][1] - node_id;
		if (node_cost[to_node_index] > node_cost[fr_node_index] + link_multiplier[m])
		{
			loop_existence = true;
			break;
		}
	}
	return loop_existence;
}

// search for the all-to-all shortest paths by the Floyd-Warshall algorithm
void floyd_warshall_shortest_path(int num_node, long int *node_id, long int ***od_predecessor_node_id, double **od_cost)
{
	for (int n = 0; n < num_node; n++)
	{
		for (int n1 = 0; n1 < num_node; n1++)
		{
			for (int n2 = 0; n2 < num_node; n2++)
			{
				if (od_cost[n1][n2] > od_cost[n1][n] + od_cost[n][n2])
				{
					od_cost[n1][n2] = od_cost[n1][n] + od_cost[n][n2];
					od_predecessor_node_id[n1][n2] = od_predecessor_node_id[n][n2];
				}
			}
		}
	}
}