/* OUTPUT FILES MAY NOT BE REPORTING EVERYTHING YOU NEED.... CHECK THIS AND STILL GO WITH ODs AND CHECKING THE DEMAND FUNCTION */

//	Iterative Subnetwork Travel Demand Modeling Process
//
//	Chi Xie, 12-20-09.
//	Copyright (c) 2009 The University of Texas at Austin
//	$Revision: 1.7$  $Date: 2011/06/15       :  :   $
//	This program includes an iterative process of elastic O-D flow estimation,
//	mode split, time-of-day split, and traffic assignment, with the input files
//	in the format generated by Bar-Gera's OBA code

//	step -7 >> read the parameter file
//	step -6 >> read the original abstract network file
//	step -5 >> read the alternative abstract network file
//	step -4 >> read the abstract origin and destination demand growth file
//	step -3 >> read the abstract O-D demand file
//	step -2 >> generate random samples of input data and parameters
//	step -1 >> compute the initial generalized O-D cost
//	step 0 >> initialization
//	step 1-5 >> O-D demand estimation, mode split, time-of-day split and traffic assignment
//	step 6 >> write the result into the output files
//	step 7 >> write the statistical result into the output file
//	Note: steps -1 to 5 are repeated for the original and alternative networks

#include <iostream>
#include <fstream>
#include <sstream>
//#include <direct.h>
#include <iomanip>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <cmath>
#include <new>
using namespace std;

#define PI 3.14159265
#define MIN_FLOW_THRESHOLD 1e-5
#define MINIMUM_CAPACITY 1
#define END_OF_ROUTE -1
#define MAX_TRANSFERS 1

void floyd_warshall_shortest_path(int, long int *, long int ***, double **);
double lognormal_generator(double, double);
void readTransitFile(int ***routeMatrix, double **routeHeadway, int *numRoutes, char *transitFile, int numArcs);
void transitZoneZoneCosts(int num_node, long int *node_id, double **od_cost, int **routeMatrix, double *routeHeadway, int num_routes, long int **link_id, double *linkTravelTime);

int debugCounter;
char outputFileString[999];
//FILE *debugFile;
//FILE *myLogFile = fopen("mylog.txt", "w");
//FILE *myDebugFile = fopen("mydebuglog.txt", "w");

const int max_num_tod = 5; // maximum number of time-of-day periods
const int max_num_node = 100; // maximum number of nodes
const int max_num_link = 500; // maximum of links
const int max_num_route = 100000; // maximum of routes
const double negative_infinity = -2e9; // a very small negative constant
const double positive_infinity = 2e9; // a very lage positive constant
const double delta = 1e-6; // a very small positive constnat

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

//	fprintf(myLogFile, "-----------------\n|| STARTING NEW RUN\n--------------\n");

	cout << "\n";
	cout << "   Elastic Equilibrium Travel Demand Modeling Tool for Sketch Planning\n";
	cout << "   Designed and Coded by Chi Xie, Ph.D.\n";
	cout << "   Version 1.8, Updated on June 14, 2012\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 -7 >> read the parameter file

	char c_parameter_row[200];
	int num_tod; // number of times of day
	int num_mode; // number of transportation modes
	int num_class; // number of traveler classes
	double growth_factor; // demand growth factor over the planning period
	double mode_scale; // scale parameter for the incremental logit model for mode choice
	double tod_scale; // scale parameter for the incremental logit model for time-of-day choice
	char *pEnd;

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

	parameter_in.getline(c_parameter_row, 200);
	pEnd = c_parameter_row + 31;
	num_tod = strtol(pEnd, &pEnd, 10);
	parameter_in.getline(c_parameter_row, 200);
	pEnd = c_parameter_row + 24;
	num_class = strtol(pEnd, &pEnd, 10);
	parameter_in.getline(c_parameter_row, 200);
	pEnd = c_parameter_row + 27;
	num_mode = strtol(pEnd, &pEnd, 10);
	double *duration_t; // ratio of the duration of time-of-day period t over the whole day
	double elasticity, elasticity_input; // demand elasticity rate
	double *ecu_t; // equivalent car unit of time-of-day period t
	double *ecu_m; // equivalent car unit of transportation mode m
	double *oc_m; // operatin cost of transportation mode m
	double *pce_m; // passenger car equivalency of transportation mode m
	double *vot_c, *vot_c_input; // value of time of traveler class c
	double *vor_c, *vor_c_input; // value of reliability of traveler class c
	double *proportion_c; // proportion of traveler class c
	duration_t = new double [num_tod];
	ecu_t = new double [num_tod];
	ecu_m = new double [num_mode];
	oc_m = new double [num_mode];
	pce_m = new double [num_mode];
	vot_c = new double [num_class];
	vot_c_input = new double [num_class];
	vor_c = new double [num_class];
	vor_c_input = new double [num_class];
	proportion_c = new double [num_class];
	double **base_mode_prob;
	base_mode_prob = new double *[num_class];
	for (int c = 0; c < num_class; c++)
	{
		base_mode_prob[c] = new double [num_mode];
	}
	int max_outer_iteration_user; // maximum number of allowable outer iterations
	double cv_trip, cv_elasticity, cv_vot, cv_vor, cv_bpr; // coefficients of variation
	int num_scenarios; // number of Monte Carlo simulation samples
	int line = 3; // row index used in reading the parameter file
	while (!parameter_in.getline(c_parameter_row, 200).eof())
	{
		line++;
		if (line == 7)
		{
			pEnd = c_parameter_row + 12;
			growth_factor = strtod(pEnd, &pEnd);
		}
		else if (line >= 9 && line <= 13)
		{
			switch (line)
			{
			case 9:
				pEnd = c_parameter_row + 33;
				duration_t[0] = strtod(pEnd, &pEnd);
			case 10:
				pEnd = c_parameter_row + 33;
				duration_t[1] = strtod(pEnd, &pEnd);
			case 11:
				pEnd = c_parameter_row + 33;
				duration_t[2] = strtod(pEnd, &pEnd);
			case 12:
				pEnd = c_parameter_row + 33;
				duration_t[3] = strtod(pEnd, &pEnd);
			case 13:
				pEnd = c_parameter_row + 33;
				duration_t[4] = strtod(pEnd, &pEnd);
			}
		}
		else if (line >= 15 && line <= 19)
		{
			pEnd = c_parameter_row + 42;
			elasticity_input = strtod(pEnd, &pEnd);
		}
		else if (line >= 21 && line <= 25)
		{
			switch (line)
			{
			case 21:
				pEnd = c_parameter_row + 30;
				ecu_m[0] = strtod(pEnd, &pEnd);
			case 22:
				pEnd = c_parameter_row + 30;
				ecu_m[1] = strtod(pEnd, &pEnd);
			case 23:
				pEnd = c_parameter_row + 30;
				ecu_m[2] = strtod(pEnd, &pEnd);
			case 24:
				pEnd = c_parameter_row + 30;
				ecu_m[3] = strtod(pEnd, &pEnd);
			case 25:
				pEnd = c_parameter_row + 30;
				ecu_m[4] = strtod(pEnd, &pEnd);
			}
		}
		else if (line >= 27 && line <= 31)
		{
			switch (line)
			{
			case 27:
				pEnd = c_parameter_row + 23;
				oc_m[0] = strtod(pEnd, &pEnd);
			case 28:
				pEnd = c_parameter_row + 23;
				oc_m[1] = strtod(pEnd, &pEnd);
			case 29:
				pEnd = c_parameter_row + 23;
				oc_m[2] = strtod(pEnd, &pEnd);
			case 30:
				pEnd = c_parameter_row + 23;
				oc_m[3] = strtod(pEnd, &pEnd);
			case 31:
				pEnd = c_parameter_row + 23;
				oc_m[4] = strtod(pEnd, &pEnd);
			}
		}
		else if (line >= 33 && line <= 37)
		{
			switch (line)
			{
			case 33:
				pEnd = c_parameter_row + 30;
				pce_m[0] = strtod(pEnd, &pEnd);
			case 34:
				pEnd = c_parameter_row + 30;
				pce_m[1] = strtod(pEnd, &pEnd);
			case 35:
				pEnd = c_parameter_row + 30;
				pce_m[2] = strtod(pEnd, &pEnd);
			case 36:
				pEnd = c_parameter_row + 30;
				pce_m[3] = strtod(pEnd, &pEnd);
			case 37:
				pEnd = c_parameter_row + 30;
				pce_m[4] = strtod(pEnd, &pEnd);
			}
		}
		else if (line >= 39 && line <= 42)
		{
			switch (line)
			{
			case 39:
				pEnd = c_parameter_row + 26;
				vot_c_input[0] = strtod(pEnd, &pEnd);
			case 40:
				pEnd = c_parameter_row + 26;
				vot_c_input[1] = strtod(pEnd, &pEnd);
			case 41:
				pEnd = c_parameter_row + 26;
				vot_c_input[2] = strtod(pEnd, &pEnd);
			case 42:
				pEnd = c_parameter_row + 26;
				vot_c_input[3] = strtod(pEnd, &pEnd);
			}
		}
		else if (line >= 44 && line <= 47)
		{
			switch (line)
			{
			case 44:
				pEnd = c_parameter_row + 26;
				vor_c_input[0] = strtod(pEnd, &pEnd);
			case 45:
				pEnd = c_parameter_row + 26;
				vor_c_input[1] = strtod(pEnd, &pEnd);
			case 46:
				pEnd = c_parameter_row + 26;
				vor_c_input[2] = strtod(pEnd, &pEnd);
			case 47:
				pEnd = c_parameter_row + 26;
				vor_c_input[3] = strtod(pEnd, &pEnd);
			}
		}
		else if (line >= 49 && line <= 52)
		{
			switch (line)
			{
			case 49:
				pEnd = c_parameter_row + 27;
				proportion_c[0] = strtod(pEnd, &pEnd);
			case 50:
				pEnd = c_parameter_row + 27;
				proportion_c[1] = strtod(pEnd, &pEnd);
			case 51:
				pEnd = c_parameter_row + 27;
				proportion_c[2] = strtod(pEnd, &pEnd);
			case 52:
				pEnd = c_parameter_row + 27;
				proportion_c[3] = strtod(pEnd, &pEnd);
			}
		}
		else if (line == 54)
		{
			pEnd = c_parameter_row + 27;
			mode_scale = strtod(pEnd, &pEnd);
		}
		else if (line >= 56 && line <= 59)
		{
			switch (line)
			{
			case 56:
				pEnd = c_parameter_row + 37;
				base_mode_prob[0][0] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[0][1] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[0][2] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[0][3] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[0][4] = strtod(pEnd, &pEnd);
			case 57:
				pEnd = c_parameter_row + 37;
				base_mode_prob[1][0] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[1][1] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[1][2] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[1][3] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[1][4] = strtod(pEnd, &pEnd);
			case 58:
				pEnd = c_parameter_row + 37;
				base_mode_prob[2][0] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[2][1] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[2][2] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[2][3] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[2][4] = strtod(pEnd, &pEnd);
			case 59:
				pEnd = c_parameter_row + 37;
				base_mode_prob[3][0] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[3][1] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[3][2] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[3][3] = strtod(pEnd, &pEnd);
				pEnd = pEnd + 1;
				base_mode_prob[3][4] = strtod(pEnd, &pEnd);
			}
		}
		else if (line == 61)
		{
			pEnd = c_parameter_row + 34;
			tod_scale = strtod(pEnd, &pEnd);
		}
		else if (line == 63)
		{
			pEnd = c_parameter_row + 25;
			max_outer_iteration_user = strtod(pEnd, &pEnd);
		}
		else if (line == 65)
		{
			pEnd = c_parameter_row + 17;
			cv_trip = strtod(pEnd, &pEnd);
		}
		else if (line == 66)
		{
			pEnd = c_parameter_row + 24;
			cv_elasticity = strtod(pEnd, &pEnd);
		}
		else if (line == 67)
		{
			pEnd = c_parameter_row + 10;
			cv_vot = strtod(pEnd, &pEnd);
		}
		else if (line == 68)
		{
			pEnd = c_parameter_row + 10;
			cv_vor = strtod(pEnd, &pEnd);
		}
		else if (line == 69)
		{
			pEnd = c_parameter_row + 21;
			cv_bpr = strtod(pEnd, &pEnd);
		}
		else if (line == 71)
		{
			pEnd = c_parameter_row + 20;
			num_scenarios = strtol(pEnd, &pEnd, 10);
		}
	}
	parameter_in.close();
	//double w = vor_c_input[0] / vot_c_input[0]; // ratio of value of reliability over value of mean (travel time)

	// step -6 >> read the original abstract network file
//	fprintf(myLogFile, "step -6\n");
	char metadata_old[200]; // temporary character variable
	char header_old[200]; // temporary character variable
	char c_link_row_old[500]; // temporary character variable
	int num_link_old; // number of links in the original network file
	long int link_id_old[max_num_link][2]; // IDs of links in the orignal network file
	long int *link_address_old[max_num_link][2]; // addresses of links in the original network file
	double link_attrbt_old_t[max_num_tod][max_num_link][16]; // attributes of links in the original file
	int k;

	for (int t = 0; t < num_tod; t++)
	{
		if (!argv[t + 2])
		{
			cout << "   argument(s) missing...";
			return 1;
		}
		ifstream network_in_old(argv[t + 2]);
		for (int i = 0; i < 7; i++)
		{
			network_in_old.getline(metadata_old, 50);
		}
		network_in_old.getline(header_old, 200);
		k = 0;
		while(!network_in_old.getline(c_link_row_old, 500).eof())
		{
			char *pEnd;
			link_id_old[k][0] = strtol(c_link_row_old, &pEnd, 10); // initial id
			link_id_old[k][1] = strtol(pEnd, &pEnd, 10); // terminal id
			link_attrbt_old_t[t][k][0] = strtod(pEnd, &pEnd); // capacity
			if (link_attrbt_old_t[t][k][0] <= 0) link_attrbt_old_t[t][k][0] = MINIMUM_CAPACITY;
			link_attrbt_old_t[t][k][1] = strtod(pEnd, &pEnd); // length
			link_attrbt_old_t[t][k][2] = strtod(pEnd, &pEnd) / 60; // free flow time
			link_attrbt_old_t[t][k][3] = strtod(pEnd, &pEnd); // base (alpha)
			link_attrbt_old_t[t][k][4] = strtod(pEnd, &pEnd); // power (beta)
			link_attrbt_old_t[t][k][5] = strtod(pEnd, &pEnd); // speed limit
			link_attrbt_old_t[t][k][6] = strtod(pEnd, &pEnd); // toll (mode 1)
			link_attrbt_old_t[t][k][7] = strtod(pEnd, &pEnd); // toll (mode 2)
			link_attrbt_old_t[t][k][8] = strtod(pEnd, &pEnd); // toll (mode 3)
			link_attrbt_old_t[t][k][9] = strtod(pEnd, &pEnd); // toll (mode 4)
			link_attrbt_old_t[t][k][10] = strtod(pEnd, &pEnd); // toll (mode 5)
			link_attrbt_old_t[t][k][11] = strtod(pEnd, &pEnd); // type
			link_attrbt_old_t[t][k][12] = strtod(pEnd, &pEnd); // free flow variance
			link_attrbt_old_t[t][k][13] = strtod(pEnd, &pEnd); // (sigma)
			link_attrbt_old_t[t][k][14] = strtod(pEnd, &pEnd); // (gamma)
			link_attrbt_old_t[t][k][15] = strtod(pEnd, &pEnd); // (tau)
			k++;
		}
		num_link_old = k;
		network_in_old.close();
	}

	long int node_id_old[max_num_node];
	int p = 0;
	bool node_existence_old = false;
	for (int l = 0; l < num_link_old; l++)
	{
		for (int n = 0; n < p; n++)
		{
			if (link_id_old[l][0] == node_id_old[n])
			{
				node_existence_old = true;
				break;
			}
		}
		if (node_existence_old == false)
		{
			node_id_old[p] = link_id_old[l][0];
			p++;
		}
		node_existence_old = false;
	}
	for (int l = 0; l < num_link_old; l++)
	{
		for (int n = 0; n < p; n++)
		{
			if (link_id_old[l][1] == node_id_old[n])
			{
				node_existence_old = true;
				break;
			}
		}
		if (node_existence_old == false)
		{
			node_id_old[p] = link_id_old[l][1];
			p++;
		}
		node_existence_old = false;
	}
	int num_node_old = p;
	for (int l = 0; l < num_link_old; l++)
	{
		for (int n = 0; n < num_node_old; n++)
		{
			if (link_id_old[l][0] == node_id_old[n])
			{
				link_address_old[l][0] = &node_id_old[n];
			}
			if (link_id_old[l][1] == node_id_old[n])
			{
				link_address_old[l][1] = &node_id_old[n];
			}
		}
	}


	// step -5 >> read the alternative abstract network file
//	fprintf(myLogFile, "step -5\n");
	char metadata[200]; // temporary character variable
	char header[200]; // temporary character variable
	char c_link_row[500]; // temporary character variable
	int num_link; // number of links in the alternative network
	long int **link_id = (long int **)malloc(max_num_link * sizeof(long int *)); for (int l = 0; l < max_num_link; l++) link_id[l] = (long int *)malloc(2 * sizeof(long int));
	//long int link_id[max_num_link][2]; // ID of links in the alternative network
	long int *link_address[max_num_link][2]; // address of links in the alternative network
	double link_attrbt_t[max_num_tod][max_num_link][16], link_attrbt_input_t[max_num_tod][max_num_link][16]; // attributes of links in the alternative network


	for (int t = 0; t < num_tod; t++)
	{
		if (!argv[t + 2 + num_tod])
		{
			cout << "   argument(s) missing...";
			return 1;
		}
		ifstream network_in(argv[t + 2 + num_tod]);
		for (int i = 0; i < 7; i++)
		{
			network_in.getline(metadata, 50);
		}
		network_in.getline(header, 200);
		k = 0;
		while(!network_in.getline(c_link_row, 500).eof())
		{
			if (strlen(c_link_row) == 0) break;
			char *pEnd;
			link_id[k][0] = strtol(c_link_row, &pEnd, 10); // initial id
			link_id[k][1] = strtol(pEnd, &pEnd, 10); // terminal id
			link_attrbt_input_t[t][k][0] = strtod(pEnd, &pEnd); // capacity
			if (link_attrbt_input_t[t][k][0] <= 0) link_attrbt_input_t[t][k][0] = MINIMUM_CAPACITY;
			link_attrbt_input_t[t][k][1] = strtod(pEnd, &pEnd); // length
			link_attrbt_input_t[t][k][2] = strtod(pEnd, &pEnd) / 60; // free flow time
			link_attrbt_input_t[t][k][3] = strtod(pEnd, &pEnd); // base (alpha)
			link_attrbt_input_t[t][k][4] = strtod(pEnd, &pEnd); // power (beta)
			link_attrbt_input_t[t][k][5] = strtod(pEnd, &pEnd); // speed limit
			link_attrbt_input_t[t][k][6] = strtod(pEnd, &pEnd); // toll
			link_attrbt_input_t[t][k][7] = strtod(pEnd, &pEnd); // type
			link_attrbt_input_t[t][k][8] = strtod(pEnd, &pEnd); // free flow variance
			link_attrbt_input_t[t][k][9] = strtod(pEnd, &pEnd); // (sigma)
			link_attrbt_input_t[t][k][10] = strtod(pEnd, &pEnd); // (gamma)
			link_attrbt_input_t[t][k][11] = strtod(pEnd, &pEnd); // (tau)
			k++;
		}
		num_link = k; // number of links
		network_in.close();
	}

	long int node_id[max_num_node];

	p = 0;

	bool node_existence = false;
	for (int l = 0; l < num_link; l++)
	{
		for (int n = 0; n < p; n++)
		{
			if (link_id[l][0] == node_id[n])
			{
				node_existence = true;
				break;
			}
		}
		if (node_existence == false)
		{
			node_id[p] = link_id[l][0];
			p++;
		}
		node_existence = false;
	}
	for (int l = 0; l < num_link; l++)
	{
		for (int n = 0; n < p; n++)
		{
			if (link_id[l][1] == node_id[n])
			{
				node_existence = true;
				break;
			}
		}
		if (node_existence == false)
		{
			node_id[p] = link_id[l][1];
			p++;
		}
		node_existence = false;
	}
	int num_node = p; // number of nodes
	for (int l = 0; l < num_link; l++)
	{
		for (int n = 0; n < num_node; n++)
		{
			if (link_id[l][0] == node_id[n])
			{
				link_address[l][0] = &node_id[n];
			}
			if (link_id[l][1] == node_id[n])
			{
				link_address[l][1] = &node_id[n];
			}
		}
	}

	// step -3.5 >> read transit route files
//	fprintf(myLogFile, "step -3.5\n"); fflush(myLogFile);
	char transitFileName[999];
	long int r;
	int totalroutes = 0;
	int ***baseRouteMatrix_t = (int ***)malloc(num_tod * sizeof(int **));
	double **baseHeadway_t = (double **)malloc(num_tod * sizeof(double *));
	int *numBaseRoutes_t = (int *)malloc(num_tod * sizeof(int));
	for (int t = 0; t < num_tod; t++) {
		strcpy(transitFileName, "Transit_");
		strcat(transitFileName, argv[t+2]);
		readTransitFile(&baseRouteMatrix_t[t], &baseHeadway_t[t], &numBaseRoutes_t[t], transitFileName, num_link);
		totalroutes += numBaseRoutes_t[t];
	}
	int ***alternateRouteMatrix_t = (int ***)malloc(num_tod * sizeof(int **));
	double **alternateHeadway_t = (double **)malloc(num_tod * sizeof(double *));
	int *numAlternateRoutes_t = (int *)malloc(num_tod * sizeof(int));
	for (int t = 0; t < num_tod; t++) {
		strcpy(transitFileName, "Transit_");
		strcat(transitFileName, argv[t+7]);
		readTransitFile(&alternateRouteMatrix_t[t], &alternateHeadway_t[t], &numAlternateRoutes_t[t], transitFileName, num_link);
		totalroutes += numAlternateRoutes_t[t];
	}
	printf("\t%d total transit routes\n", totalroutes);
	//if (totalroutes == 0) num_mode--;

	double ****baseRouteTravelTimeMatrix_t = (double ****)malloc(num_tod * sizeof(double ***));
	for (int t = 0; t < num_tod; t++) {
		baseRouteTravelTimeMatrix_t[t] = (double ***)malloc(numBaseRoutes_t[t] * sizeof(double **));
		for (int r = 0; r < numBaseRoutes_t[r]; r++) {
			baseRouteTravelTimeMatrix_t[t][r] = (double **)malloc(num_node * sizeof(double *));
			for (int n1 = 0; n1 < num_node; n1++) {
				baseRouteTravelTimeMatrix_t[t][r][n1] = (double *)malloc(num_node * sizeof(double));
			}
		}
	}
	double ****alternateRouteTravelTimeMatrix_t = (double ****)malloc(num_tod * sizeof(double ***));
	for (int t = 0; t < num_tod; t++) {
		alternateRouteTravelTimeMatrix_t[t] = (double ***)malloc(numAlternateRoutes_t[t] * sizeof(double **));
		for (int r = 0; r < numAlternateRoutes_t[r]; r++) {
			alternateRouteTravelTimeMatrix_t[t][r] = (double **)malloc(num_node * sizeof(double *));
			for (int n1 = 0; n1 < num_node; n1++) {
				alternateRouteTravelTimeMatrix_t[t][r][n1] = (double *)malloc(num_node * sizeof(double));
			}
		}
	}
	double ****baseTransitTravelTime_t = (double ****)malloc(num_tod * sizeof(double ***));
	for (int t = 0; t < num_tod; t++) {
		baseTransitTravelTime_t[t] = (double ***)malloc(num_node * sizeof(double **));
		for (int n1 = 0; n1 < num_node; n1++) {
			baseTransitTravelTime_t[t][n1] = (double **)malloc(num_node * sizeof(double *));
			for (int n2 = 0; n2 < num_node; n2++) {
				baseTransitTravelTime_t[t][n1][n2] = (double *)malloc((MAX_TRANSFERS + 1) * sizeof(double));
			}
		}
	}
	double ****alternateTransitTravelTime_t = (double ****)malloc(num_tod * sizeof(double ***));
	for (int t = 0; t < num_tod; t++) {
		alternateTransitTravelTime_t[t] = (double ***)malloc(num_node * sizeof(double **));
		for (int n1 = 0; n1 < num_node; n1++) {
			alternateTransitTravelTime_t[t][n1] = (double **)malloc(num_node * sizeof(double *));
			for (int n2 = 0; n2 < num_node; n2++) {
				alternateTransitTravelTime_t[t][n1][n2] = (double *)malloc((MAX_TRANSFERS + 1) * sizeof(double));
			}

		}
	}

	// step -4 >> read the abstract origin and destination demand growth file and fixed costs
//	fprintf(myLogFile, "step -4\n"); fflush(myLogFile);
	if (!argv[2 + num_tod + num_tod])
	{
		cout << "   argument(s) missing...";
		return 1;
	}
	ifstream growth_in(argv[2 + num_tod + num_tod]);

	double growth_list[max_num_node][2]; // list of growth factors
	double **growth_table; // O-D table of growth factors
	growth_table = new double *[num_node];
	for (int n = 0; n < num_node; n++)
	{
		growth_table[n] = new double [num_node];
	}
	char c_growth_row[200];

	for (int i = 0; i < 6; i++)
	{
		growth_in.getline(metadata, 50);
	}
	growth_in.getline(header, 200);
	int temp_id;
	p = 0;
	while(!growth_in.getline(c_growth_row, 200).eof())
	{
		char *pEnd;
		temp_id = strtol(c_growth_row, &pEnd, 10); // temporary node IDs
		growth_list[p][0] = strtod(pEnd, &pEnd); // origin growth factors
		growth_list[p][1] = strtod(pEnd, &pEnd); // destination growth factors
		p++;
	}
	growth_in.close();
	for (int n1 = 0; n1 < num_node; n1++)
	{
		for (int n2 = 0; n2 < num_node; n2++)
		{
			growth_table[n1][n2] = 0.5 * (growth_list[n1][0] + growth_list[n2][1]);
		}
	}

	if (!argv[3 + num_tod + num_tod])
	{
		cout << "   argument(s) missing...";
		return 1;
	}
	//ifstream fixed_in(argv[3 + num_tod + num_tod]);

	// Read base fixed costs
//	fprintf(myLogFile, "Reading fixed costs.\n"); fflush(myLogFile);
	ifstream fixed_in(argv[3 + num_tod + num_tod]);

	double ****fixed_cost = (double ****)malloc(max_num_node * sizeof(double ***));
	double ****fixed_cost_base = (double ****)malloc(max_num_node * sizeof(double ***));
	double ****fixed_cost_alt = (double ****)malloc(max_num_node * sizeof(double ***));
	for (int n1 = 0; n1 < max_num_node; n1++) {
		fixed_cost[n1] = (double ***)malloc(num_tod * sizeof(double **));
		fixed_cost_base[n1] = (double ***)malloc(num_tod * sizeof(double **));
		fixed_cost_alt[n1] = (double ***)malloc(num_tod * sizeof(double **));
		for (int t = 0; t < num_tod; t++) {
			fixed_cost[n1][t] = (double **)malloc(num_mode * sizeof(double *));
			fixed_cost_base[n1][t] = (double **)malloc(num_mode * sizeof(double *));
			fixed_cost_alt[n1][t] = (double **)malloc(num_mode * sizeof(double *));
			for (int m = 0; m < num_mode; m++) {
				fixed_cost[n1][t][m] = (double *)malloc(2 * sizeof(double));
				fixed_cost_base[n1][t][m] = (double *)malloc(2 * sizeof(double));
				fixed_cost_alt[n1][t][m] = (double *)malloc(2 * sizeof(double));
				for (int od = 0; od <= 1; od++) {
					fixed_cost[n1][t][m][od] = 0;
					fixed_cost_base[n1][t][m][od] = 0;
					fixed_cost_alt[n1][t][m][od] = 0;
				}
			}
		}
	}

	/*
	double ****fixed_cost, ****fixed_cost_base;
	fixed_cost = new double ***[max_num_node];
	fixed_cost_base = new double ***[max_num_node];
	for (int n1 = 0; n1 < max_num_node; n1++) {
		fixed_cost[n1] = new double **[num_tod];
		fixed_cost_base[n1] = new double **[num_tod];
		for (int t = 0; t < num_tod; t++) {
			fixed_cost[n1][t] = new double *[num_mode];
			fixed_cost_base[n1][t] = new double *[num_mode];
			for (int m = 0; m < num_mode; m++) {
				fixed_cost[n1][t][m] = new double [2];
				fixed_cost_base[n1][t][m] = new double [2];
			}
		}
	}
	*/

	char c_fixed_row[9999];
	// Four header rows...
	for (int ctr = 0; ctr < 4; ctr++) {
		fixed_in.getline(c_fixed_row, 999);
	}
	// ...then data
	while (!fixed_in.getline(c_fixed_row, 9990).eof()) {
		if (strlen(c_fixed_row) == 0) break;
		char *pEnd;
		p = strtol(c_fixed_row, &pEnd, 10); // temporary node IDs
		for (int od = 0; od <= 1; od++) {
			for (int t = 0; t < num_tod; t++) {
				for (int m = 0; m < num_mode; m++) {
               fixed_cost_base[p][t][m][od] = strtod(pEnd, &pEnd); // destination growth factors
				}
			}
		}
	}
	fixed_in.close();

	// Read alternate fixed costs
	fixed_in.open(argv[4 + num_tod + num_tod]);
	//double fixed_cost_alt[max_num_node][num_tod][num_mode][2]; // list of fixed costs... [0] is for origin and [1] is for destination
	// Four header rows...
	for (int ctr = 0; ctr < 4; ctr++) {
		fixed_in.getline(c_fixed_row, 999);
	}
	// ...then data
	while (!fixed_in.getline(c_fixed_row, 9990).eof()) {
		if (strlen(c_fixed_row) == 0) break;
		//sscanf(c_fixed_row, " %d", &p);
		char *pEnd;
		p = strtol(c_fixed_row, &pEnd, 10); // temporary node IDs
		for (int od = 0; od <= 1; od++) {
			for (int t = 0; t < num_tod; t++) {
				for (int m = 0; m < num_mode; m++) {
               fixed_cost_alt[p][t][m][od] = strtod(pEnd, &pEnd); // destination growth factors
				}
			}
		}
	}
	fixed_in.close();

	// step -3 >> read the abstract O-D demand file
//	fprintf(myLogFile, "step -3\n"); fflush(myLogFile);
	char c_od_row[200]; // temporary character variable
	double ***base_od_flow_t, ***base_od_flow_t_input; // initial O-D flow rate by time of day
	base_od_flow_t = new double **[num_tod];
	base_od_flow_t_input = new double **[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		base_od_flow_t[t] = new double *[num_node];
		base_od_flow_t_input[t] = new double *[num_node];
		for (int n1 = 0; n1 < num_node; n1++)
		{
			base_od_flow_t[t][n1] = new double [num_node];
			base_od_flow_t_input[t][n1] = new double [num_node];
			for (int n2 = 0; n2 < num_node; n2++)
			{
				base_od_flow_t[t][n1][n2] = 0;
				base_od_flow_t_input[t][n1][n2] = 0;
			}
		}
	}
	long int s;

	for (int t = 0; t < num_tod; t++)
	{
		if (!argv[t + 5 + num_tod + num_tod])
		{
			cout << "   argument(s) missing...";
			return 1;
		}
//		fprintf(myLogFile, "reading %s\n", argv[t + 5 + num_tod + num_tod]); fflush(myLogFile);
		ifstream demand_in(argv[t + 5 + num_tod + num_tod]);
		for (int i = 0; i < 5; i++)
		{
			demand_in.getline(metadata, 50);
		}
		while (!demand_in.getline(c_od_row, 200).eof())
		{
			char *pEnd;
			pEnd = c_od_row + 8;
			//r = strtod(pEnd, &pEnd);
			sscanf(c_od_row, "Origin %ld", &r);
			do
			{
				demand_in.getline(c_od_row, 200);
				if (strlen(c_od_row) > 2)
				{
					pEnd = c_od_row;
					for (int num_col = 0; num_col < 5; num_col++)
					{
						s = strtol(pEnd, &pEnd, 10);
						if (s > 0)
						{
							pEnd = pEnd + 2;
							base_od_flow_t_input[t][r - 1][s - 1] = strtod(pEnd, &pEnd) * (1 + growth_table[r - 1][s - 1]);
							pEnd = pEnd + 1;
						}
						else
						{
							break;
						}
					}
				}
				else
				{
					break;
				}
			}
			while (1);
		}
		demand_in.close();
//		fprintf(myLogFile, "done.\n");
		for (r = 0; r < num_node; r++) base_od_flow_t_input[t][r][r] = 0;
	}


	// create a subfolder for individual scenario results
//	fprintf(myLogFile, "creating subfolder\n"); fflush(myLogFile);
	char current_directory[100], original_directory[100];
//	getcwd(current_directory, 100);
	for (int c = 0; c < 100; c++)
	{
		original_directory[c] = current_directory[c];
	}
	if (num_scenarios > 1)
	{
		strcat(current_directory, "/Individual_Output");
//		mkdir(current_directory);
//		chdir(current_directory);
	}
//	fprintf(myLogFile, "done\n"); fflush(myLogFile);
	double *user_surplus_series; user_surplus_series = new double [num_scenarios]; // set of sample user surplus
	double *travel_cost_series; travel_cost_series = new double [num_scenarios]; // set of sample travel costs
	double *travel_time_series; travel_time_series = new double [num_scenarios]; // set of sample travel times
	double *vht_series; vht_series = new double [num_scenarios]; // set of vehicle hours traveled
	double *vmt_series; vmt_series = new double [num_scenarios]; // set of vehicle miles traveled
	double **total_trip_rate_series; total_trip_rate_series = new double* [num_class];
	for (int c = 0; c < num_class; c++)
	{
		total_trip_rate_series[c] = new double [num_tod];
	}

	for (int scenario = 0; scenario < num_scenarios; scenario++)
	{

	// step -2 >> generate random samples of input data and parameters
//	fprintf(myLogFile, "step -2\n"); fflush(myLogFile);
	if (scenario < num_scenarios - 1)
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					base_od_flow_t[t][n1][n2] = lognormal_generator(base_od_flow_t_input[t][n1][n2], cv_trip * base_od_flow_t_input[t][n1][n2]);
				}
			}
		}
		elasticity = (-1) * lognormal_generator((-1) * elasticity_input, cv_elasticity * (-1) * elasticity_input);
		for (int c = 0; c < num_class; c++)
		{
			vot_c[c] = lognormal_generator(vot_c_input[c], cv_vot * vot_c_input[c]);
			vor_c[c] = lognormal_generator(vor_c_input[c], cv_vor * vor_c_input[c]);
		}
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_attrbt_t[t][l][0] = link_attrbt_input_t[t][l][0]; // capacity
				link_attrbt_t[t][l][1] = link_attrbt_input_t[t][l][1]; // length
				link_attrbt_t[t][l][2] = link_attrbt_input_t[t][l][2]; // free flow time
				link_attrbt_t[t][l][3] = link_attrbt_input_t[t][l][3]; // base (alpha)
				link_attrbt_t[t][l][4] = link_attrbt_input_t[t][l][4]; // power (beta)
				//link_attrbt_t[t][l][3] = lognormal_generator(link_attrbt_input_t[t][l][3], cv_bpr * link_attrbt_input_t[t][l][3]); // base (alpha)
				//link_attrbt_t[t][l][4] = lognormal_generator(link_attrbt_input_t[t][l][4], cv_bpr * link_attrbt_input_t[t][l][4]); // power (beta)
				link_attrbt_t[t][l][5] = link_attrbt_input_t[t][l][5]; // speed limit
				link_attrbt_t[t][l][6] = link_attrbt_input_t[t][l][6]; // toll
				link_attrbt_t[t][l][7] = link_attrbt_input_t[t][l][7]; // toll
				link_attrbt_t[t][l][8] = link_attrbt_input_t[t][l][8]; // toll
				link_attrbt_t[t][l][9] = link_attrbt_input_t[t][l][9]; // toll
				link_attrbt_t[t][l][10] = link_attrbt_input_t[t][l][10]; // toll
				link_attrbt_t[t][l][11] = link_attrbt_input_t[t][l][11]; // type
				link_attrbt_t[t][l][12] = link_attrbt_input_t[t][l][12]; // free flow variance
				link_attrbt_t[t][l][13] = link_attrbt_input_t[t][l][13]; // (sigma)
				link_attrbt_t[t][l][14] = link_attrbt_input_t[t][l][14]; // (gamma)
				link_attrbt_t[t][l][15] = link_attrbt_input_t[t][l][15]; // (tau)
				//link_attrbt_t[t][l][13] = lognormal_generator(link_attrbt_input_t[t][l][13], cv_bpr * link_attrbt_input_t[t][l][13]); // (sigma)
				//link_attrbt_t[t][l][14] = lognormal_generator(link_attrbt_input_t[t][l][14], cv_bpr * link_attrbt_input_t[t][l][14]); // (gamma)
				//link_attrbt_t[t][l][15] = lognormal_generator(link_attrbt_input_t[t][l][15], cv_bpr * link_attrbt_input_t[t][l][15]); // (tau)
			}
		}
	}
	else // the last sample uses the input data without randomness
	{

		for (int t = 0; t < num_tod; t++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					base_od_flow_t[t][n1][n2] = base_od_flow_t_input[t][n1][n2];
				}
			}
		}
		elasticity = elasticity_input;
		for (int c = 0; c < num_class; c++)
		{
			vot_c[c] = vot_c_input[c];
			vor_c[c] = vor_c_input[c];
		}
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_attrbt_t[t][l][0] = link_attrbt_input_t[t][l][0]; // capacity
				link_attrbt_t[t][l][1] = link_attrbt_input_t[t][l][1]; // length
				link_attrbt_t[t][l][2] = link_attrbt_input_t[t][l][2]; // free flow time
				link_attrbt_t[t][l][3] = link_attrbt_input_t[t][l][3]; // base (alpha)
				link_attrbt_t[t][l][4] = link_attrbt_input_t[t][l][4]; // power (beta)
				link_attrbt_t[t][l][5] = link_attrbt_input_t[t][l][5]; // speed limit
				link_attrbt_t[t][l][6] = link_attrbt_input_t[t][l][6]; // toll
				link_attrbt_t[t][l][7] = link_attrbt_input_t[t][l][7]; // toll
				link_attrbt_t[t][l][8] = link_attrbt_input_t[t][l][8]; // toll
				link_attrbt_t[t][l][9] = link_attrbt_input_t[t][l][9]; // toll
				link_attrbt_t[t][l][10] = link_attrbt_input_t[t][l][10]; // toll
				link_attrbt_t[t][l][11] = link_attrbt_input_t[t][l][11]; // type
				link_attrbt_t[t][l][12] = link_attrbt_input_t[t][l][12]; // free flow variance
				link_attrbt_t[t][l][13] = link_attrbt_input_t[t][l][13]; // (sigma)
				link_attrbt_t[t][l][14] = link_attrbt_input_t[t][l][14]; // (gamma)
				link_attrbt_t[t][l][15] = link_attrbt_input_t[t][l][15]; // (tau)
			}
		}
	}

	double **od_cost;
	od_cost = new double *[num_node];
	for (int n = 0; n < num_node; n++)
	{
		od_cost[n] = new double [num_node];
	}
	//
	double *****base_od_flow_t_c_m, *****base_od_cost_t_c_m; // initial O-D flow rates by time of day, user class, transportation mode
	base_od_flow_t_c_m = new double ****[num_tod];
	base_od_cost_t_c_m = new double ****[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		base_od_flow_t_c_m[t] = new double ***[num_class];
		base_od_cost_t_c_m[t] = new double ***[num_class];
		for (int c = 0; c < num_class; c++)
		{
			base_od_flow_t_c_m[t][c] = new double **[num_mode];
			base_od_cost_t_c_m[t][c] = new double **[num_mode];
			for (int m = 0; m < num_mode; m++)
			{
				base_od_flow_t_c_m[t][c][m] = new double *[num_node];
				base_od_cost_t_c_m[t][c][m] = new double *[num_node];
				for (int n1 = 0; n1 < num_node; n1++)
				{
					base_od_flow_t_c_m[t][c][m][n1] = new double [num_node];
					base_od_cost_t_c_m[t][c][m][n1] = new double [num_node];
				}
			}
		}
	}
	//
	double *****base_tod_prob; // initial time-of-day choice probability
	base_tod_prob = new double ****[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		base_tod_prob[t] = new double ***[num_class];
		for (int c = 0; c < num_class; c++)
		{
			base_tod_prob[t][c] = new double **[num_mode];
			for (int m = 0; m < num_mode; m++)
			{
				base_tod_prob[t][c][m] = new double *[num_node];
				for (int n1 = 0; n1 < num_node; n1++)
				{
					base_tod_prob[t][c][m][n1] = new double [num_node];
				}
			}
		}
	}
	double base_od_flow; // initial O-D flow rates
	//
	double **link_flow_t, **pre_link_flow_t, **left_link_flow_t, **right_link_flow_t, **mid_link_flow_t, **link_flow_t_pce,
		**link_time_t, **link_var_t, link_time, link_var; // link flow rate, link travel time, link travel time variance by time of day
	link_flow_t = new double *[num_tod];
	pre_link_flow_t = new double *[num_tod];
	left_link_flow_t = new double *[num_tod];
	right_link_flow_t = new double *[num_tod];
	mid_link_flow_t = new double *[num_tod];
	link_flow_t_pce = new double *[num_tod];
	link_time_t = new double *[num_tod];
	link_var_t = new double *[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		link_flow_t[t] = new double [num_link];
		pre_link_flow_t[t] = new double [num_link];
		left_link_flow_t[t] = new double [num_link];
		right_link_flow_t[t] = new double [num_link];
		mid_link_flow_t[t] = new double [num_link];
		link_flow_t_pce[t] = new double [num_link];
		link_time_t[t] = new double [num_link];
		link_var_t[t] = new double [num_link];
	}
	//
	double **link_length_t; // link length by time of day
	link_length_t = new double *[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		link_length_t[t] = new double [num_link];
	}
	//
	double ***link_toll_t_m; // link toll by time of day
	link_toll_t_m = new double **[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		link_toll_t_m[t] = new double *[num_mode];
		for (int m = 0; m < num_mode; m++)
		{
			link_toll_t_m[t][m] = new double [num_link];
		}
	}
	//
	double ****link_cost_t_c_m, ****link_flow_t_c_m, ****pre_link_flow_t_c_m, ****left_link_flow_t_c_m, ****right_link_flow_t_c_m, ****mid_link_flow_t_c_m;
	// link travel cost, link flow rate by time of day, user class, transportation mode
	link_cost_t_c_m = new double ***[num_tod];
	link_flow_t_c_m = new double ***[num_tod];
	left_link_flow_t_c_m = new double ***[num_tod];
	right_link_flow_t_c_m = new double ***[num_tod];
	mid_link_flow_t_c_m = new double ***[num_tod];
	pre_link_flow_t_c_m = new double ***[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		link_cost_t_c_m[t] = new double **[num_class];
		link_flow_t_c_m[t] = new double **[num_class];
		left_link_flow_t_c_m[t] = new double **[num_class];
		right_link_flow_t_c_m[t] = new double **[num_class];
		mid_link_flow_t_c_m[t] = new double **[num_class];
		pre_link_flow_t_c_m[t] = new double **[num_class];
		for (int c = 0; c < num_class; c++)
		{
			link_cost_t_c_m[t][c] = new double *[num_mode];
			link_flow_t_c_m[t][c] = new double *[num_mode];
			left_link_flow_t_c_m[t][c] = new double *[num_mode];
			right_link_flow_t_c_m[t][c] = new double *[num_mode];
			mid_link_flow_t_c_m[t][c] = new double *[num_mode];
			pre_link_flow_t_c_m[t][c] = new double *[num_mode];
			for (int m = 0; m < num_mode; m++)
			{
				link_cost_t_c_m[t][c][m] = new double [num_link];
				link_flow_t_c_m[t][c][m] = new double [num_link];
				left_link_flow_t_c_m[t][c][m] = new double [num_link];
				right_link_flow_t_c_m[t][c][m] = new double [num_link];
				mid_link_flow_t_c_m[t][c][m] = new double [num_link];
				pre_link_flow_t_c_m[t][c][m] = new double [num_link];
			}
		}
	}
	//
	long int ***od_predecessor_node_id; // predecessor node in constructing an O-D pair
	od_predecessor_node_id = new long int **[num_node];
	for (int n = 0; n < num_node; n++)
	{
		od_predecessor_node_id[n] = new long int *[num_node];
	}
	//
	int fr_node_index, to_node_index; // index of "from" node, index of "to" node
	int num_intervals = 20; // number of intervals used for line search of the traffic assignment procedure
	double derivative; // derivative used for line search of the traffic assignment procedure
	double alpha, beta, sigma, gamma, tau; // link parameters
	double accml_error, avg_error; // accumulated and average convergence errors
	double upsilon = 0.001, epsilon = 0.001; // travel demand procedure's and traffic assignment procedure's convergence criterion
	//double upsilon = 1.001, epsilon = 0.001; // travel demand procedure's and traffic assignment procedure's convergence criterion
	double outer_iteration, inner_iteration; // number of travel demand procedure's iteratrion and traffic assignment procedure's iteration
	double max_outer_iteration = 50, max_inner_iteration = 50; // maximum number of travel demand procedure's iteratrion and traffic assignment procedure's iteration
	double denominator;
	//
	double *****od_flow_t_c_m, *****od_cost_t_c_m, *****diff_od_cost_t_c_m; // O-D flow rate, O-D travel cost by time of day, user class, transportation mode in the alternative network
	od_flow_t_c_m = new double ****[num_tod];
	od_cost_t_c_m = new double ****[num_tod];
	diff_od_cost_t_c_m = new double ****[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		od_flow_t_c_m[t] = new double ***[num_class];
		od_cost_t_c_m[t] = new double ***[num_class];
		diff_od_cost_t_c_m[t] = new double ***[num_class];
		for (int c = 0; c < num_class; c++)
		{
			od_flow_t_c_m[t][c] = new double **[num_mode];
			od_cost_t_c_m[t][c] = new double **[num_mode];
			diff_od_cost_t_c_m[t][c] = new double **[num_mode];
			for (int m = 0; m < num_mode; m++)
			{
				od_flow_t_c_m[t][c][m] = new double *[num_node];
				od_cost_t_c_m[t][c][m] = new double *[num_node];
				diff_od_cost_t_c_m[t][c][m] = new double *[num_node];
				for (int n1 = 0; n1 < num_node; n1++)
				{
					od_flow_t_c_m[t][c][m][n1] = new double [num_node];
					od_cost_t_c_m[t][c][m][n1] = new double [num_node];
					diff_od_cost_t_c_m[t][c][m][n1] = new double [num_node];
				}
			}
		}
	}
	//
	double *****od_flow_t_c_m_old, *****od_cost_t_c_m_old; // O-D flow rate, O-D travel cost by time of day, user class, transportation mode in the original network
	od_flow_t_c_m_old = new double ****[num_tod];
	od_cost_t_c_m_old = new double ****[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		od_flow_t_c_m_old[t] = new double ***[num_class];
		od_cost_t_c_m_old[t] = new double ***[num_class];
		for (int c = 0; c < num_class; c++)
		{
			od_flow_t_c_m_old[t][c] = new double **[num_mode];
			od_cost_t_c_m_old[t][c] = new double **[num_mode];
			for (int m = 0; m < num_mode; m++)
			{
				od_flow_t_c_m_old[t][c][m] = new double *[num_node];
				od_cost_t_c_m_old[t][c][m] = new double *[num_node];
				for (int n1 = 0; n1 < num_node; n1++)
				{
					od_flow_t_c_m_old[t][c][m][n1] = new double [num_node];
					od_cost_t_c_m_old[t][c][m][n1] = new double [num_node];
				}
			}
		}
	}
	//
	double ***od_flow_t, ***od_flow_t_old; // O-D flow rate, O-D travel cost in the alternative network
	od_flow_t = new double **[num_tod];
	od_flow_t_old = new double **[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		od_flow_t[t] = new double *[num_node];
		od_flow_t_old[t] = new double *[num_node];
		for (int n1 = 0; n1 < num_node; n1++)
		{
			od_flow_t[t][n1] = new double [num_node];
			od_flow_t_old[t][n1] = new double [num_node];
		}
	}
	//
	double *total_od_flow_t, *total_od_flow_t_old;  // sum of O-D flow rates, sum of O-D travel costs in the original network
	total_od_flow_t = new double [num_tod];
	total_od_flow_t_old = new double [num_tod];
	//
	double ****base_od_flow_t_c, ****base_od_cost_t_c; // O-D flow rate, O-D travel cost by time of day, user class in the alternative network
	base_od_flow_t_c = new double ***[num_tod];
	base_od_cost_t_c = new double ***[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		base_od_flow_t_c[t] = new double **[num_class];
		base_od_cost_t_c[t] = new double **[num_class];
		for (int c = 0; c < num_class; c++)
		{
			base_od_flow_t_c[t][c] = new double *[num_node];
			base_od_cost_t_c[t][c] = new double *[num_node];
			for (int n1 = 0; n1 < num_node; n1++)
			{
				base_od_flow_t_c[t][c][n1] = new double [num_node];
				base_od_cost_t_c[t][c][n1] = new double [num_node];
			}
		}
	}
	//
	double ****base_od_flow_c_m, ****base_od_cost_c_m; // O-D flow rate, O-D travel cost by user class, transportation mode in the alternative network
	base_od_flow_c_m = new double ***[num_class];
	base_od_cost_c_m = new double ***[num_class];
	for (int c = 0; c < num_class; c++)
	{
		base_od_flow_c_m[c] = new double **[num_mode];
		base_od_cost_c_m[c] = new double **[num_mode];
		for (int m = 0; m < num_mode; m++)
		{
			base_od_flow_c_m[c][m] = new double *[num_node];
			base_od_cost_c_m[c][m] = new double *[num_node];
			for (int n1 = 0; n1 < num_node; n1++)
			{
				base_od_flow_c_m[c][m][n1] = new double [num_node];
				base_od_cost_c_m[c][m][n1] = new double [num_node];
			}
		}
	}
	//
	double ****od_flow_t_c, ****pre_od_flow_t_c, ****od_cost_t_c, ****pre_od_cost_t_c; // O-D flow rate, O-D travel cost by time of day, user class in the alternative network
	od_flow_t_c = new double ***[num_tod];
	pre_od_flow_t_c = new double ***[num_tod];
	od_cost_t_c = new double ***[num_tod];
	pre_od_cost_t_c = new double ***[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		od_flow_t_c[t] = new double **[num_class];
		pre_od_flow_t_c[t] = new double **[num_class];
		od_cost_t_c[t] = new double **[num_class];
		pre_od_cost_t_c[t] = new double **[num_class];
		for (int c = 0; c < num_class; c++)
		{
			od_flow_t_c[t][c] = new double *[num_node];
			pre_od_flow_t_c[t][c] = new double *[num_node];
			od_cost_t_c[t][c] = new double *[num_node];
			pre_od_cost_t_c[t][c] = new double *[num_node];
			for (int n1 = 0; n1 < num_node; n1++)
			{
				od_flow_t_c[t][c][n1] = new double [num_node];
				pre_od_flow_t_c[t][c][n1] = new double [num_node];
				od_cost_t_c[t][c][n1] = new double [num_node];
				pre_od_cost_t_c[t][c][n1] = new double [num_node];
			}
		}
	}
	//
	double ***od_trip_c, ***od_ecu_c; // O-D trip rate, O-D equivalent car unit by user class in the alternative network
	od_trip_c = new double **[num_class];
	od_ecu_c = new double **[num_class];
	for (int c = 0; c < num_class; c++)
	{
		od_trip_c[c] = new double *[num_node];
		od_ecu_c[c] = new double *[num_node];
		for (int n1 = 0; n1 < num_node; n1++)
		{
			od_trip_c[c][n1] = new double [num_node];
			od_ecu_c[c][n1] = new double [num_node];
		}
	}
	double ****od_flow_c_m, ****diff_od_flow_c_m, ****od_cost_c_m, ****diff_od_cost_c_m; // O-D flow rates, O-D travel costs by user class, transportation mode in the alternative network
	od_flow_c_m = new double ***[num_class];
	diff_od_flow_c_m = new double ***[num_class];
	od_cost_c_m = new double ***[num_class];
	diff_od_cost_c_m = new double ***[num_class];
	for (int c = 0; c < num_class; c++)
	{
		od_flow_c_m[c] = new double **[num_mode];
		diff_od_flow_c_m[c] = new double **[num_mode];
		od_cost_c_m[c] = new double **[num_mode];
		diff_od_cost_c_m[c] = new double **[num_mode];
		for (int m = 0; m < num_mode; m++)
		{
			od_flow_c_m[c][m] = new double *[num_node];
			diff_od_flow_c_m[c][m] = new double *[num_node];
			od_cost_c_m[c][m] = new double *[num_node];
			diff_od_cost_c_m[c][m] = new double *[num_node];
			for (int n1 = 0; n1 < num_node; n1++)
			{
				od_flow_c_m[c][m][n1] = new double [num_node];
				diff_od_flow_c_m[c][m][n1] = new double [num_node];
				od_cost_c_m[c][m][n1] = new double [num_node];
				diff_od_cost_c_m[c][m][n1] = new double [num_node];
			}
		}
	}
	//
	double ***link_cost_t_c, ***link_flow_t_c; // link travel cost, link flow rate by time of day, user class in the alternative network
	link_cost_t_c = new double **[num_tod];
	link_flow_t_c = new double **[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		link_cost_t_c[t] = new double *[num_class];
		link_flow_t_c[t] = new double *[num_class];
		for (int c = 0; c < num_class; c++)
		{
			link_cost_t_c[t][c] = new double [num_link];
			link_flow_t_c[t][c] = new double [num_link];
		}
	}
	//
	double ***link_cost_c_m, ***link_flow_c_m; // link travel cost, link flow rate by user class, transportation mode in the alternative network
	link_cost_c_m = new double **[num_class];
	link_flow_c_m = new double **[num_class];
	for (int c = 0; c < num_class; c++)
	{
		link_cost_c_m[c] = new double *[num_mode];
		link_flow_c_m[c] = new double *[num_mode];
		for (int m = 0; m < num_mode; m++)
		{
			link_cost_c_m[c][m] = new double [num_link];
			link_flow_c_m[c][m] = new double [num_link];
		}
	}
	//
	double ****mode_prob; // transportation mode choice probability
	mode_prob = new double ***[num_class];
	for (int c = 0; c < num_class; c++)
	{
		mode_prob[c] = new double **[num_mode];
		for (int m = 0; m < num_mode; m++)
		{
			mode_prob[c][m] = new double *[num_node];
			for (int n1 = 0; n1 < num_node; n1++)
			{
				mode_prob[c][m][n1] = new double [num_node];
			}
		}
	}
	//
	double *****tod_prob; // time-of-day choice probability
	tod_prob = new double ****[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		tod_prob[t] = new double ***[num_class];
		for (int c = 0; c < num_class; c++)
		{
			tod_prob[t][c] = new double **[num_mode];
			for (int m = 0; m < num_mode; m++)
			{
				tod_prob[t][c][m] = new double *[num_node];
				for (int n1 = 0; n1 < num_node; n1++)
				{
					tod_prob[t][c][m][n1] = new double [num_node];
				}
			}
		}
	}
	//
	double **last_link_flow_t; // link flow rate in the last iteration
	last_link_flow_t = new double *[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		last_link_flow_t[t] = new double [num_link];
	}
	double ***last_link_flow_t_c;
	last_link_flow_t_c = new double **[num_tod];
	for (int t = 0; t < num_tod; t++)
	{
		last_link_flow_t_c[t] = new double *[num_class];
		for (int c = 0; c < num_class; c++)
		{
			last_link_flow_t_c[t][c] = new double [num_link];
		}
	}

	for (int i = 0; i < 2; i++) // the i values of 0 and 1 indicate the original and alternative networks
	{

//	fprintf(myLogFile, "***************Starting %s network\n", (i == 0 ? "Base" : "Alternative"));
//	fprintf(myDebugFile, "***************Starting %s network\n", (i == 0 ? "Base" : "Alternative"));

	cout << "\n";
	if (i == 0)
	{
		cout << "   Estimate the travel demand pattern in the original network\n";
	}
	else
	{
		cout << "   Estimate the travel demand pattern in the alternative network\n";
	}

	// step -1 >> compute the initial generalized O-D cost
//	fprintf(myLogFile, "step -1\n"); fflush(myLogFile);
if (i == 0) {
	cout << "   ====================================================\n";
	cout << "   Iteration (outer loop): " << "0\n";
	cout << "   Estimate base generalized O-D costs ...\n";
	for (int t = 0; t < num_tod; t++)
	{
		ecu_t[t] = 0;
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				ecu_t[t] += proportion_c[c] * base_mode_prob[c][m] * ecu_m[m]; // calculate the aggregate ECU value
			}
		}
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						base_od_flow_t_c_m[t][c][m][n1][n2] = (base_od_flow_t[t][n1][n2] / ecu_t[t]) * proportion_c[c] * base_mode_prob[c][m] * ecu_m[m];
					}
				}
			}
		}
	}
	// calculate the base time-of-day choice probability
	for (int n1 = 0; n1 < num_node; n1++)
	{
		for (int n2 = 0; n2 < num_node; n2++)
		{
			base_od_flow = 0;
			for (int t = 0; t < num_tod; t++)
			{
				base_od_flow += base_od_flow_t[t][n1][n2];
			}
			for (int t = 0; t < num_tod; t++)
			{
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						if (base_od_flow * proportion_c[c] * base_mode_prob[c][m] > 0)
						{
							base_tod_prob[t][c][m][n1][n2] = base_od_flow_t_c_m[t][c][m][n1][n2] / (base_od_flow * proportion_c[c] * base_mode_prob[c][m]);
						}
						else
						{
							base_tod_prob[t][c][m][n1][n2] = 0;
						}
					}
				}
			}
		}
	}
}
	//
	if (i == 0)
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_flow_t[t][l] = 0;
				link_flow_t_pce[t][l] = 0;
				link_time_t[t][l] = link_attrbt_old_t[t][l][2];
				link_var_t[t][l] = link_attrbt_old_t[t][l][12];
			}
		}
		for (int p = 0; p < max_num_node; p++) {
			for (int od = 0; od <= 1; od++) {
				for (int t = 0; t < num_tod; t++) {
					for (int m = 0; m < num_mode; m++) {
						fixed_cost[p][t][m][od] = fixed_cost_base[p][t][m][od];
//						fprintf(myDebugFile, "Fixed cost for %d %d %d %d is %f\n", p, t, m, od, fixed_cost[p][t][m][od]);
					}
				}
			}
		}
	}
	else
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_flow_t[t][l] = 0;
				link_flow_t_pce[t][l] = 0;
				link_time_t[t][l] = link_attrbt_t[t][l][2];
				link_var_t[t][l] = link_attrbt_t[t][l][12];
			}
		}
		for (int p = 0; p < max_num_node; p++) {
			for (int od = 0; od <= 1; od++) {
				for (int t = 0; t < num_tod; t++) {
					for (int m = 0; m < num_mode; m++) {
						fixed_cost[p][t][m][od] = fixed_cost_alt[p][t][m][od];
//						fprintf(myDebugFile, "Fixed cost for %d %d %d %d is %f\n", p, t, m, od, fixed_cost[p][t][m][od]);
					}
				}
			}
		}

	}
	//
	if (i == 0)
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_length_t[t][l] = link_attrbt_old_t[t][l][1];
			}
		}
	}
	else
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_length_t[t][l] = link_attrbt_t[t][l][1];
			}
		}
	}
	//
	if (i == 0)
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int l = 0; l < num_link; l++)
				{
					link_toll_t_m[t][m][l] = link_attrbt_old_t[t][l][m + 6];
				}
			}
		}
	}
	else
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int l = 0; l < num_link; l++)
				{
					link_toll_t_m[t][m][l] = link_attrbt_t[t][l][m + 6];
				}
			}
		}
	}
	//
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int l = 0; l < num_link; l++)
				{
					link_cost_t_c_m[t][c][m][l] = vot_c[c] * link_time_t[t][l] + oc_m[m] * link_length_t[t][l] + link_toll_t_m[t][m][l];
					link_flow_t_c_m[t][c][m][l] = 0;
				}
			}
		}
	}
	// step (-1).0 >> initialization
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						if (n1 == n2)
						{
							od_cost[n1][n2] = 0;
							od_predecessor_node_id[n1][n2] = &node_id[n1];
						}
						else
						{
							od_cost[n1][n2] = positive_infinity * num_link;
						}
					}
				}
				for (int l = 0; l < num_link; l++)
				{
					fr_node_index = link_address[l][0] - node_id;
					to_node_index = link_address[l][1] - node_id;
					od_cost[fr_node_index][to_node_index] = link_cost_t_c_m[t][c][m][l];
					od_predecessor_node_id[fr_node_index][to_node_index] = link_address[l][0];
				}
				floyd_warshall_shortest_path(num_node, node_id, od_predecessor_node_id, od_cost);
				for (int n1 = 0; n1 < num_node; n1++) {
					for (int n2 = 0; n2 < num_node; n2++) {
					   //fprintf(myLogFile, "OD cost %d -> %d BEFORE fixed cost is %f", n1, n2, od_cost[n1][n2]);
						od_cost[n1][n2] += fixed_cost[n1][t][m][0] + fixed_cost[n2][t][m][1];
						//fprintf(myLogFile, " and after (+%f+%f) is %f\n", fixed_cost[n1][t][m][0], fixed_cost[n2][t][m][1], od_cost[n1][n2]); fflush(myLogFile);
					}
				}
				for (int l = 0; l < num_link; l++)
				{
					link_flow_t_c_m[t][c][m][l] = 0;
				}
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						if (base_od_flow_t_c_m[t][c][m][n1][n2] > 0)
						{
							to_node_index = n2;
							while (to_node_index != n1)
							{
								fr_node_index = od_predecessor_node_id[n1][to_node_index] - node_id;
								for (int l = 0; l < num_link; l++)
								{
									if (node_id[fr_node_index] == link_id[l][0] && node_id[to_node_index] == link_id[l][1])
									{
										link_flow_t_c_m[t][c][m][l] += base_od_flow_t_c_m[t][c][m][n1][n2]; // initial O-D flow rates
									}
								}
								to_node_index = fr_node_index;
							}
						}
					}
				}
			}
		}
		for (int l = 0; l < num_link; l++)
		{
			link_flow_t[t][l] = 0;
			link_flow_t_pce[t][l] = 0;
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					link_flow_t[t][l] += link_flow_t_c_m[t][c][m][l];
					link_flow_t_pce[t][l] += link_flow_t_c_m[t][c][m][l] * pce_m[m];
				}
			}
		}
	}
	// steps (-1).1-.5 >> cost update, direction finding, line search, solution update and convergence test
	for (int t = 0; t < num_tod; t++)
	{
		cout << "   ";
		inner_iteration = 0;
		do
		{
			inner_iteration++;
			// step (-1).1 >> cost update
			for (int l = 0; l < num_link; l++)
			{
				pre_link_flow_t[t][l] = link_flow_t[t][l];
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						pre_link_flow_t_c_m[t][c][m][l] = link_flow_t_c_m[t][c][m][l];
					}
				}
			}
			for (int l = 0; l < num_link; l++)
			{
				if (i == 0)
				{
					if (link_flow_t_pce[t][l] < MIN_FLOW_THRESHOLD) {
						link_time_t[t][l] = link_attrbt_old_t[t][l][2];
						link_var_t[t][l] = link_attrbt_old_t[t][l][12];
					} else {
						alpha = link_attrbt_old_t[t][l][3];
						beta = link_attrbt_old_t[t][l][4];
						link_time_t[t][l] = link_attrbt_old_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], beta));
						sigma = link_attrbt_old_t[t][l][13];
						gamma = link_attrbt_old_t[t][l][14];
						tau = link_attrbt_old_t[t][l][15];
						link_var_t[t][l] = link_attrbt_old_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], tau));
					}
				}
				else
				{
					if (link_flow_t_pce[t][l] < MIN_FLOW_THRESHOLD) {
						link_time_t[t][l] = link_attrbt_t[t][l][2];
						link_var_t[t][l] = link_attrbt_t[t][l][12];
					} else {
						alpha = link_attrbt_t[t][l][3];
						beta = link_attrbt_t[t][l][4];
						link_time_t[t][l] = link_attrbt_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], beta));
						sigma = link_attrbt_t[t][l][13];
						gamma = link_attrbt_t[t][l][14];
						tau = link_attrbt_t[t][l][15];
						link_var_t[t][l] = link_attrbt_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], tau));
					}
				}
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						link_cost_t_c_m[t][c][m][l] = vot_c[c] * link_time_t[t][l] + oc_m[m] * link_length_t[t][l] + link_toll_t_m[t][m][l];
					}
				}
			}
			// step (-1).2 >> direction finding
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							if (n1 == n2)
							{
								od_cost[n1][n2] = 0;
								od_predecessor_node_id[n1][n2] = &node_id[n1];
							}
							else
							{
								od_cost[n1][n2] = positive_infinity * num_link;
							}
						}
					}
					for (int l = 0; l < num_link; l++)
					{
						fr_node_index = link_address[l][0] - node_id;
						to_node_index = link_address[l][1] - node_id;
						od_cost[fr_node_index][to_node_index] = link_cost_t_c_m[t][c][m][l];
						od_predecessor_node_id[fr_node_index][to_node_index] = link_address[l][0];
					}
					floyd_warshall_shortest_path(num_node, node_id, od_predecessor_node_id, od_cost);
					for (int n1 = 0; n1 < num_node; n1++) {
						for (int n2 = 0; n2 < num_node; n2++) {
							od_cost[n1][n2] += fixed_cost[n1][t][m][0] + fixed_cost[n2][t][m][1];
						}
					}
					for (int l = 0; l < num_link; l++)
					{
						link_flow_t_c_m[t][c][m][l] = 0;
					}
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							if (base_od_flow_t_c_m[t][c][m][n1][n2] > 0)
							{
								to_node_index = n2;
								while (to_node_index != n1)
								{
									fr_node_index = od_predecessor_node_id[n1][to_node_index] - node_id;
									for (int l = 0; l < num_link; l++)
									{
										if (node_id[fr_node_index] == link_id[l][0] && node_id[to_node_index] == link_id[l][1])
										{
											link_flow_t_c_m[t][c][m][l] += base_od_flow_t_c_m[t][c][m][n1][n2]; // initial O-D flow rates
										}
									}
									to_node_index = fr_node_index;
								}
							}
						}
					}
				}
			}
			// step (-1).3 >> line search
			for (int l = 0; l < num_link; l++)
			{
				left_link_flow_t[t][l] = 0;
				right_link_flow_t[t][l] = 0;
			}
			for (int l = 0; l < num_link; l++)
			{
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						left_link_flow_t_c_m[t][c][m][l] = pre_link_flow_t_c_m[t][c][m][l];
						left_link_flow_t[t][l] += left_link_flow_t_c_m[t][c][m][l];
						right_link_flow_t_c_m[t][c][m][l] = link_flow_t_c_m[t][c][m][l];
						right_link_flow_t[t][l] += right_link_flow_t_c_m[t][c][m][l];
					}
				}
			}
			for (int n = 0; n < num_intervals; n++)
			{
				derivative = 0;
				for (int l = 0; l < num_link; l++)
				{
					if (i == 0)
					{
						alpha = link_attrbt_old_t[t][l][3];
						beta = link_attrbt_old_t[t][l][4];
						link_time = link_attrbt_old_t[t][l][2] * (1 +
							alpha * pow(0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_old_t[t][l][0], beta));
						sigma = link_attrbt_old_t[t][l][13];
						gamma = link_attrbt_old_t[t][l][14];
						tau = link_attrbt_old_t[t][l][15];
						link_var = link_attrbt_old_t[t][l][12] * (1 +
							sigma * pow(gamma + 0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_old_t[t][l][0], tau));
					}
					else
					{
						alpha = link_attrbt_t[t][l][3];
						beta = link_attrbt_t[t][l][4];
						link_time = link_attrbt_t[t][l][2] * (1 +
							alpha * pow(0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_t[t][l][0], beta));
						sigma = link_attrbt_t[t][l][13];
						gamma = link_attrbt_t[t][l][14];
						tau = link_attrbt_t[t][l][15];
						link_var = link_attrbt_t[t][l][12] * (1 +
							sigma * pow(gamma + 0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_t[t][l][0], tau));
					}
					derivative += link_time * (right_link_flow_t[t][l] - left_link_flow_t[t][l]);
					for (int c = 0; c < num_class; c++)
					{
						for (int m = 0; m < num_mode; m++)
						{
							derivative += (right_link_flow_t_c_m[t][c][m][l] - left_link_flow_t_c_m[t][c][m][l]) * oc_m[m] * link_length_t[t][l] / vot_c[c];
							derivative += (right_link_flow_t_c_m[t][c][m][l] - left_link_flow_t_c_m[t][c][m][l]) * link_toll_t_m[t][m][l] / vot_c[c];
							mid_link_flow_t_c_m[t][c][m][l] = 0.5 * (left_link_flow_t_c_m[t][c][m][l] + right_link_flow_t_c_m[t][c][m][l]);
						}
					}
				}
				if (derivative == 0 || abs(derivative) < 0.001)
				{
					break;
				}
				else if (derivative < 0)
				{
					for (int l = 0; l < num_link; l++)
					{
						left_link_flow_t[t][l] = 0;
						for (int c = 0; c < num_class; c++)
						{
							for (int m = 0; m < num_mode; m++)
							{
								left_link_flow_t_c_m[t][c][m][l] = mid_link_flow_t_c_m[t][c][m][l];
								left_link_flow_t[t][l] += left_link_flow_t_c_m[t][c][m][l];
							}
						}
					}
				}
				else if (derivative > 0)
				{
					for (int l = 0; l < num_link; l++)
					{
						right_link_flow_t[t][l] = 0;
						for (int c = 0; c < num_class; c++)
						{
							for (int m = 0; m < num_mode; m++)
							{
								right_link_flow_t_c_m[t][c][m][l] = mid_link_flow_t_c_m[t][c][m][l];
								right_link_flow_t[t][l] += right_link_flow_t_c_m[t][c][m][l];
							}
						}
					}
				}
			}
			// step (-1).4 >> solution update
			for (int l = 0; l < num_link; l++)
			{
				link_flow_t[t][l] = 0;
				link_flow_t_pce[t][l] = 0;
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						link_flow_t_c_m[t][c][m][l] = mid_link_flow_t_c_m[t][c][m][l];
						link_flow_t[t][l] += link_flow_t_c_m[t][c][m][l];
						link_flow_t_pce[t][l] += link_flow_t_c_m[t][c][m][l] * pce_m[m];
					}
				}
			}
			// step (-1).5 >> convergence test
			accml_error = 0;
			for (int l = 0; l < num_link; l++)
			{
				//fprintf(myLogFile, "%d %d %f %f\n", t, l, link_flow_t[t][l], pre_link_flow_t[t][l]);
				if (link_flow_t[t][l] > 0)
				{
					accml_error += abs(link_flow_t[t][l] - pre_link_flow_t[t][l]) / link_flow_t[t][l];
				}
			}
			avg_error = accml_error / (double) num_link;
			cout << ".";
			if (avg_error < epsilon || inner_iteration >= 30)
			{
				cout << "\n";
				break;
			}
		}
		while (1);
	}

	// step 0 >> initialization

	for (int t = 0; t < num_tod; t++)
	{
		for (int l = 0; l < num_link; l++)
		{
			if (i == 0)
			{
				alpha = link_attrbt_old_t[t][l][3];
				beta = link_attrbt_old_t[t][l][4];
				link_time_t[t][l] = link_attrbt_old_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], beta));
				sigma = link_attrbt_old_t[t][l][13];
				gamma = link_attrbt_old_t[t][l][14];
				tau = link_attrbt_old_t[t][l][15];
				link_var_t[t][l] = link_attrbt_old_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], tau));
			}
			else
			{
				alpha = link_attrbt_t[t][l][3];
				beta = link_attrbt_t[t][l][4];
				link_time_t[t][l] = link_attrbt_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], beta));
				sigma = link_attrbt_t[t][l][13];
				gamma = link_attrbt_t[t][l][14];
				tau = link_attrbt_t[t][l][15];
				link_var_t[t][l] = link_attrbt_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], tau));
			}
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					link_cost_t_c_m[t][c][m][l] = vot_c[c] * link_time_t[t][l] + oc_m[m] * link_length_t[t][l] + link_toll_t_m[t][m][l];
				}
			}
		}
	}
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						if (n1 == n2)
						{
							od_cost[n1][n2] = 0;
							od_predecessor_node_id[n1][n2] = &node_id[n1];
						}
						else
						{
							od_cost[n1][n2] = positive_infinity * num_link;
						}
					}
				}
				for (int l = 0; l < num_link; l++)
				{
					fr_node_index = link_address[l][0] - node_id;
					to_node_index = link_address[l][1] - node_id;
					od_cost[fr_node_index][to_node_index] = link_cost_t_c_m[t][c][m][l];
					od_predecessor_node_id[fr_node_index][to_node_index] = link_address[l][0];
				}
				floyd_warshall_shortest_path(num_node, node_id, od_predecessor_node_id, od_cost);
				for (int n1 = 0; n1 < num_node; n1++) {
					for (int n2 = 0; n2 < num_node; n2++) {
						od_cost[n1][n2] += fixed_cost[n1][t][m][0] + fixed_cost[n2][t][m][1];
					}
				}
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						base_od_cost_t_c_m[t][c][m][n1][n2] = od_cost[n1][n2];
					}
				}
			}
		}
	}
	//
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						od_flow_t_c_m[t][c][m][n1][n2] = base_od_flow_t_c_m[t][c][m][n1][n2];
					}
				}
			}
		}
	}
	//
if (i == 0) {
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					base_od_flow_t_c[t][c][n1][n2] = base_od_flow_t[t][n1][n2] * proportion_c[c];
				}
			}
		}
	}
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					if (base_od_flow_t_c[t][c][n1][n2] > 0)
					{
						base_od_cost_t_c[t][c][n1][n2] = 0;
						for (int m = 0; m < num_mode; m++)
						{
							base_od_cost_t_c[t][c][n1][n2] += base_od_flow_t_c_m[t][c][m][n1][n2] * base_od_cost_t_c_m[t][c][m][n1][n2];
						}
						base_od_cost_t_c[t][c][n1][n2] = base_od_cost_t_c[t][c][n1][n2] / base_od_flow_t_c[t][c][n1][n2];
					}
				}
			}
		}
	}
	//
	for (int c = 0; c < num_class; c++)
	{
		for (int m = 0; m < num_mode; m++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					base_od_flow_c_m[c][m][n1][n2] = 0;
					for (int t = 0; t < num_tod; t++)
					{
						base_od_flow_c_m[c][m][n1][n2] += base_od_flow_t_c_m[t][c][m][n1][n2];
					}
				}
			}
		}
	}
	for (int c = 0; c < num_class; c++)
	{
		for (int m = 0; m < num_mode; m++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					if (base_od_flow_c_m[c][m][n1][n2] > 0)
					{
						base_od_cost_c_m[c][m][n1][n2] = 0;
						for (int t = 0; t < num_tod; t++)
						{
							base_od_cost_c_m[c][m][n1][n2] += base_od_flow_t_c_m[t][c][m][n1][n2] * base_od_cost_t_c_m[t][c][m][n1][n2];
						}
						base_od_cost_c_m[c][m][n1][n2] = base_od_cost_c_m[c][m][n1][n2] / base_od_flow_c_m[c][m][n1][n2];
					}
					else
					{
						base_od_cost_c_m[c][m][n1][n2] = positive_infinity * num_link;
					}
				}
			}
		}
	}
	//
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					pre_od_flow_t_c[t][c][n1][n2] = base_od_flow_t_c[t][c][n1][n2];
					pre_od_cost_t_c[t][c][n1][n2] = base_od_cost_t_c[t][c][n1][n2];
				}
			}
		}
	}
}
	//
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_cost_t_c[t][c][l] = 0;
				link_flow_t_c[t][c][l] = 0;
				for (int m = 0; m < num_mode; m++)
				{
					link_cost_t_c[t][c][l] += link_cost_t_c_m[t][c][m][l] * link_flow_t_c_m[t][c][m][l];
					link_flow_t_c[t][c][l] += link_flow_t_c_m[t][c][m][l];
				}
				if (link_flow_t_c[t][c][l] > 0)
				{
					link_cost_t_c[t][c][l] = link_cost_t_c[t][c][l] / link_flow_t_c[t][c][l];
				}
				else
				{
					link_cost_t_c[t][c][l] = positive_infinity; //*****
				}
			}
		}
	}
	//
	for (int c = 0; c < num_class; c++)
	{
		for (int m = 0; m < num_mode; m++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_cost_c_m[c][m][l] = 0;
				link_flow_c_m[c][m][l] = 0;
				for (int t = 0; t < num_tod; t++)
				{
					link_cost_c_m[c][m][l] += link_cost_t_c_m[t][c][m][l] * link_flow_t_c_m[t][c][m][l];
					link_flow_c_m[c][m][l] += link_flow_t_c_m[t][c][m][l];
				}
				if (link_flow_c_m[c][m][l] > 0)
				{
					link_cost_c_m[c][m][l] = link_cost_c_m[c][m][l] / link_flow_c_m[c][m][l];
				}
				else
				{
					link_cost_c_m[c][m][l] = positive_infinity; //*****
				}
			}
		}
	}

	// steps 1-5 >> O-D matrix estimation, mode split, time-of-day split and traffic assignment

	outer_iteration = 0;
	do
	{
		outer_iteration++;
//		fprintf(myLogFile, "++++Beginning outer iteration %f\n", outer_iteration);
//		fprintf(myDebugFile, "++++Beginning outer iteration %f\n", outer_iteration);
		cout << "   ====================================================\n";
		cout << "   Iteration (outer loop): " << outer_iteration << "\n";
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				last_link_flow_t[t][l] = link_flow_t[t][l];
			}
		}
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				for (int l = 0; l < num_link; l++)
				{
					last_link_flow_t_c[t][c][l] = link_flow_t_c[t][c][l];
				}
			}
		}
		// step 1 >> estimate the elastic O-D matrix (trip generation and distribution)
		cout << "   Estimate the elastic origin-destination matrix ..." << "\n";
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				if (i == 0)
				{
					alpha = link_attrbt_old_t[t][l][3];
					beta = link_attrbt_old_t[t][l][4];
					link_time_t[t][l] = link_attrbt_old_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], beta));
					sigma = link_attrbt_old_t[t][l][13];
					gamma = link_attrbt_old_t[t][l][14];
					tau = link_attrbt_old_t[t][l][15];
					link_var_t[t][l] = link_attrbt_old_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], tau));
				}
				else
				{
					alpha = link_attrbt_t[t][l][3];
					beta = link_attrbt_t[t][l][4];
					link_time_t[t][l] = link_attrbt_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], beta));
					sigma = link_attrbt_t[t][l][13];
					gamma = link_attrbt_t[t][l][14];
					tau = link_attrbt_t[t][l][15];
					link_var_t[t][l] = link_attrbt_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], tau));
					if (isnan(link_time_t[t][l])) printf("ISNAN with flow %f and cap %f and FFT %f\n", link_flow_t_pce[t][l], link_attrbt_t[t][l][0], link_attrbt_t[t][l][2]);
				}
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						link_cost_t_c_m[t][c][m][l] = vot_c[c] * link_time_t[t][l] + oc_m[m] * link_length_t[t][l] + link_toll_t_m[t][m][l];
					}
				}
			}
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							if (n1 == n2)
							{
								od_cost[n1][n2] = 0;
								od_predecessor_node_id[n1][n2] = &node_id[n1];
							}
							else
							{
								od_cost[n1][n2] = positive_infinity;
							}
						}
					}
					for (int l = 0; l < num_link; l++)
					{
						fr_node_index = link_address[l][0] - node_id;
						to_node_index = link_address[l][1] - node_id;
						od_cost[fr_node_index][to_node_index] = link_cost_t_c_m[t][c][m][l];
						od_predecessor_node_id[fr_node_index][to_node_index] = link_address[l][0];
						if (isnan(od_cost[fr_node_index][to_node_index])) printf("NaN OD cost.\n");
					}
					if (m == 3 && totalroutes > 0)
						if (i == 0)
							transitZoneZoneCosts(num_node, node_id, od_cost, baseRouteMatrix_t[t], baseHeadway_t[t], numBaseRoutes_t[t], link_id, link_cost_t_c_m[t][c][0]);
						else
							transitZoneZoneCosts(num_node, node_id, od_cost, alternateRouteMatrix_t[t], alternateHeadway_t[t], numAlternateRoutes_t[t], link_id, link_cost_t_c_m[t][c][0]);
					else {
						floyd_warshall_shortest_path(num_node, node_id, od_predecessor_node_id, od_cost);
					}
					for (int n1 = 0; n1 < num_node; n1++) {
						for (int n2 = 0; n2 < num_node; n2++) {
							od_cost[n1][n2] += fixed_cost[n1][t][m][0] + fixed_cost[n2][t][m][1];
						}
					}
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							od_cost_t_c_m[t][c][m][n1][n2] = od_cost[n1][n2];
						}
					}
				}
			}
		}
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						od_cost_c_m[c][m][n1][n2] = delta;
						od_flow_c_m[c][m][n1][n2] = delta;
						if (base_od_flow_c_m[c][m][n1][n2] > 0)
						{
							for (int t = 0; t < num_tod; t++)
							{
								od_cost_c_m[c][m][n1][n2] += od_flow_t_c_m[t][c][m][n1][n2] * od_cost_t_c_m[t][c][m][n1][n2];
								if (isnan(od_cost_c_m[c][m][n1][n2])) printf("NaN after adding %f * %f\n", od_flow_t_c_m[t][c][m][n1][n2], od_cost_t_c_m[t][c][m][n1][n2]);
								od_flow_c_m[c][m][n1][n2] += od_flow_t_c_m[t][c][m][n1][n2];
							}
							if (od_flow_c_m[c][m][n1][n2] > 0)
							{
								od_cost_c_m[c][m][n1][n2] = od_cost_c_m[c][m][n1][n2] / od_flow_c_m[c][m][n1][n2];
							}
							else
							{
								od_cost_c_m[c][m][n1][n2] = base_od_cost_c_m[c][m][n1][n2];
							}
						}
						else
						{
							od_cost_c_m[c][m][n1][n2] = positive_infinity * num_link;
							if (isnan(od_cost_c_m[c][m][n1][n2])) { printf("It's NAN through infinity."); getchar(); }
						}
						if (isnan(od_cost_c_m[c][m][n1][n2])) { printf("Unrecoverable error."); exit(1); }
					}
				}
			}
		}
		// DEMAND FUNCTION
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						if (base_od_flow_c_m[c][m][n1][n2] > 0)
						{
							if (isnan(base_od_flow_c_m[c][m][n1][n2])) { printf("Unrecoverable error in base flow between %d and %d, class %d, mode %d", n1, n2, c, m); exit(1); }
							od_flow_c_m[c][m][n1][n2] = (1 - 1.0 / outer_iteration) * od_flow_c_m[c][m][n1][n2] +
							(1.0 / outer_iteration) * base_od_flow_c_m[c][m][n1][n2] *
								exp(elasticity * log(od_cost_c_m[c][m][n1][n2] / base_od_cost_c_m[c][m][n1][n2]));
						}
						else
						{
							od_flow_c_m[c][m][n1][n2] = 0;
						}
						if (isnan(od_flow_c_m[c][m][n1][n2])) { printf("Unrecoverable error in OD flow between %d and %d, class %d, mode %d.\n OD cost and base cost %f %f", n1, n2, c, m, od_cost_c_m[c][m][n1][n2], base_od_cost_c_m[c][m][n1][n2]); exit(1); }
					}
				}
//				fprintf(myDebugFile, "Flow and cost between 20 and 30 class %d mode %d is %f and %f with fixed costs %f %f\n", c, m, od_flow_c_m[c][m][20][30], od_cost_c_m[c][m][20][30], fixed_cost[20][0][m][0], fixed_cost[30][0][m][1]);
			}
		}


		// step 2 >> split traffic by mode (mode split)
		cout << "   Split travel demand by transportation mode ..." << "\n";
		// calculate the O-D trip rate for each class
		for (int c = 0; c < num_class; c++)
		{
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					od_trip_c[c][n1][n2] = 0;
					for (int m = 0; m < num_mode; m++)
					{
						od_trip_c[c][n1][n2] += od_flow_c_m[c][m][n1][n2] / ecu_m[m];
					}
				}
			}
		}
		// calculate the mode choice probability
		for (int n1 = 0; n1 < num_node; n1++)
		{
			for (int n2 = 0; n2 < num_node; n2++)
			{
				for (int c = 0; c < num_class; c++)
				{
					denominator = 0;
					double tempmodeprob;
					for (int m = 0; m < num_mode; m++)
					{
						diff_od_cost_c_m[c][m][n1][n2] = od_cost_c_m[c][m][n1][n2] - base_od_cost_c_m[c][m][n1][n2];
						denominator += base_mode_prob[c][m] * exp(-1 * mode_scale * diff_od_cost_c_m[c][m][n1][n2]);
						//if (c == 0 && n1 == 25 && n2 == 0) printf("Denominator is %f based on %f - %f\n", denominator, od_cost_c_m[c][m][n1][n2], base_od_cost_c_m[c][m][n1][n2]);
					}
					if (isinf(denominator)) {
						int mincost = positive_infinity;
						int minalt = -1;
						for (int m = 0; m < num_mode; m++) {
							if (od_cost_c_m[c][m][n1][n2] < mincost) {
								mincost = od_cost_c_m[c][m][n1][n2];
								minalt = m;
							}
						}
						for (int m = 0; m < num_mode; m++) {
							tempmodeprob = (m == minalt) ? 1 : 0;
							mode_prob[c][m][n1][n2] = tempmodeprob; // (1 - 1.0 / outer_iteration) * mode_prob[c][m][n1][n2] + (1.0 / outer_iteration) * tempmodeprob;
						}
					} else {
						if (denominator > 0)
						{
							for (int m = 0; m < num_mode; m++)
							{
								tempmodeprob = base_mode_prob[c][m] *
									exp(-1 * mode_scale * diff_od_cost_c_m[c][m][n1][n2]) / denominator;
								mode_prob[c][m][n1][n2] = tempmodeprob; // (1 - 1.0 / outer_iteration) * mode_prob[c][m][n1][n2] + (1.0 / outer_iteration) * tempmodeprob;
								if (isnan(mode_prob[c][m][n1][n2])) {  printf("Mode prob %d %d %d %d is NaN based on exp(-%f * %f)/%f\n", c, m, n1, n2, mode_scale, diff_od_cost_c_m[c][m][n1][n2], denominator); getchar(); }
							}
						}
						else
						{
							for (int m = 0; m < num_mode; m++)
							{
								tempmodeprob = base_mode_prob[c][m];
								mode_prob[c][m][n1][n2] = tempmodeprob; // (1 - 1.0 / outer_iteration) * mode_prob[c][m][n1][n2] + (1.0 / outer_iteration) * tempmodeprob;
							}
						}
					}
				}
			}
		}
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						od_flow_c_m[c][m][n1][n2] = (1 - 1.0 / outer_iteration) * od_flow_c_m[c][m][n1][n2] + (1.0 / outer_iteration) * od_trip_c[c][n1][n2] * mode_prob[c][m][n1][n2] * ecu_m[m];
						if (isnan(od_flow_c_m[c][m][n1][n2])) { printf("Flow  %d %d %d %d is NaN based on %f * %f * %f\n", c, m, n1, n2, od_trip_c[c][n1][n2], mode_prob[c][m][n1][n2], ecu_m[m]);  }
					}
				}
//				fprintf(myDebugFile, "Flow and cost between 20 and 30 class %d mode %d is %f and %f [mode split %f]\n", c, m, od_flow_c_m[c][m][20][30], od_cost_c_m[c][m][20][30], mode_prob[c][m][20][30]);
			}
		}
		// step 3 >> split traffic by time of day (time-of-day split)
		cout << "   Split travel demand by time of day ..." << "\n";
		for (int n1 = 0; n1 < num_node; n1++)
		{
			for (int n2 = 0; n2 < num_node; n2++)
			{
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						denominator = 0;
						for (int t = 0; t < num_tod; t++)
						{
							diff_od_cost_t_c_m[t][c][m][n1][n2] = od_cost_t_c_m[t][c][m][n1][n2] - base_od_cost_t_c_m[t][c][m][n1][n2];
							denominator += base_tod_prob[t][c][m][n1][n2] * exp(-1 * tod_scale * diff_od_cost_t_c_m[t][c][m][n1][n2]);
						}
						if (denominator > 0)
						{
							for (int t = 0; t < num_tod; t++)
							{
								tod_prob[t][c][m][n1][n2] = base_tod_prob[t][c][m][n1][n2] *
									exp(-1 * tod_scale * diff_od_cost_t_c_m[t][c][m][n1][n2]) / denominator;
							}
						}
						else
						{
							for (int t = 0; t < num_tod; t++)
							{
								tod_prob[t][c][m][n1][n2] = base_tod_prob[t][c][m][n1][n2];
							}
						}
					}
				}
			}
		}
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int t = 0; t < num_tod; t++)
				{
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							od_flow_t_c_m[t][c][m][n1][n2] = (1 - 1.0 / outer_iteration) * od_flow_t_c_m[t][c][m][n1][n2] + (1.0 / outer_iteration) * od_flow_c_m[c][m][n1][n2] * tod_prob[t][c][m][n1][n2];
							if (isnan(od_flow_t_c_m[t][c][m][n1][n2])) { printf("Flow %d %d %d is NaN based on %f * %f\n", t, c, m, od_flow_c_m[c][m][n1][n2], tod_prob[t][c][m][n1][n2]); getchar(); }
						}
					}
//					fprintf(myDebugFile, "Flow and cost between 20 and 30 class %d mode %d time %d is %f [time split %f]\n", c, m, t, od_flow_t_c_m[t][c][m][20][30], tod_prob[t][c][m][20][30]);
				}
			}
		}
		// step 4 >> split traffic by route (traffic assignment)
		cout << "   Split travel demand by route under equilibrium ..." << "\n";
		//num_mode--; // Transit is treated separately, easiest way to implement.  Be sure to reset after step 4 is done
		// step 4.0 >> initialization
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							if (n1 == n2)
							{
								od_cost[n1][n2] = 0;
								od_predecessor_node_id[n1][n2] = &node_id[n1];
							}
							else
							{
								od_cost[n1][n2] = positive_infinity * num_link;
							}
						}
					}
					for (int l = 0; l < num_link; l++)
					{
						fr_node_index = link_address[l][0] - node_id;
						to_node_index = link_address[l][1] - node_id;
						od_cost[fr_node_index][to_node_index] = link_cost_t_c_m[t][c][m][l];
						od_predecessor_node_id[fr_node_index][to_node_index] = link_address[l][0];
					}
					floyd_warshall_shortest_path(num_node, node_id, od_predecessor_node_id, od_cost);
					for (int n1 = 0; n1 < num_node; n1++) {
						for (int n2 = 0; n2 < num_node; n2++) {
							od_cost[n1][n2] += fixed_cost[n1][t][m][0] + fixed_cost[n2][t][m][1];
						}
					}
					for (int l = 0; l < num_link; l++)
					{
						link_flow_t_c_m[t][c][m][l] = 0;
					}
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							if (base_od_flow_t_c_m[t][c][m][n1][n2] > 0)
							{
								to_node_index = n2;
								while (to_node_index != n1)
								{
									fr_node_index = od_predecessor_node_id[n1][to_node_index] - node_id;
									for (int l = 0; l < num_link; l++)
									{
										if (node_id[fr_node_index] == link_id[l][0] && node_id[to_node_index] == link_id[l][1])
										{
											link_flow_t_c_m[t][c][m][l] += od_flow_t_c_m[t][c][m][n1][n2]; // initial O-D flow rates
										}
									}
									to_node_index = fr_node_index;
								}
							}
						}
					}
				}
			}
			for (int l = 0; l < num_link; l++)
			{
				link_flow_t[t][l] = 0;
				link_flow_t_pce[t][l] = 0;
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						if (isnan(link_flow_t_pce[t][l])) printf("NaN link flow on %d at %d.\n", l, t);
						link_flow_t[t][l] += link_flow_t_c_m[t][c][m][l];
						link_flow_t_pce[t][l] += link_flow_t_c_m[t][c][m][l] * pce_m[m];
						if (isnan(link_flow_t_pce[t][l])) printf("NaN link flow on %d at %d by adding %f * %f.\n", l, t, link_flow_t_c_m[t][c][m][l], pce_m[m]);
					}
				}
			}
		}

		// steps 4.1-.5 >> cost update, direction finding, line search, solution update and convergence test
		for (int t = 0; t < num_tod; t++)
		{
			cout << "   ";
			inner_iteration = 0;
			do
			{
				inner_iteration++;
				// step 4.1 >> cost update

				for (int l = 0; l < num_link; l++)
				{
					pre_link_flow_t[t][l] = link_flow_t[t][l];
					for (int c = 0; c < num_class; c++)
					{
						for (int m = 0; m < num_mode; m++)
						{
							pre_link_flow_t_c_m[t][c][m][l] = link_flow_t_c_m[t][c][m][l];
						}
					}
				}
				for (int l = 0; l < num_link; l++)
				{
					if (i == 0)
					{
						if (link_flow_t_pce[t][l] < MIN_FLOW_THRESHOLD) {
							link_time_t[t][l] = link_attrbt_old_t[t][l][2];
							link_var_t[t][l] = link_attrbt_old_t[t][l][12];
						} else {
							alpha = link_attrbt_old_t[t][l][3];
							beta = link_attrbt_old_t[t][l][4];
							link_time_t[t][l] = link_attrbt_old_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], beta));
							sigma = link_attrbt_old_t[t][l][13];
							gamma = link_attrbt_old_t[t][l][14];
							tau = link_attrbt_old_t[t][l][15];
							link_var_t[t][l] = link_attrbt_old_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_old_t[t][l][0], tau));
						}
					}
					else
					{
						if (link_flow_t_pce[t][l] < MIN_FLOW_THRESHOLD) {
							link_time_t[t][l] = link_attrbt_t[t][l][2];
							link_var_t[t][l] = link_attrbt_t[t][l][12];
						} else {
							alpha = link_attrbt_t[t][l][3];
							beta = link_attrbt_t[t][l][4];
							link_time_t[t][l] = link_attrbt_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], beta));
							sigma = link_attrbt_t[t][l][13];
							gamma = link_attrbt_t[t][l][14];
							tau = link_attrbt_t[t][l][15];
							link_var_t[t][l] = link_attrbt_t[t][l][12] * (1 + sigma * pow(gamma + link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], tau));
						}
					}
					for (int c = 0; c < num_class; c++)
					{
						for (int m = 0; m < num_mode; m++)
						{
							link_cost_t_c_m[t][c][m][l] = vot_c[c] * link_time_t[t][l] + oc_m[m] * link_length_t[t][l] + link_toll_t_m[t][m][l];

						}
					}

				}
				// step 4.2 >> direction finding

				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						for (int n1 = 0; n1 < num_node; n1++)
						{
							for (int n2 = 0; n2 < num_node; n2++)
							{
								if (n1 == n2)
								{
									od_cost[n1][n2] = 0;
									od_predecessor_node_id[n1][n2] = &node_id[n1];
								}
								else
								{
									od_cost[n1][n2] = positive_infinity;
								}
							}
						}
						for (int l = 0; l < num_link; l++)
						{
							fr_node_index = link_address[l][0] - node_id;
							to_node_index = link_address[l][1] - node_id;
							od_cost[fr_node_index][to_node_index] = link_cost_t_c_m[t][c][m][l];
							od_predecessor_node_id[fr_node_index][to_node_index] = link_address[l][0];
						}
						floyd_warshall_shortest_path(num_node, node_id, od_predecessor_node_id, od_cost);
						for (int n1 = 0; n1 < num_node; n1++) {
							for (int n2 = 0; n2 < num_node; n2++) {
								od_cost[n1][n2] += fixed_cost[n1][t][m][0] + fixed_cost[n2][t][m][1];
							}
						}
						for (int l = 0; l < num_link; l++)
						{
							link_flow_t_c_m[t][c][m][l] = 0;
						}
						for (int n1 = 0; n1 < num_node; n1++)
							{
							for (int n2 = 0; n2 < num_node; n2++)
							{
								if (base_od_flow_t_c_m[t][c][m][n1][n2] > 0)
								{
									to_node_index = n2;
									while (to_node_index != n1)

									{
										fr_node_index = od_predecessor_node_id[n1][to_node_index] - node_id;
										for (int l = 0; l < num_link; l++)
										{
											if (node_id[fr_node_index] == link_id[l][0] && node_id[to_node_index] == link_id[l][1])
											{
												link_flow_t_c_m[t][c][m][l] += od_flow_t_c_m[t][c][m][n1][n2];
											}
										}
										to_node_index = fr_node_index;

									}

								}
							}
						}
					}
				}
				// step 4.3 >> line search

				for (int l = 0; l < num_link; l++)
				{
					left_link_flow_t[t][l] = 0;
					right_link_flow_t[t][l] = 0;
				}
				for (int l = 0; l < num_link; l++)
				{
					for (int c = 0; c < num_class; c++)
					{
						for (int m = 0; m < num_mode; m++)
						{
							left_link_flow_t_c_m[t][c][m][l] = pre_link_flow_t_c_m[t][c][m][l];
							left_link_flow_t[t][l] += left_link_flow_t_c_m[t][c][m][l];
							right_link_flow_t_c_m[t][c][m][l] = link_flow_t_c_m[t][c][m][l];
							right_link_flow_t[t][l] += right_link_flow_t_c_m[t][c][m][l];
						}
					}
				}
				for (int n = 0; n < num_intervals; n++)
				{
					derivative = 0;
					for (int l = 0; l < num_link; l++)
					{
						if (i == 0)
						{
							alpha = link_attrbt_old_t[t][l][3];
							beta = link_attrbt_old_t[t][l][4];
							link_time = link_attrbt_old_t[t][l][2] * (1 +
								alpha * pow(0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_old_t[t][l][0], beta));
							sigma = link_attrbt_old_t[t][l][13];
							gamma = link_attrbt_old_t[t][l][14];
							tau = link_attrbt_old_t[t][l][15];
							link_var = link_attrbt_old_t[t][l][12] * (1 +
								sigma * pow(gamma + 0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_old_t[t][l][0], tau));
						}
						else
						{
							alpha = link_attrbt_t[t][l][3];
							beta = link_attrbt_t[t][l][4];
							link_time = link_attrbt_t[t][l][2] * (1 +
								alpha * pow(0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_t[t][l][0], beta));
							sigma = link_attrbt_t[t][l][13];
							gamma = link_attrbt_t[t][l][14];
							tau = link_attrbt_t[t][l][15];
							link_var = link_attrbt_t[t][l][12] * (1 +
								sigma * pow(gamma + 0.5 * (right_link_flow_t[t][l] + left_link_flow_t[t][l]) / link_attrbt_t[t][l][0], tau));
						}
						derivative += link_time * (right_link_flow_t[t][l] - left_link_flow_t[t][l]);
						for (int c = 0; c < num_class; c++)
						{
							for (int m = 0; m < num_mode; m++)
							{
								derivative += (right_link_flow_t_c_m[t][c][m][l] - left_link_flow_t_c_m[t][c][m][l]) * oc_m[m] * link_length_t[t][l] / vot_c[c];
								derivative += (right_link_flow_t_c_m[t][c][m][l] - left_link_flow_t_c_m[t][c][m][l]) * link_toll_t_m[t][m][l] / vot_c[c];
								mid_link_flow_t_c_m[t][c][m][l] = 0.5 * (left_link_flow_t_c_m[t][c][m][l] + right_link_flow_t_c_m[t][c][m][l]);
							}
						}
					}
					if (derivative == 0 || abs(derivative) < 0.001)
					{
						break;
					}
					else if (derivative < 0)
					{
						for (int l = 0; l < num_link; l++)
						{
							left_link_flow_t[t][l] = 0;
							for (int c = 0; c < num_class; c++)
							{
								for (int m = 0; m < num_mode; m++)
								{
									left_link_flow_t_c_m[t][c][m][l] = mid_link_flow_t_c_m[t][c][m][l];
									left_link_flow_t[t][l] += left_link_flow_t_c_m[t][c][m][l];
								}
							}
						}
					}
					else if (derivative > 0)
					{
						for (int l = 0; l < num_link; l++)
						{
							right_link_flow_t[t][l] = 0;
							for (int c = 0; c < num_class; c++)
							{
								for (int m = 0; m < num_mode; m++)
								{
									right_link_flow_t_c_m[t][c][m][l] = mid_link_flow_t_c_m[t][c][m][l];
									right_link_flow_t[t][l] += right_link_flow_t_c_m[t][c][m][l];
								}
							}
						}
					}
				}
				// step 4.4 >> solution update
				for (int l = 0; l < num_link; l++)
				{
					link_flow_t[t][l] = 0;
					link_flow_t_pce[t][l] = 0;
					for (int c = 0; c < num_class; c++)
					{
						link_flow_t_c[t][c][l] = 0;
						for (int m = 0; m < num_mode; m++)
						{
							link_flow_t_c_m[t][c][m][l] = mid_link_flow_t_c_m[t][c][m][l];
							link_flow_t[t][l] += link_flow_t_c_m[t][c][m][l];
							link_flow_t_pce[t][l] += link_flow_t_c_m[t][c][m][l] * pce_m[m];
							link_flow_t_c[t][c][l] += link_flow_t_c_m[t][c][m][l];
						}
					}
				}
				// step 4.5 >> convergence test
				accml_error = 0;
				for (int l = 0; l < num_link; l++)
				{
					if (link_flow_t[t][l] > 0)
					{
						accml_error += abs(link_flow_t[t][l] - pre_link_flow_t[t][l]) / link_flow_t[t][l];
					}
				}
				avg_error = accml_error / (double) num_link;
				cout << ".";
				//cout << "    Average link flow difference (inner loop): " << avg_error << "\n";
				if (avg_error < epsilon || inner_iteration == max_inner_iteration)
				{
					cout << "\n";
					break;
				}
			}
			while (1);
		}

		// step 5 >> convergence test
		{
			double CVMT = 0;
			for (int t = 0; t < num_tod; t++)
			{
				double VMT = 0;
				for (int l = 0; l < num_link; l++)
				{
					VMT += 24 * duration_t[t] * link_attrbt_t[t][l][1] * link_flow_t[t][l];
				}
//				fprintf(myDebugFile, "VMT in TOD %d (%f) before averaging is %f\n", t, duration_t[t], VMT);
//				printf("VMT in TOD %d (%f) before averaging is %f\n", t, duration_t[t], VMT);
				CVMT += VMT;
			}
//			fprintf(myDebugFile, "TOTAL annual VMT (before averaging) is %f\n", 365 * CVMT / 1e6);
//			printf("TOTAL annual VMT (before averaging) is %f\n", 365 * CVMT / 1e6);
		}
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				link_flow_t[t][l] = last_link_flow_t[t][l] +
					(1 / (double) outer_iteration) * (link_flow_t[t][l] - last_link_flow_t[t][l]);
			}
		}
		{
			double CVMT = 0;
			for (int t = 0; t < num_tod; t++)
			{
				double VMT = 0;
				for (int l = 0; l < num_link; l++)
				{
					VMT += 24 * duration_t[t] * link_attrbt_t[t][l][1] * link_flow_t[t][l];
				}
//				fprintf(myDebugFile, "VMT in TOD %d (%f) after averaging is %f\n", t, duration_t[t], VMT);
//				printf("VMT in TOD %d (%f) after averaging is %f\n", t, duration_t[t], VMT);
				CVMT += VMT;
			}
//			fprintf(myDebugFile, "TOTAL annual VMT (after averaging) is %f\n", 365 * CVMT / 1e6);
//			printf("TOTAL annual VMT (after averaging) is %f\n", 365 * CVMT/ 1e6);
		}
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				for (int l = 0; l < num_link; l++)
				{
					link_flow_t_c[t][c][l] = last_link_flow_t_c[t][c][l] +
						(1 / (double) outer_iteration) * (link_flow_t_c[t][c][l] - last_link_flow_t_c[t][c][l]);
				}
			}
		}
		accml_error = 0;
		for (int t = 0; t < num_tod; t++)
		{
			for (int l = 0; l < num_link; l++)
			{
				if (link_flow_t[t][l] > 0)
				{
					accml_error += abs(link_flow_t[t][l] - last_link_flow_t[t][l]) / link_flow_t[t][l];
				}
				//fprintf(myLogFile, "%f %d %d %f %f\n", accml_error, t, l, link_flow_t[t][l], last_link_flow_t[t][l]);
			}
		}
		avg_error = accml_error / ((double) num_link * (double) num_tod);
		cout << "   Average link flow difference (outer loop): " << avg_error << "\n";
//		fprintf(myLogFile, "Average: %f\n", avg_error);
		if (avg_error < upsilon || outer_iteration == max_outer_iteration || outer_iteration == max_outer_iteration_user)
		{
			break;
		}
		//num_mode++; // Reset action at start of step 4
	}
	while (1);
	if (i == 0)
	{
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					for (int n1 = 0; n1 < num_node; n1++)
					{
						for (int n2 = 0; n2 < num_node; n2++)
						{
							od_flow_t_c_m_old[t][c][m][n1][n2] = od_flow_t_c_m[t][c][m][n1][n2];
							od_cost_t_c_m_old[t][c][m][n1][n2] = od_cost_t_c_m[t][c][m][n1][n2];

						}
					}
				}
			}
		}
	}
	if (i == 0)
	{
		for (int t = 0; t < num_tod; t++)
		{
			total_od_flow_t_old[t] = 0;
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					od_flow_t_old[t][n1][n2] = 0;
					for (int c = 0; c < num_class; c++)
					{
						for (int m = 0; m < num_mode; m++)
						{
							od_flow_t_old[t][n1][n2] += od_flow_t_c_m_old[t][c][m][n1][n2];
						}
					}
					total_od_flow_t_old[t] += od_flow_t_old[t][n1][n2];
				}
			}
		}
	}
	else
	{
		for (int t = 0; t < num_tod; t++)
		{
			total_od_flow_t[t] = 0;
			for (int n1 = 0; n1 < num_node; n1++)
			{
				for (int n2 = 0; n2 < num_node; n2++)
				{
					od_flow_t[t][n1][n2] = 0;
					for (int c = 0; c < num_class; c++)
					{
						for (int m = 0; m < num_mode; m++)
						{
							od_flow_t[t][n1][n2] += od_flow_t_c_m[t][c][m][n1][n2];
						}
					}
					total_od_flow_t[t] += od_flow_t[t][n1][n2];
				}
			}
		}
	}

	}

	// step 6 >> write the result into the output files

	// write the link flow rates of the alternative network
	if (!argv[5 + num_tod + num_tod + num_tod])
	{
		strcpy(argv[5 + num_tod + num_tod + num_tod], "default_link_flow.txt");
		endtime = clock() / CLOCKS_PER_SEC;
		cout << "   Clock" << "\t\t" << endtime << " : Argument(s) missing...Default output file name used" << "\n";
	}
	/*
	if (num_scenarios == 1)
	{
		ofstream linkflow_out(argv[3 + num_tod + num_tod + num_tod]);
		linkflow_out << "<NUMBER OF ZONES> " << num_node << "\n";
		linkflow_out << "<NUMBER OF LINKS> " << num_link << "\n";
		linkflow_out << "<END OF METADATA>" << "\n\n\n";
		linkflow_out.width(10); linkflow_out << left << "~";
		linkflow_out.width(10); linkflow_out << left << "Tail";
		linkflow_out.width(10); linkflow_out << left << "Head";
		linkflow_out.width(10); linkflow_out << left << ":";
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				linkflow_out.width(3); linkflow_out << left << "TOD";
				linkflow_out.width(1); linkflow_out << left << t + 1;
				linkflow_out.width(6); linkflow_out << left << "_Class";
				linkflow_out.width(4); linkflow_out << left << c + 1;
			}
		}
		linkflow_out << ";\n";
		for (int l = 0; l < num_link; l++)
		{
			linkflow_out.width(10); linkflow_out << left << "";
			linkflow_out.width(10); linkflow_out << left << link_id[l][0];
			linkflow_out.width(10); linkflow_out << left << link_id[l][1];
			linkflow_out.width(10); linkflow_out << left << ":";
			for (int t = 0; t < num_tod; t++)
			{
				for (int c = 0; c < num_class; c++)
				{
					linkflow_out.width(10); linkflow_out << left << setiosflags(ios::fixed) << setprecision(1) << link_flow_t_c[t][c][l];
				}
			}
			linkflow_out << left << ";\n";
		}
		linkflow_out.close();
	}
	else
	{
		string argv_series = argv[3 + num_tod + num_tod + num_tod];
		int q = argv_series.find(".", 0);
		argv_series.erase(q, 4);
		stringstream argv_stream;
		argv_stream << scenario + 1;
		argv_series += "_";
		argv_series += argv_stream.str();
		argv_series += ".txt";
		const char *argv_series_c = argv_series.c_str();
		ofstream linkflow_out(argv_series_c);
		linkflow_out << "<NUMBER OF ZONES> " << num_node << "\n";
		linkflow_out << "<NUMBER OF LINKS> " << num_link << "\n";
		linkflow_out << "<END OF METADATA>" << "\n\n\n";
		linkflow_out.width(10); linkflow_out << left << "~";
		linkflow_out.width(10); linkflow_out << left << "Tail";
		linkflow_out.width(10); linkflow_out << left << "Head";
		linkflow_out.width(10); linkflow_out << left << ":";
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				linkflow_out.width(3); linkflow_out << left << "TOD";
				linkflow_out.width(1); linkflow_out << left << t + 1;
				linkflow_out.width(6); linkflow_out << left << "_Class";
				linkflow_out.width(4); linkflow_out << left << c + 1;
			}
		}
		linkflow_out << ";\n";
		for (int l = 0; l < num_link; l++)
		{
			linkflow_out.width(10); linkflow_out << left << "";
			linkflow_out.width(10); linkflow_out << left << link_id[l][0];
			linkflow_out.width(10); linkflow_out << left << link_id[l][1];
			linkflow_out.width(10); linkflow_out << left << ":";
			for (int t = 0; t < num_tod; t++)
			{
				for (int c = 0; c < num_class; c++)
				{
					linkflow_out.width(10); linkflow_out << left << setiosflags(ios::fixed) << setprecision(1) << link_flow_t_c[t][c][l];
				}
			}
			linkflow_out << left << ";\n";
		}
		linkflow_out.close();
	}
	*/
	if (num_scenarios == 1)
	{
		ofstream linkflow_out(argv[5 + num_tod + num_tod + num_tod]);
		linkflow_out << "<NUMBER OF ZONES> " << num_node << "\n";
		linkflow_out << "<NUMBER OF LINKS> " << num_link << "\n";
		linkflow_out << "<END OF METADATA>" << "\n\n\n";
		linkflow_out.width(8); linkflow_out << left << "~";
		linkflow_out.width(8); linkflow_out << left << "Tail";
		linkflow_out.width(8); linkflow_out << left << "Head";
		linkflow_out.width(8); linkflow_out << left << ":";
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					linkflow_out.width(3); linkflow_out << left << "TOD";
					linkflow_out.width(1); linkflow_out << left << t + 1;
					linkflow_out.width(4); linkflow_out << left << "_CLS";
					linkflow_out.width(1); linkflow_out << left << c + 1;
					linkflow_out.width(4); linkflow_out << left << "_MOD";
					linkflow_out.width(4); linkflow_out << left << m + 1;
				}
			}
		}
		linkflow_out << ";\n";
		for (int l = 0; l < num_link; l++)
		{
			linkflow_out.width(8); linkflow_out << left << "";
			linkflow_out.width(8); linkflow_out << left << link_id[l][0];
			linkflow_out.width(8); linkflow_out << left << link_id[l][1];
			linkflow_out.width(8); linkflow_out << left << ":";
			for (int t = 0; t < num_tod; t++)
			{
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						linkflow_out.width(8); linkflow_out << left << setiosflags(ios::fixed) << setprecision(1) << link_flow_t_c_m[t][c][m][l];
					}
				}
			}
			linkflow_out << left << ";\n";
		}
		linkflow_out.close();
	}
	else
	{
		string argv_series = argv[5 + num_tod + num_tod + num_tod];
		int q = argv_series.find(".", 0);
		argv_series.erase(q, 4);
		stringstream argv_stream;
		argv_stream << scenario + 1;
		argv_series += "_";
		argv_series += argv_stream.str();
		argv_series += ".txt";
		const char *argv_series_c = argv_series.c_str();
		ofstream linkflow_out(argv_series_c);
		linkflow_out << "<NUMBER OF ZONES> " << num_node << "\n";
		linkflow_out << "<NUMBER OF LINKS> " << num_link << "\n";
		linkflow_out << "<END OF METADATA>" << "\n\n\n";
		linkflow_out.width(8); linkflow_out << left << "~";
		linkflow_out.width(8); linkflow_out << left << "Tail";
		linkflow_out.width(8); linkflow_out << left << "Head";
		linkflow_out.width(8); linkflow_out << left << ":";
		for (int t = 0; t < num_tod; t++)
		{
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_node; m++)
				{
					linkflow_out.width(3); linkflow_out << left << "TOD";
					linkflow_out.width(1); linkflow_out << left << t + 1;
					linkflow_out.width(4); linkflow_out << left << "_CLS";
					linkflow_out.width(1); linkflow_out << left << c + 1;
					linkflow_out.width(4); linkflow_out << left << "_MOD";
					linkflow_out.width(4); linkflow_out << left << m + 1;
				}
			}
		}
		linkflow_out << ";\n";
		for (int l = 0; l < num_link; l++)
		{
			linkflow_out.width(8); linkflow_out << left << "";
			linkflow_out.width(8); linkflow_out << left << link_id[l][0];
			linkflow_out.width(8); linkflow_out << left << link_id[l][1];
			linkflow_out.width(8); linkflow_out << left << ":";
			for (int t = 0; t < num_tod; t++)
			{
				for (int c = 0; c < num_class; c++)
				{
					for (int m = 0; m < num_mode; m++)
					{
						linkflow_out.width(8); linkflow_out << left << setiosflags(ios::fixed) << setprecision(1) << link_flow_t_c_m[t][c][m][l];
					}
				}
			}
			linkflow_out << left << ";\n";
		}
		linkflow_out.close();
	}

	// calculate the summary performance measures for the alternative network
	for (int t = 0; t < num_tod; t++)
	{
		for (int l = 0; l < num_link; l++)
		{
			alpha = link_attrbt_t[t][l][3];
			beta = link_attrbt_t[t][l][4];
			link_flow_t[t][l] = 0;
			link_flow_t_pce[t][l] = 0;
			for (int c = 0; c < num_class; c++)
			{
				link_flow_t_c[t][c][l] = 0;
				for (int m = 0; m < num_mode; m++)
				{
					link_flow_t[t][l] += link_flow_t_c_m[t][c][m][l];
					link_flow_t_pce[t][l] += link_flow_t_c_m[t][c][m][l] * pce_m[m];
					link_flow_t_c[t][c][l] += link_flow_t_c_m[t][c][m][l];
				}
			}
			link_time_t[t][l] = link_attrbt_t[t][l][2] * (1 + alpha * pow(link_flow_t_pce[t][l] / link_attrbt_t[t][l][0], beta));
		}
	}
	//double user_surplus_old = 0;
	double user_surplus = 0;
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						if (base_od_flow_t_c_m[t][c][m][n1][n2] > 0)
						{
							//user_surplus_old += (base_od_flow_t_c_old[t][c][n1][n2] / od_ecu_c[c][n1][n2]) * (base_od_cost_t_c_old[t][c][n1][n2] - od_cost_t_c[t][c][n1][n2])
							//	+ elasticity * (2 / od_ecu_c[c][n1][n2] * base_od_cost_t_c_old[t][c][n1][n2]) * (pow(od_flow_t_c[t][c][n1][n2], 1 / elasticity + 1) - pow(base_od_flow_t_c[t][c][n1][n2], 1 / elasticity + 1))
							//	/ (2 * (1 + elasticity)) / pow(od_flow_t_c[t][c][n1][n2], 1 / elasticity)
							//	+ (base_od_flow_t_c_old[t][c][n1][n2] - od_flow_t_c[t][c][n1][n2]) / od_ecu_c[c][n1][n2] * od_cost_t_c[t][c][n1][n2];
							user_surplus += 24 * duration_t[t] * 0.5 * (od_flow_t_c_m_old[t][c][m][n1][n2] + od_flow_t_c_m[t][c][m][n1][n2])
								* (od_cost_t_c_m_old[t][c][m][n1][n2] - od_cost_t_c_m[t][c][m][n1][n2]) / ecu_m[m];
//fprintf(myLogFile, "Tod %d class %d mode %d OD %d -> %d: %f from %f from %f %f %f %f\n", t, c, m, n1, n2, user_surplus, 24 * duration_t[t] * 0.5 * (od_flow_t_c_m_old[t][c][m][n1][n2] + od_flow_t_c_m[t][c][m][n1][n2])
//		* (od_cost_t_c_m_old[t][c][m][n1][n2] - od_cost_t_c_m[t][c][m][n1][n2]) / ecu_m[m], od_flow_t_c_m_old[t][c][m][n1][n2], od_flow_t_c_m[t][c][m][n1][n2], od_cost_t_c_m_old[t][c][m][n1][n2], od_cost_t_c_m[t][c][m][n1][n2]);

						}
					}
				}
			}
		}
	}
	double travel_cost = 0, travel_time = 0, vht = 0, vmt = 0;
	for (int t = 0; t < num_tod; t++)
	{
		for (int l = 0; l < num_link; l++)
		{
			vht += 24 * duration_t[t] * link_time_t[t][l] * link_flow_t[t][l];
			vmt += 24 * duration_t[t] * link_attrbt_t[t][l][1] * link_flow_t[t][l];
			for (int c = 0; c < num_class; c++)
			{
				for (int m = 0; m < num_mode; m++)
				{
					travel_cost += 24 * duration_t[t] * link_cost_t_c_m[t][c][m][l] * link_flow_t_c_m[t][c][m][l] / ecu_m[m];
					travel_time += 24 * duration_t[t] * link_time_t[t][l] * link_flow_t_c_m[t][c][m][l] / ecu_m[m];
				}
			}
		}
	}
	user_surplus_series[scenario] = user_surplus;
	travel_cost_series[scenario] = travel_cost;
	travel_time_series[scenario] = travel_time;
	vht_series[scenario] = vht;
	vmt_series[scenario] = vmt;
	for (int c = 0; c < num_class; c++) {
		for (int t = 0; t < num_tod; t++) {
			total_trip_rate_series[c][t] = 0;
			for (int m = 0; m < num_mode; m++) {
				for (int n1 = 0; n1 < num_node; n1++) {
					for (int n2 = 0; n2 < num_node; n2++) {
						total_trip_rate_series[c][t] += od_flow_t_c_m[t][c][m][n1][n2];
					}
				}
			}
		}
	}

	// write the summary performance measures for the alternative network
	if (!argv[6 + num_tod + num_tod + num_tod])
	{
		strcpy(argv[6 + num_tod + num_tod + num_tod], "default_summary.txt");
		endtime = clock() / CLOCKS_PER_SEC;
		cout << "   Clock" << "\t\t" << endtime << " : Argument(s) missing...Default output file name used" << "\n";
	}
	if (num_scenarios == 1)
	{
		ofstream summary_out(argv[6 + num_tod + num_tod + num_tod]);
		summary_out << "<NUMBER OF ZONES> " << num_node << "\n";
		summary_out << "<NUMBER OF LINKS> " << num_link << "\n";
		summary_out << "<END OF METADATA>" << "\n\n\n";
		summary_out << "Traveler surplus change: " << user_surplus << " dollars" << "\n";
		summary_out << "Traveler travel cost: " << travel_cost << " dollars" << "\n";
		summary_out << "Traveler travel time: " << travel_time << " hours" << "\n";
		summary_out << "Vehicle hours traveled: " << vht << " hours" << "\n";
		summary_out << "Vehicle miles traveled: " << vmt << " miles" << "\n";
		summary_out << "\n";
		for (int c = 0; c < num_class; c++) {
			for (int t = 0; t < num_tod; t++) {
				summary_out << " Total trips for class " << c+1 << " during TOD " << t+1 << ": " << total_trip_rate_series[c][t] << "\n";
			}
			summary_out << "\n";
		}
		summary_out.close();
	}
	else
	{
		string argv_series = argv[6 + num_tod + num_tod + num_tod];
		int g = argv_series.find(".", 0);
		argv_series.erase(g, 4);
		stringstream argv_stream;
		argv_stream << scenario + 1;
		argv_series += "_";
		argv_series += argv_stream.str();
		argv_series += ".txt";
		const char *argv_series_c = argv_series.c_str();
		ofstream summary_out(argv_series_c);
		summary_out << "<NUMBER OF ZONES> " << num_node << "\n";
		summary_out << "<NUMBER OF LINKS> " << num_link << "\n";
		summary_out << "<END OF METADATA>" << "\n\n\n";
		summary_out << "Traveler surplus change: " << user_surplus << " dollars" << "\n";
		summary_out << "Traveler travel cost: " << travel_cost << " dollars" << "\n";
		summary_out << "Traveler travel time: " << travel_time << " hours" << "\n";
		summary_out << "Vehicle hours traveled: " << vht << " hours" << "\n";
		summary_out << "Vehicle miles traveled: " << vmt << " miles" << "\n";
		summary_out << "\n";
		for (int c = 0; c < num_class; c++) {
			for (int t = 0; t < num_tod; t++) {
				summary_out << " Total trips for class " << c+1 << " during TOD " << t+1 << ": " << total_trip_rate_series[c][t] << "\n";
			}
			summary_out << "\n";
		}
		summary_out.close();
	}

	// write the O-D trip rates of the alternative network
	if (argv[7 + num_tod + num_tod + num_tod])
	{
		if (num_scenarios == 1)
		{
			for (int t = 0; t < num_tod; t++)
			{
				ofstream odflow_out(argv[t + 7 + num_tod + num_tod + num_tod]);
				odflow_out << "<NUMBER OF ZONES> " << node_id[num_node - 1] << "\n";
				odflow_out << "<TOTAL OD FLOW> " << total_od_flow_t[t] << "\n";
				odflow_out << "<END OF METADATA>" << "\n\n\n";
				for (int fr = 0; fr < num_node; fr++)
				{
					odflow_out << "Origin" << " " << "\t" << node_id[fr] << " \n";
					int count = 0;
					for (int to = 0; to < num_node; to++)
					{
						odflow_out.width(5); odflow_out << right << node_id[to] << " :";
						odflow_out.width(9); odflow_out << right << setiosflags(ios::fixed) << setprecision(1) << od_flow_t[t][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();
			}
		}
		else
		{
			for (int t = 0; t < num_tod; t++)
			{
				string argv_series = argv[t + 7 + num_tod + num_tod + num_tod];
				int q = argv_series.find(".", 0);
				argv_series.erase(q, 4);
				stringstream argv_stream;
				argv_stream << scenario + 1;
				argv_series += "_";
				argv_series += argv_stream.str();
				argv_series += ".txt";
				const char *argv_series_c = argv_series.c_str();
				ofstream odflow_out(argv_series_c);
				odflow_out << "<NUMBER OF ZONES> " << node_id[num_node - 1] << "\n";
				odflow_out << "<TOTAL OD FLOW> " << total_od_flow_t[t] << "\n";
				odflow_out << "<END OF METADATA>" << "\n\n\n";
				for (int fr = 0; fr < num_node; fr++)
				{
					odflow_out << "Origin" << " " << "\t" << node_id[fr] << " \n";
					int count = 0;
					for (int to = 0; to < num_node; to++)
					{
						odflow_out.width(5); odflow_out << right << node_id[to] << " :";
						odflow_out.width(9); odflow_out << right << setiosflags(ios::fixed) << setprecision(1) << od_flow_t[t][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();
			}
		}
	}

	/*
	if (!argv[6 + num_tod + num_tod + num_tod + num_tod + num_tod])
	{
		argv[6 + num_tod + num_tod + num_tod + num_tod + num_tod] = "test_OD_pair_changes.txt";
		endtime = clock() / CLOCKS_PER_SEC;
	}
	ofstream odpair_out(argv[6 + num_tod + num_tod + num_tod + num_tod + num_tod]);
	for (int t = 0; t < num_tod; t++)
	{
		for (int c = 0; c < num_class; c++)
		{
			for (int m = 0; m < num_mode; m++)
			{
				for (int n1 = 0; n1 < num_node; n1++)
				{
					for (int n2 = 0; n2 < num_node; n2++)
					{
						if (base_od_flow_t_c_m[t][c][m][n1][n2] > 0)
						{
							odpair_out << t + 1 << "\t" << c + 1 << "\t" << m + 1 << "\t" << n1 + 1 << "\t" << n2 + 1 << "\t" << duration_t[t] << "\t"
								<< od_flow_t_c_m_old[t][c][m][n1][n2] << "\t" << od_flow_t_c_m[t][c][m][n1][n2] << "\t"
								<< od_cost_t_c_m_old[t][c][m][n1][n2] << "\t" << od_cost_t_c_m[t][c][m][n1][n2] << "\t" << ecu_m[m] << "\n";
						}
					}
				}
			}
		}
	}
	odpair_out.close();
	*/



	cout << "\n";
	cout << "   Traveler surplus change: " << user_surplus << " dollars" << "\n";
	cout << "   Traveler travel cost: " << travel_cost << " dollars" << "\n";
	cout << "   Traveler travel time: " << travel_time << " hours" << "\n";
	cout << "   Vehicle hours traveled: " << vht << " hours" << "\n";
	cout << "   Vehicle miles traveled: " << vmt << " miles" << "\n";
	cout << "\n";
//	for (int t = 0; t < num_tod; t++) {
//		cout << " Total trips during TOD " << t+1 << ": " << total_trip_rate_series[c][t] << "\n";
//	}
	endtime = clock() / CLOCKS_PER_SEC;
	cputime = endtime - starttime;
	cout << "   Clock" << "\t\t" << endtime << " : Finishing" << "\n";
	cout << "   CPU time elapsed: " << cputime << " sec" << "\n";
	cout << "\n";

	}

	// step 7 >> write the statistical result into the output file

	if (num_scenarios > 1)
	{
//		chdir(original_directory);
		double user_surplus_mean, travel_cost_mean, travel_time_mean, vht_mean, vmt_mean;
		double user_surplus_std, travel_cost_std, travel_time_std, vht_std, vmt_std;
		double user_surplus_cov, travel_cost_cov, travel_time_cov, vht_cov, vmt_cov;
		double user_surplus_sum = 0, travel_cost_sum = 0, travel_time_sum = 0, vht_sum = 0, vmt_sum = 0;
		double user_surplus_sum_square = 0, travel_cost_sum_square = 0, travel_time_sum_square = 0, vht_sum_square = 0, vmt_sum_square = 0;
		for (k = 0; k < num_scenarios; k++)
		{
			user_surplus_sum += user_surplus_series[k];
			travel_cost_sum += travel_cost_series[k];
			travel_time_sum += travel_time_series[k];
			vht_sum += vht_series[k];
			vmt_sum += vmt_series[k];
		}
		user_surplus_mean = user_surplus_sum / num_scenarios;
		travel_cost_mean = travel_cost_sum / num_scenarios;
		travel_time_mean = travel_time_sum / num_scenarios;
		vht_mean = vht_sum / num_scenarios;
		vmt_mean = vmt_sum / num_scenarios;
		for (k = 0; k < num_scenarios; k++)
		{
			user_surplus_sum_square += pow(user_surplus_series[k] - user_surplus_mean, 2);
			travel_cost_sum_square += pow(travel_cost_series[k] - travel_cost_mean, 2);
			travel_time_sum_square += pow(travel_time_series[k] - travel_time_mean, 2);
			vht_sum_square += pow(vht_series[k] - vht_mean, 2);
			vmt_sum_square += pow(vmt_series[k] - vmt_mean, 2);
		}
		user_surplus_std = pow(user_surplus_sum_square / num_scenarios, 0.5);
		travel_cost_std = pow(travel_cost_sum_square / num_scenarios, 0.5);
		travel_time_std = pow(travel_time_sum_square / num_scenarios, 0.5);
		vht_std = pow(vht_sum_square / num_scenarios, 0.5);
		vmt_std = pow(vmt_sum_square / num_scenarios, 0.5);
		user_surplus_cov = user_surplus_std / user_surplus_mean;
		travel_cost_cov = travel_cost_std / travel_cost_mean;
		travel_time_cov = travel_time_std / travel_time_mean;
		vht_cov = vht_std / vht_mean;
		vmt_cov = vmt_std / vmt_mean;
		ofstream summary_out(argv[6 + num_tod + num_tod + num_tod]);
		summary_out << "<NUMBER OF ZONES> " << num_node << "\n";
		summary_out << "<NUMBER OF LINKS> " << num_link << "\n";
		summary_out << "<NUMBER OF SCENARIOS> " << num_scenarios << "\n";
		summary_out << "<END OF METADATA>" << "\n\n\n";
		summary_out << "Traveler surplus change: " << user_surplus_mean << " dollars (mean); " << user_surplus_std << " dollars (std); " << user_surplus_cov << " (cov)" << "\n";
		summary_out << "Traveler travel cost: " << travel_cost_mean << " dollars (mean); " << travel_cost_std << " dollars (std); " << travel_cost_cov << " (cov)" << "\n";
		summary_out << "Traveler travel time: " << travel_time_mean << " hours (mean); " << travel_time_std << " hours (std); " << travel_time_cov << " (cov)" << "\n";
		summary_out << "Vehicle hours traveled: " << vht_mean << " hours (mean); " << vht_std << " hours (std); " << vht_cov << " (cov)" << "\n";
		summary_out << "Vehicle miles traveled: " << vmt_mean << " miles (mean); " << vmt_std << " miles (std); " << vmt_cov << " (cov)" << "\n";
		summary_out << "\n";
		//summary_out << "Trips in time period " << t << ": " << total_trip_rate_series[t] << "\n";
		for (int c = 0; c < num_class; c++) {
			for (int t = 0; t < num_tod; t++) {
				summary_out << " Total trips for class " << c+1 << " during TOD " << t+1 << ": " << total_trip_rate_series[c][t] << "\n";
			}
			summary_out << "\n";
		}
		summary_out.close();
	}

	return 0;

}

// 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];
				}
			}
		}
	}
}

// generate lognormal random variables by the Box-Muller method
double lognormal_generator(double mean_l, double std_l)
{
	srand(time(NULL));
	//double random_num0 = (double) rand() / RAND_MAX;
	double random_num1 = (double) rand() / RAND_MAX;
	double random_num2 = (double) rand() / RAND_MAX;
	double normal_num = pow(-2 * log(random_num1), 0.5) * sin(2 * PI * random_num2);
	double lognormal_num;
	if (std_l == 0)
	{
		lognormal_num = mean_l;
	}
	else
	{
		double mean = log(pow(mean_l, 2) / pow(pow(mean_l, 2) + pow(std_l, 2), 0.5));
		double std = pow(log((pow(mean_l, 2) + pow(std_l, 2)) / pow(mean_l, 2)), 0.5);
		lognormal_num = exp(mean + std * normal_num);
	}

	return lognormal_num;
}

void readTransitFile(int ***routeMatrix, double **routeHeadway, int *numRoutes, char *transitFileName, int numArcs) {

	FILE *transitFile = fopen(transitFileName, "r");
	if (transitFile == NULL) { *numRoutes = 0; return; }
	char fullLine[999];
	fgets(fullLine, 999, transitFile);
	sscanf(fullLine, "<NUMBER OF ROUTES> %d", numRoutes);
	fgets(fullLine, 999, transitFile); // End of metadata
	fgets(fullLine, 999, transitFile);
	int r, check, arc, arcno;

	*routeMatrix = (int **)malloc(*numRoutes * sizeof(int *));
	for (r = 0; r < *numRoutes; r++) (*routeMatrix)[r] = (int *)malloc(numArcs * sizeof(int));
	*routeHeadway = (double *)malloc(*numRoutes * sizeof(double));

	for (r = 0; r < *numRoutes; r++) {
		fgets(fullLine, 999, transitFile);
		sscanf(fullLine, "Route %d : Headway %lf", &check, &((*routeHeadway)[r]));
		//if (check != r) exit(EXIT_FAILURE); // TODO : Could include route # index for reporting ridership?

		//fgets(fullLine, 999, transitFile);
		//routeHeadway[r] = atof(fullLine);

		arcno = 0;
		do {
			fgets(fullLine, 999, transitFile);
			arc = atoi(fullLine);
			if (arcno >= numArcs || arc > numArcs || (arc < 0 && arc != END_OF_ROUTE )) {
				printf("Bad transit route %d contains out-of-range arc %d.\n", r+1, arc);
				exit(EXIT_FAILURE);
			}
			(*routeMatrix)[r][arcno++] = (arc == END_OF_ROUTE) ? arc : arc - 1;
		} while (arc != END_OF_ROUTE);
		fgets(fullLine, 999, transitFile);
	}
}

void transitZoneZoneCosts(int num_node, long int *node_id, double **od_cost, int **routeMatrix, double *routeHeadway, int num_routes, long int **link_id, double *linkTravelTime) {


	// Allocate
	int o, d, r, a, b, i, j, t;
	double **old_od_cost = (double **)malloc(num_node * sizeof(double *));
	for (o = 0; o < num_node; o++) old_od_cost[o] = (double *)malloc(num_node * sizeof(double));
	double ***route_cost = (double ***)malloc(num_node * sizeof(double **));
	for (o = 0; o < num_node; o++) {
		route_cost[o] = (double **)malloc(num_node * sizeof(double *));
		for (d = 0; d < num_node; d++) {
			route_cost[o][d] = (double *)malloc(num_routes * sizeof(double));
		}
	}

	// Initialize
	for (o = 0; o < num_node; o++) {
		for (d = 0; d < num_node; d++) {
			od_cost[o][d] = (o == d) ? 0 : positive_infinity;
			old_od_cost[o][d] = positive_infinity;
			for (r = 0; r < num_routes; r++) {
				route_cost[o][d][r] = positive_infinity;
			}
		}
	}
	if (num_routes == 0) goto done;

	// Calculate zone-zone travel times for each route
	for (r = 0; r < num_routes; r++) {
		a = -1;
		while (routeMatrix[r][++a] != END_OF_ROUTE) {
			i = link_id[routeMatrix[r][a]][0] - 1;
			route_cost[i][i][r] = routeHeadway[r] / 2;
			b = a;
			while (routeMatrix[r][b] != END_OF_ROUTE) {
				j = link_id[routeMatrix[r][b]][1] - 1;
				route_cost[i][j][r] = route_cost[i][link_id[routeMatrix[r][b]][0] - 1][r] + linkTravelTime[b];
				b++;
			}
		}
	}

	// Iteration

	for (t = 0; t <= MAX_TRANSFERS; t++) {
		for (o = 0; o < num_node; o++) {
			for (d = 0; d < num_node; d++) {
				old_od_cost[o][d] = od_cost[o][d];
			}
		}
		for (o = 0; o < num_node; o++) {
			for (d = 0; d < num_node; d++) {
				for (r = 0; r < num_routes; r++) {
					for (a = 0; a < num_node; a++) {
//						if (od_cost[o][d] > old_od_cost[o][a] + route_cost[a][d][r]) fprintf(myLogFile, "Set trasnfer %d stop time from %d to %d by transferring at %d to %d: new time %f=%f+%f\n", t, o, d, a, r, old_od_cost[o][a] + route_cost[a][d][r], old_od_cost[o][a], route_cost[a][d][r]);
						od_cost[o][d] = min(od_cost[o][d], old_od_cost[o][a] + route_cost[a][d][r]);
					}
				}
			}
		}
	}

done:
	// Here call frees
	for (o = 0; o < num_node; o++) {
		for (d = 0; d < num_node; d++) {
			free(route_cost[o][d]);
		}
		free(route_cost[o]);
		free(old_od_cost[o]);
	}
	free(route_cost);
	free(old_od_cost);
}