main.py

"""
HARMONY_LUTI_ATH
main.py

May - August 2020
Author: Eleni Kalantzi, Msc Smart Cities and Urban Analytics, Centre for Advanced Spatial Analysis, University College London
Supervisor: Fulvio D. Lopane, Research Fellow, Centre for Advanced Spatial Analysis, University College London
https://www.casa.ucl.ac.uk

Repository: https://github.com/elenikalantzi/Dissertation---CASA0010

Developed from Richard Milton's QUANT_RAMP And HARMONY_LUTI_OXFORDSHIRE
"""
import os
import time
import pandas as pd
from geojson import dump

from globals import *
from utils import loadQUANTMatrix, loadMatrix, saveMatrix
from analytics import graphProbabilities, flowArrowsGeoJSON
from quantlhmodel import QUANTLHModel
import numpy as np

################################################################################
# Initialisation                                                               #
################################################################################

# NOTE: this section provides the base data for the models that come later. This
# will only be run on the first run of the program to assemble all the tables
# required from the original sources. After that, if the file exists in the
# directory, then nothing new is created and this section is effectively
# skipped, up until the model run section at the end.

# make a model-runs dir if we need it
if not os.path.exists(modelRunsDir):
    os.makedirs(modelRunsDir)

################################################################################
# Journey to work Model                                                        #
################################################################################

"""
runJourneyToWorkModel
Same as runRetailModel, except that now Origins: workplaces, Destinations: households' population
"""
# Journey to work model with households (HH) floorspace as attractor
def runJourneyToWorkModel(Scenario='2019', Beta_calibrated=None):
    print("Running Journey to Work ", Scenario, " model.")
    start = time.perf_counter()
    # Singly constrained model:
    # We conserve the number of jobs and predict the working population residing in Athens zones
    # journeys to work generated by jobs
    # Origins: workplaces
    # Destinations: Zones' households
    # Attractor: floorspace of housing

    """
                    Journey to work       |   Retail model
     Origins:       workplaces            |   households
     Destinations:  households            |   supermarkets
     conserved:     jobs                  |   income
     predicted:     population of zones   |   expenditure @ supermarkets
     attractor:     HH floorspace density |   supermarket floorspace
    """

    # load jobs data for residential zones
    dfEi = pd.read_csv(data_employment, usecols=['zone','employment'], index_col='zone')
    dfEi.astype({'employment': 'int64'})

    # df_pop = pd.read_csv(data_census_pop, usecols=['zone','pop_tot'], index_col='zone')

    df_floorspace = pd.read_csv(HH_floorspace, usecols=['zone','hh_floorspace_density'], index_col='zone')
    # Need to substitute 0 values in floorspace dataframe with very low values (e.g. 1) to avoid division by zero:
    df_floorspace.replace(0, 1, inplace=True)

    # if we are running a scenario, update the Zones with the new attractors and the numbers of jobs
    if Scenario == '2019':
        # This is the base scenario, so we don't have to modify the attractors/jobs
        pass

    elif Scenario == 'Elliniko_2030':
        dfEi_30 = pd.read_csv(data_employment_2030, usecols=['zone', 'employment'], index_col='zone')

        # df_pop = pd.read_csv(data_census_pop, usecols=['zone', 'pop_tot'], index_col='zone')

        df_floorspace = pd.read_csv(HH_floorspace, usecols=['zone', 'hh_floorspace_density'], index_col='zone')
        # Need to substitute 0 values in floorspace dataframe with very low values (e.g. 1) to avoid division by zero:
        df_floorspace.replace(0, 1, inplace=True)

    elif Scenario == 'Elliniko_2045':
        dfEi_45 = pd.read_csv(data_employment_2045, usecols=['zone', 'employment'], index_col='zone')
        dfEi_45.astype({'employment': 'int64'})

        # Scenario_pop_table = pd.read_csv(pop_2045, usecols=['zone', 'pop_tot'], index_col='zone')

        Scenario_floorspace = pd.read_csv(HH_floorspace_2045, usecols=['zone', 'hh_floorspace_density'], index_col='zone')
        # Need to substitute 0 values in floorspace dataframe with very low values (e.g. 1) to avoid division by zero:
        Scenario_floorspace.replace(0, 1, inplace=True)

    # threshold = 2.7 *(1./1265)  # The threshold below which to ignore low probability trips - for geojson flows files (lines)

    if Scenario == '2019':
        # Load observed data for model calibration:
        # OD_Pu, OD_Pr = QUANTLHModel.loadODData()

        # Use cij as cost matrix (zone to zone)
        m, n = cij_pu.shape
        model = QUANTLHModel(m, n)

        # OBS Cbar public is 47 minutes, according to Public transit facts & statistics for Athens | Moovit Public Transit Index (moovitapp.com)
        # OBS Cbar private is 37 minutes, according to Traffic in Athens (numbeo.com)
        model.setObsCbar(47, 37)

        # model.setODMatrix(OD_Pu, OD_Pr)
        model.setAttractorsAj(df_floorspace, 'zone', 'hh_floorspace_density')
        model.setPopulationEi(dfEi, 'zone', 'employment')
        model.setCostMatrixCij(cij_pu, cij_pr)

        Tij, beta_k, cbar_k_pred, cbar_k_obs = model.runTwoModes()

        #NOTE: these are saved later as csvs, but not with an easy to read formatter
        np.savetxt("debug_Tij_public_2019.txt", Tij[0], delimiter=",", fmt="%i")
        np.savetxt("debug_Tij_private_2019.txt", Tij[1], delimiter=",", fmt="%i")
        #now an Oj Dj table
        # DfEi is employment - really hope these come out in the right order
        dfEi['DjPred_pu'] = Tij[0].sum(axis=1)
        dfEi['DjPred_pr'] = Tij[1].sum(axis=1)
        dfEi['DjPred'] = Tij[0].sum(axis=1) + Tij[1].sum(axis=1)
        dfEi['OiPred_pu'] = Tij[0].sum(axis=0)
        dfEi['OiPred_pr'] = Tij[1].sum(axis=0)
        dfEi['OiPred'] = Tij[0].sum(axis=0) + Tij[1].sum(axis=0)
        dfEi.to_csv(os.path.join(modelRunsDir,"EiDjOi_2019.csv"))

        #print("Computing probabilities...")
        # Compute the probability of a flow from a zone to any (i.e. all) of the possible point zones.
        jobs_probTij = model.computeProbabilities2modes(Tij)

        # Save output matrices
        print("Saving output matrices...")

        # Jobs accessibility:
        # Jobs_accessibility = number of jobs / d
        # d = average TT travelled by workers = ( cij * Tij ) / Tij_tot
        d_pu = (cij_pu * Tij[0]).sum(axis=1) / (Tij[0].sum(axis=1) + np.ones(len(Tij[0].sum(axis=1)))) # To avoid division by zero
        d_pr = (cij_pr * Tij[1]).sum(axis=1) / (Tij[0].sum(axis=1) + np.ones(len(Tij[0].sum(axis=1))))
        Jobs_accessibility_pu = Tij[0].sum(axis=1) / (d_pu + np.ones(len(Tij[0].sum(axis=1))))
        Jobs_accessibility_pr = Tij[0].sum(axis=1) / (d_pr + np.ones(len(Tij[0].sum(axis=1))))
        # Normalise the two vectors:
        Jobs_accessibility_pu_norm = Jobs_accessibility_pu / np.sqrt(np.sum(Jobs_accessibility_pu ** 2))
        Jobs_accessibility_pr_norm = Jobs_accessibility_pr / np.sqrt(np.sum(Jobs_accessibility_pr ** 2))
        Jobs_accessibility_df = pd.DataFrame(
            {'Jobs_accessibility_pu': Jobs_accessibility_pu_norm, 'Jobs_accessibility_pr': Jobs_accessibility_pr_norm})
        Jobs_accessibility_df.to_csv(os.path.join(modelRunsDir, "Jobs_accessibility_2019.csv"))

        # Probabilities:
        saveMatrix(jobs_probTij[0], data_jobs_probTij_public_2019)
        saveMatrix(jobs_probTij[1], data_jobs_probTij_private_2019)
        np.savetxt(data_jobs_probTij_public_2019_csv, jobs_probTij[0], delimiter=",")
        np.savetxt(data_jobs_probTij_private_2019_csv, jobs_probTij[1], delimiter=",")

        # People flows
        saveMatrix(Tij[0], data_jobs_Tij_pu_2019)
        saveMatrix(Tij[1], data_jobs_Tij_pr_2019)
        np.savetxt(data_jobs_Tij_pu_2019_csv, Tij[0], delimiter=",")
        np.savetxt(data_jobs_Tij_pr_2019_csv, Tij[1], delimiter=",")

        '''
        # Geojson flows files - lines
        dfPointsPopulation_2019 = pd.read_csv(os.path.join(datadir_ex, 'zones_data_coordinates.csv'),
                                              usecols=["zonei", "zone", "pop_tot_2019", "employment_2019",
                                                       "Greek_Grid_east", "Greek_Grid_north"])  # create df from csv
        geojson_flows_2019_pr = graphProbabilities(threshold, dfPointsPopulation_2019, jobs_probTij[1])
        with open(os.path.join(modelRunsDir, 'flows_2019_pr.geojson'), 'w') as f:
            dump(geojson_flows_2019_pr, f)

        geojson_flows_2019_pu = graphProbabilities(threshold, dfPointsPopulation_2019, jobs_probTij[0])
        with open(os.path.join(modelRunsDir, 'flows_2019_pu.geojson'), 'w') as f:
            dump(geojson_flows_2019_pu, f)
        '''

        # Geojson flows files - arrows
        # I need my own zone codes file containing the zonei and greek grid indexes as
        # ZoneCodes_ATH does not contain the information
        flow_zonecodes = pd.read_csv('external-data/zones_data_coordinates.csv')
        flow_pu = flowArrowsGeoJSON(Tij[0], flow_zonecodes)
        with open(os.path.join(modelRunsDir, 'flows_2019_pu.geojson'), 'w') as f:
            dump(flow_pu, f)
        flow_pr = flowArrowsGeoJSON(Tij[1], flow_zonecodes)
        with open(os.path.join(modelRunsDir, 'flows_2019_pr.geojson'), 'w') as f:
            dump(flow_pr, f)

        # cbar = model.computeCBar(Tij, cij_pu)
        print("JtW model", Scenario, "observed cbar [public, private] = ", cbar_k_obs)
        print("JtW model", Scenario, "predicted cbar [public, private] = ", cbar_k_pred)
        print("JtW model", Scenario, "beta [public, private] = ", beta_k)

        # Calculate predicted population
        DjPred = np.zeros(n)
        for k in range(len(Tij)):
            DjPred += Tij[k].sum(axis=1)
        # Create a dataframe with Zone and people count
        DjPred = pd.DataFrame(DjPred, columns=['population'])
        DjPred['zone'] = zonecodes_ATH_list

        end = time.perf_counter()
        # print("Journey to work model", Scenario," run elapsed time (secs)=", end - start)
        print("Journey to work model run elapsed time (minutes)=", (end - start) / 60)
        print()

        return beta_k, DjPred

    elif Scenario == 'Elliniko_2030':
        # Use cij as cost matrix
        m, n = cij_pu.shape
        model = QUANTLHModel(m, n)
        model.setAttractorsAj(df_floorspace, 'zone', 'hh_floorspace_density')
        model.setPopulationEi(dfEi_30, 'zone', 'employment')
        model.setCostMatrixCij(cij_pu, cij_pr)

        Tij, cbar_k = model.run2modes_NoCalibration(Beta_calibrated)
        # Compute the probability of a flow from a zone to any (i.e. all) of the possible point zones.
        jobs_probTij = model.computeProbabilities2modes(Tij)

        # NOTE: these are saved later as csvs, but not with an easy to read formatter
        np.savetxt("debug_Tij_public_2030.txt", Tij[0], delimiter=",", fmt="%i")
        np.savetxt("debug_Tij_private_2030.txt", Tij[1], delimiter=",", fmt="%i")
        # now an Oj Dj table
        # DfEi is employment - really hope these come out in the right order
        dfEi['DjPred_pu'] = Tij[0].sum(axis=1)
        dfEi['DjPred_pr'] = Tij[1].sum(axis=1)
        dfEi['DjPred'] = Tij[0].sum(axis=1) + Tij[1].sum(axis=1)
        dfEi['OiPred_pu'] = Tij[0].sum(axis=0)
        dfEi['OiPred_pr'] = Tij[1].sum(axis=0)
        dfEi['OiPred'] = Tij[0].sum(axis=0) + Tij[1].sum(axis=0)
        dfEi.to_csv(os.path.join(modelRunsDir,"EiDjOi_2030.csv"))

        # Save output matrices
        print("Saving output matrices...")

        # Jobs accessibility:
        # Jobs_accessibility = number of jobs / d
        # d = average TT travelled by workers = ( cij * Tij ) / Tij_tot
        d_pu = (cij_pu * Tij[0]).sum(axis=1) / (
                    Tij[0].sum(axis=1) + np.ones(len(Tij[0].sum(axis=1))))  # To avoid division by zero
        d_pr = (cij_pr * Tij[1]).sum(axis=1) / (Tij[0].sum(axis=1) + np.ones(len(Tij[0].sum(axis=1))))
        Jobs_accessibility_pu = Tij[0].sum(axis=1) / (d_pu + np.ones(len(Tij[0].sum(axis=1))))
        Jobs_accessibility_pr = Tij[0].sum(axis=1) / (d_pr + np.ones(len(Tij[0].sum(axis=1))))
        # Normalise the two vectors:
        Jobs_accessibility_pu_norm = Jobs_accessibility_pu / np.sqrt(np.sum(Jobs_accessibility_pu ** 2))
        Jobs_accessibility_pr_norm = Jobs_accessibility_pr / np.sqrt(np.sum(Jobs_accessibility_pr ** 2))
        Jobs_accessibility_df = pd.DataFrame(
            {'Jobs_accessibility_pu': Jobs_accessibility_pu_norm, 'Jobs_accessibility_pr': Jobs_accessibility_pr_norm})
        Jobs_accessibility_df.to_csv(os.path.join(modelRunsDir, "Jobs_accessibility_2030.csv"))

        # Probabilities:
        saveMatrix(jobs_probTij[0], data_jobs_probTij_public_2030)
        saveMatrix(jobs_probTij[1], data_jobs_probTij_private_2030)
        np.savetxt(data_jobs_probTij_public_2030_csv, Tij[0], delimiter=",")
        np.savetxt(data_jobs_probTij_private_2030_csv, Tij[1], delimiter=",")

        # People flows
        saveMatrix(Tij[0], data_jobs_Tij_pu_2030)
        saveMatrix(Tij[1], data_jobs_Tij_pr_2030)
        np.savetxt(data_jobs_Tij_pu_2030_csv, Tij[0], delimiter=",")
        np.savetxt(data_jobs_Tij_pr_2030_csv, Tij[1], delimiter=",")

        '''
        # Geojson files - lines
        dfPointsPopulation_2030 = pd.read_csv(os.path.join(datadir_ex, 'zones_data_coordinates.csv'),
                                              usecols=["zonei", "zone", "pop_tot_2030", "employment_2030",
                                                       "Greek_Grid_east", "Greek_Grid_north"])  # create df from csv
        geojson_flows_2030_pr = graphProbabilities(threshold, dfPointsPopulation_2030, jobs_probTij[1])
        with open(os.path.join(modelRunsDir, 'flows_2030_pr.geojson'), 'w') as f:
            dump(geojson_flows_2030_pr, f)

        geojson_flows_2030_pu = graphProbabilities(threshold, dfPointsPopulation_2030, jobs_probTij[0])
        with open(os.path.join(modelRunsDir, 'flows_2030_pu.geojson'), 'w') as f:
            dump(geojson_flows_2030_pu, f)
        '''

        # Geojson flows files - arrows
        # I need my own zone codes file containing the zonei and greek grid indexes as
        # ZoneCodes_ATH does not contain the information
        flow_zonecodes = pd.read_csv('external-data/zones_data_coordinates.csv')
        flow_pu = flowArrowsGeoJSON(Tij[0], flow_zonecodes)
        with open(os.path.join(modelRunsDir, 'flows_2030_pu.geojson'), 'w') as f:
            dump(flow_pu, f)
        flow_pr = flowArrowsGeoJSON(Tij[1], flow_zonecodes)
        with open(os.path.join(modelRunsDir, 'flows_2030_pr.geojson'), 'w') as f:
            dump(flow_pr, f)

        print("JtW model", Scenario, " cbar [public, private] = ", cbar_k)

        DjPred = np.zeros(n)
        for k in range(len(Tij)):
            DjPred += Tij[k].sum(axis=1)
        # Create a dataframe with Zone and people count
        DjPred = pd.DataFrame(DjPred, columns=['population'])
        DjPred['zone'] = zonecodes_ATH_list

        end = time.perf_counter()
        # print("Journey to work model run elapsed time (secs) =", end - start)
        print("Journey to work model run elapsed time (minutes) =", (end - start)/60)
        print()

        return DjPred

    elif Scenario == 'Elliniko_2045':
        # Use cij as cost matrix (zone to zone)
        m, n = cij_pu.shape
        model = QUANTLHModel(m, n)
        model.setAttractorsAj(Scenario_floorspace, 'zone', 'hh_floorspace_density')
        model.setPopulationEi(dfEi_45, 'zone', 'employment')
        model.setCostMatrixCij(cij_pu, cij_pr)

        Tij, cbar_k = model.run2modes_NoCalibration(Beta_calibrated)
        # Compute the probability of a flow from a zone to any (i.e. all) of the possible point zones.
        jobs_probTij = model.computeProbabilities2modes(Tij)

        # NOTE: these are saved later as csvs, but not with an easy to read formatter
        np.savetxt("debug_Tij_public_2045.txt", Tij[0], delimiter=",", fmt="%i")
        np.savetxt("debug_Tij_private_2045.txt", Tij[1], delimiter=",", fmt="%i")
        # now an Oj Dj table
        # DfEi is employment - really hope these come out in the right order
        dfEi['DjPred_pu'] = Tij[0].sum(axis=1)
        dfEi['DjPred_pr'] = Tij[1].sum(axis=1)
        dfEi['DjPred'] = Tij[0].sum(axis=1) + Tij[1].sum(axis=1)
        dfEi['OiPred_pu'] = Tij[0].sum(axis=0)
        dfEi['OiPred_pr'] = Tij[1].sum(axis=0)
        dfEi['OiPred'] = Tij[0].sum(axis=0) + Tij[1].sum(axis=0)
        dfEi.to_csv(os.path.join(modelRunsDir,"EiDjOi_2045.csv"))

        # Save output matrices
        print("Saving output matrices...")

        # Jobs accessibility:
        # Jobs_accessibility = number of jobs / d
        # d = average TT travelled by workers = ( cij * Tij ) / Tij_tot
        d_pu = (cij_pu * Tij[0]).sum(axis=1) / (
                    Tij[0].sum(axis=1) + np.ones(len(Tij[0].sum(axis=1))))  # To avoid division by zero
        d_pr = (cij_pr * Tij[1]).sum(axis=1) / (Tij[0].sum(axis=1) + np.ones(len(Tij[0].sum(axis=1))))
        Jobs_accessibility_pu = Tij[0].sum(axis=1) / (d_pu + np.ones(len(Tij[0].sum(axis=1))))
        Jobs_accessibility_pr = Tij[0].sum(axis=1) / (d_pr + np.ones(len(Tij[0].sum(axis=1))))
        # Normalise the two vectors:
        Jobs_accessibility_pu_norm = Jobs_accessibility_pu / np.sqrt(np.sum(Jobs_accessibility_pu ** 2))
        Jobs_accessibility_pr_norm = Jobs_accessibility_pr / np.sqrt(np.sum(Jobs_accessibility_pr ** 2))
        Jobs_accessibility_df = pd.DataFrame(
            {'Jobs_accessibility_pu': Jobs_accessibility_pu_norm, 'Jobs_accessibility_pr': Jobs_accessibility_pr_norm})
        Jobs_accessibility_df.to_csv(os.path.join(modelRunsDir, "Jobs_accessibility_2045.csv"))

        # Probabilities:
        saveMatrix(jobs_probTij[0], data_jobs_probTij_public_2045)
        saveMatrix(jobs_probTij[1], data_jobs_probTij_private_2045)
        np.savetxt(data_jobs_probTij_public_2045_csv, Tij[0], delimiter=",")
        np.savetxt(data_jobs_probTij_private_2045_csv, Tij[1], delimiter=",")

        # People flows:
        saveMatrix(Tij[0], data_jobs_Tij_pu_2045)
        saveMatrix(Tij[1], data_jobs_Tij_pr_2045)
        np.savetxt(data_jobs_Tij_pu_2045_csv, Tij[0], delimiter=",")
        np.savetxt(data_jobs_Tij_pr_2045_csv, Tij[1], delimiter=",")

        '''
        # Geojson  files - lines
        dfPointsPopulation_2045 = pd.read_csv(os.path.join(datadir_ex, 'zones_data_coordinates.csv'),
                                              usecols=["zonei", "zone", "pop_tot_2045", "employment_2045",
                                                       "Greek_Grid_east", "Greek_Grid_north"])  # create df from csv
        geojson_flows_2045_pr = graphProbabilities(threshold, dfPointsPopulation_2045, jobs_probTij[1])
        with open(os.path.join(modelRunsDir, 'flows_2045_pr.geojson'), 'w') as f:
            dump(geojson_flows_2045_pr, f)

        geojson_flows_2045_pu = graphProbabilities(threshold, dfPointsPopulation_2045, jobs_probTij[0])
        with open(os.path.join(modelRunsDir, 'flows_2045_pu.geojson'), 'w') as f:
            dump(geojson_flows_2045_pu, f)
        '''

        # Geojson flows files - arrows
        # I need my own zone codes file containing the zonei and greek grid indexes as
        # ZoneCodes_ATH does not contain the information
        flow_zonecodes = pd.read_csv('external-data/zones_data_coordinates.csv')
        flow_pu = flowArrowsGeoJSON(Tij[0], flow_zonecodes)
        with open(os.path.join(modelRunsDir, 'flows_2045_pu.geojson'), 'w') as f:
            dump(flow_pu, f)
        flow_pr = flowArrowsGeoJSON(Tij[1], flow_zonecodes)
        with open(os.path.join(modelRunsDir, 'flows_2045_pr.geojson'), 'w') as f:
            dump(flow_pr, f)

        print("JtW model", Scenario, " cbar [public, private] = ", cbar_k)

        # Calculate predicted population
        DjPred = np.zeros(n)
        for k in range(len(Tij)):
            DjPred += Tij[k].sum(axis=1)
        # Create a dataframe with Zone and people count
        DjPred = pd.DataFrame(DjPred, columns=['population'])
        DjPred['zone'] = zonecodes_ATH_list

        end = time.perf_counter()
        # print("Journey to work model run elapsed time (secs)=", end - start)
        print("Journey to work model run elapsed time (mins)=", (end - start) / 60)
        print()

        return DjPred

################################################################################
def runEllinikoScenarios():
    # First run the base model to calibrate it with 2011 observed trip data:
    # Run Journey to work model:
    beta_2019, DjPred_JtW_2019 = runJourneyToWorkModel()

    # Create population variables for 2030 and 2045:
    # Elliniko_2030 = pd.read_csv(data_census_pop, usecols=['zone', 'pop_tot'])
    # Elliniko_2045 = pd.read_csv(pop_2045, usecols=['zone', 'pop_tot'])

    # Now run the JtW model with 2011 beta and 2019 pop (without calibration this time)
    DjPred_JtW_2030 = runJourneyToWorkModel('Elliniko_2030', beta_2019)
    DjPred_JtW_2045 = runJourneyToWorkModel('Elliniko_2045', beta_2019)

################################################################################
# Now on to the model run section
################################################################################

#Zone Codes for Athens
print()
print("Importing ATH cij matrices")

zonecodes_ATH = pd.read_csv(os.path.join(modelRunsDir,ZoneCodesFilename))
zonecodes_ATH.set_index('zone')
zonecodes_ATH_list = zonecodes_ATH['zone'].tolist()
# print('zonecodes_ATH:')
# print(zonecodes_ATH)
# print()
#_____________________________________________________________________________________
# Import cij data

# Fisrt, import intrazone distances
intrazone_dist_df = pd.read_csv(os.path.join(datadir_ex,intrazone_dist_csv), usecols=["Intrazone_dist"]) # create df from csv
intrazone_dist_list = intrazone_dist_df.Intrazone_dist.values.tolist() # save the intrazone distances into a list

# Import the csv cij private as Pandas DataFrame:
cij_pr_df = pd.read_csv(os.path.join(datadir_ex,CijPrivateMinFilename_csv), header=None)

# Convert the dataframe to a nupmy matrix:
cij_pr = cij_pr_df.to_numpy()
cij_pr[cij_pr < 0] = 120 # upper limit of two hours - change missing vlaues (set as -1) to very high numbers to simulate impossible routes
cij_pr[cij_pr < 1] = 1 #lower limit of 1 minute links

# Change values in the main diagonal from 0 to intrazone impedance:
av_speed_pr = 13 # define average speed in km/h - based on Athens vehicles
average_TT_pr = [] # average intrazonal travel time
for i in intrazone_dist_list:
    average_TT_pr.append((i/1000) / (av_speed_pr/60)) # result in minutes
np.fill_diagonal(cij_pr, average_TT_pr) # save the average TT list as the main diagonal of the cij matrix

# Print the dimensions of the matrix to check that everything is ok:
print('cij private shape: ', cij_pr.shape)
cij_pr[cij_pr<1] = 1 # repeat of above after leading diagonal zeros added
cij_pr = cij_pr * 2.5 # HACK! Make driving less attractive - are costs wrong? --> probably speeds from OASA not to be trusted?

# Import the csv cij public as Pandas DataFrame:
cij_pu_df = pd.read_csv(os.path.join(datadir_ex,CijPublicMinFilename_csv), header=None)

# Convert the dataframe to a nupmy matrix:
cij_pu = cij_pu_df.to_numpy()
cij_pu[cij_pu < 0] = 120 # upper limit of two hours - change missing vlaues (set as -1) to very high numbers to simulate impossible routes
cij_pu[cij_pu < 1] = 1 # lower limit of 1 minute links
# Add waiting time: 16min for Athen from Moovit data:
waiting_time_proportion = 16.0/47.0 # 16 = average waiting time; 47 = average travel time
cij_pu[cij_pu < 9999] *= (1 + waiting_time_proportion)

# Change values in the main diagonal from 0 to intrazone impedance:
av_speed_pu = 10 # define average speed in km/h - based on Athens buses
average_TT_pu = [] # average intrazonal travel time
for i in intrazone_dist_list:
    average_TT_pu.append((i/1000) / (av_speed_pu/60)) # result in minutes
np.fill_diagonal(cij_pu, average_TT_pu) # save the average TT list as the main diagonal of the cij matrix

# Export cij matrices for checking
# np.savetxt(os.path.join(modelRunsDir,'debug_cij_pr.csv'), cij_pr, delimiter=',', fmt='%i')
# np.savetxt(os.path.join(modelRunsDir,'debug_cij_pu.csv'), cij_pu, delimiter=',', fmt='%i')
#RWM there are zeros in the private one and 9999s in the public one...

# Print the dimensions of the matrix to check that everything is ok:
print('cij public shape: ', cij_pu.shape)
#_____________________________________________________________________________________
'''
# Import OD data
# Import the csv OD private as Pandas DataFrame:
OD_pr_df = pd.read_csv(os.path.join(datadir_ex,ODPrFilename_csv), header=None)
# Convert the dataframe to a nupmy matrix:
OD_pr = OD_pr_df.to_numpy()
print('OD private shape: ', OD_pr.shape)

# Import the csv OD public as Pandas DataFrame:
OD_pu_df = pd.read_csv(os.path.join(datadir_ex,ODPuFilename_csv), header=None)
# Convert the dataframe to a nupmy matrix:
OD_pu = OD_pu_df.to_numpy()
print('OD public shape: ', OD_pu.shape)
'''
#_____________________________________________________________________________________

# now run the relevant models to produce the outputs
runEllinikoScenarios()

################################################################################
# END OF MAIN PROGRAM                                                          #
################################################################################