#-------------------------------------------------------------------------------------------------------------------------------
# Name:        Basic-Data
# Purpose:
#
# Author:      Gaoy
#
# Created:     09/10/2015
# Copyright:   (c) Gaoy 2015
# Licence:     <your licence>
#-------------------------------------------------------------------------------------------------------------------------------

# This file consists of two parts. The first part includes basic parameters, such as station number, services number,
# link running times, station dwelling times, extra acc/dec times, headway, passenger arrival rate, OD matrix and so on.
# The second part simulates the passenger flow in the disruption period, which generates the numbers of passengers
# stranded in stations and in trains in the end of the disruption.

# Note that in Python,  array index starts from "0"!!!

import numpy as np
import scipy as sp
import xlwt


stationNum=13 #Station amount
M=5000.0 #A large value
M1=5000.0
M2=5000.0

START=0 #the beginning of the recovery is set as the "START"
Maxdt=60
Delta1=300
Delta2=45 # slack time in the model

W=900 #900 seconds or 15 minutes
Gap=0.001 #Gurobi Gap

SK1=1
SK2=2
SK3=3
# number of stations skipped by a service

xlength=5000 # This is the length of x axis, for plotting

H1=0 #disruption starting time: 7:00+H1*15
H2=2 #recovery starting time:  7:00+H2*15


Disruption=(H2-H1)*W #disruption period


trainNum=20 # number of trains
existingService=20 # there are 10 upstream existing services and 10 downstream existing services

serviceNum=35 # we consider (serviceNum-10)=25 services in each direction, that is, 10 existing services and 15 future services


sectionlen=[2631,1275,2366,1982,993,1728,1090,1355,2337,2301,2055,1281] # the length of each section

accExtra=25
decExtra=25 # extra acceleration/deceleration time

headway=90 # minimum headway
backturning=120 # minimum back-turning time
linkrun10=[190,108,157,135,90,114,103,104,166,150,140,102] # pure running times in upstream direction
linkrun20=[100,141,150,162,103,101,111,90,135,157,105,195] # pure running times in downstream direction

stationdt=[30,35,30,30,35,30,30,30,30,30,30,35,45]
stationdt.extend([40,35,30,30,30,35,30,30,35,30,30,35,30]) # minimum dwelling time in each station


N=stationNum # N is short for "stationNum"
Capacity=1200 # train capacity

#----------------------------------------------Generate passenger arrival rate-------------------------------------------

#------------------1----2---3---4---5---6---7---8---9--10--11--12---13 7:00-10:00
PsnVol0=np.array([[315,420,543,630,720,795,618,540,471,420,360,300,300],             #1
                  [495,600,750,732,810,930,939,780,663,540,444,390,402],             #2
                  [486,540,570,630,600,618,570,420,390,354,282,300,288],              #3
                  [1170,1440,1818,2190,2280,2187,1920,1812,1710,1470,1500,1200,972], #4
                  [756,900,1014,1080,1125,1170,1032,1020,930,780,720,612,510],       #5
                  [315,357,354,390,435,552,495,441,390,381,312,279,228],             #6
                  [273,330,393,480,549,663,510,438,378,339,303,255,216],             #7
                  [222,303,360,420,516,555,450,390,360,339,303,288,243],             #8
                  [957,1020,1170,1344,1416,1518,1440,1341,1170,1059,906,756,642],    #9
                  [2355,2556,2766,2862,2796,2835,2850,2715,2541,2196,1857,1590,1314],#10
                  [585,672,744,870,846,744,660,591,603,510,426,303,267],             #11
                  [180,222,249,171,147,216,264,159,138,240,168,96,90],               #12
                  [816,897,1077,1170,1134,1209,1245,1119,891,840,702,648,510]])      #13
# number of passengers arriving in 15 minitues

PsnVol0=(1./W)*PsnVol0# This is the arrival rate in per second

#----------------------Construction of OD matrix------------------------

#              1  2  3  4  5  6  7  8  9 10 11 12 13
OD0=np.array([[0, 5, 5, 5, 5, 5, 5, 5, 5,20,20,20,5],             #1
              [5, 0, 5, 5, 5, 5, 5, 5, 5,20,20,20,5],             #2
              [10,5, 0, 5, 5, 5, 5, 5, 5,20,20,20,5],             #3
              [15,5, 5, 0, 5, 5, 5, 5, 5,20,20,20,5],             #4
              [15,5, 5, 5, 0, 5, 5, 5, 5,20,20,20,5],             #5
              [15,5, 5, 5, 5, 0, 5, 5, 5,15,20,20,5],             #6
              [15,10, 5, 5, 5, 5, 0,5, 5,15,15,20,5],             #7
              [20,10, 5, 5, 5, 5, 5,0, 5,10,15,20,5],             #8
              [20,15,5, 5, 5, 5, 5, 5, 0, 5,10,20,5],             #9
              [20,10,10,5, 5, 5, 5, 5, 5, 0,5,15,10],             #10
              [20,15, 5,5, 5, 5, 5, 5, 5, 5, 0, 5,5],             #11
              [20,10, 5,5, 5, 5, 5, 5, 5, 5, 5, 0,5],             #12
              [20,10,10,5, 5, 5, 5, 5, 5, 5, 5, 5,0]])            #13
#heterogenous OD


#OD0=np.ones([N,N])
#for i in range(N):
#    OD0[i,i]=0
#homogenous OD

temp=OD0.sum(1)
OD1=np.zeros([N,N])
for i in range(N):
    OD1[i,]=np.divide(OD0[i,],float(temp[i]))

OD2=np.zeros([2*N,2*N])
for i in range(N):
    for j in range(i,N):
        OD2[i,j]=OD1[i,j]

for i in range(N,2*N):
    for j in range(i,2*N):
        OD2[i,j]=OD1[2*N-1-i,2*N-1-j]

#-------------------------------------OD2 is the OD matrix we will use------------------------------------


#---------------------------------------------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------------------------------------------
#--------------------------Simulation: generate the value of PsnInTrain and PsnAtStation at the end of disruption-----------
#---------------------------------------------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------------------------------------------

#------------------------------------the planned timetable before the disruption------------------------------------

OringinDepInt=210 # departure interval before the disruption

PlanedTimetableArr=np.zeros((40,2*stationNum))
PlanedTimetableDep=np.zeros((40,2*stationNum))

StdTimetableArr=np.zeros(2*stationNum)
StdTimetableDep=np.zeros(2*stationNum)

StdTimetableDep[0]=StdTimetableArr[0]+stationdt[0]
for i in range(1,N):
    StdTimetableArr[i]=StdTimetableDep[i-1]+linkrun10[i-1]
    StdTimetableDep[i]=StdTimetableArr[i]+stationdt[i]

StdTimetableArr[N]=StdTimetableDep[N-1]+backturning
StdTimetableDep[N]=StdTimetableArr[N]+stationdt[N]
for i in range(N+1,2*N):
    StdTimetableArr[i]=StdTimetableDep[i-1]+linkrun20[i-1-N]
    StdTimetableDep[i]=StdTimetableArr[i]+stationdt[i]

for k in range(40):
    for i in range(2*N):
        PlanedTimetableArr[k,i]=k*210+StdTimetableArr[i]
        PlanedTimetableDep[k,i]=k*210+StdTimetableDep[i]

#----------------------------passenger flow simulation before the disruption------------------------

nWantOn=np.zeros([2*N,2*N])# numbers of passengers arriving at station i with destination j during OringinDepInt(210s) in the normal situation before the disruption
nWantOn2=np.zeros([2*N,2*N])#passengers stranded in station i with destination j in the beginning of the recovery period
PsnAtStationStart=np.zeros(2*N) # passenger number stranded in stations

PsnInTrain=np.zeros([2*N,2*N])# the number of passengers from i to j in a service when the service is departing from a station in the normal situation before the disruption
PsnInTrain0=np.zeros([existingService,2*N]) # passengers in existing services in the beginning of the recovery period

for i in range(N):
    for j in range(i+1,N):
        nWantOn[i,j]=PsnVol0[i,H1]*OringinDepInt*OD2[i,j]

for i in range(N,2*N):
    for j in range(i+1,2*N):
        nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*OringinDepInt*OD2[i,j]

nCanOn=min(Capacity,sum(nWantOn[0]))#passenger number in service when the service departs from station 0
nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[0])# the ratio between "nCanOn" and the sum of "nWantOn"
PsnInTrain[0]=[round(nWantOn[0,j]*nCanOn_VS_nWantOn) for j in range(2*N)]# the number of passengers in a service when the service is departing from station 0

for i in range(1,N-1):
    nOff=PsnInTrain[i-1,i]
    nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i]))
    nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i])
    PsnInTrain[i]=[0*j for j in range(i+1)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+1,2*N)]

nCanOn=min(Capacity,sum(nWantOn[N]))#passenger number in train when the train departs from station N
nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[N])
PsnInTrain[N]=[round(nWantOn[N,j]*nCanOn_VS_nWantOn) for j in range(2*N)]

for i in range(N+1,2*N-1):
    nOff=PsnInTrain[i-1,i]
    nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i]))
    nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i])
    PsnInTrain[i]=[0*j for j in range(i+1)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+1,2*N)]

# In the normal situation, a passenger can always get on the first coming service.

#-----------------------------passenger flow simulation during disruption------------------------
startStatus=np.zeros((existingService))
startStatus[:existingService]=[ 25, 24, 23, 22, 21, 19, 18, 16, 15, 13,     12, 11,  10,  9,  8,  6,  5,  3,  2,  1]
# 'startStatus' is the locations (station or next station) of the existing services when the disruption is over

timeInstance=np.ones(existingService)*(-15)
# timeInstrance[k]<0, means that the train has arrived at corresponding station abs(timeInstrance[k]) ago.
# In the experiments, it is assumed that all the existing services are dwelling at stations during the disruption period.
# So, timeInstrance[k] can be any negative value.


#---------------------------downstream existing services---------------------------
temp=2*N-2 # station
for k in range(10):
    if startStatus[int(k)]==25:
        continue;
    for i in range(temp,int(startStatus[k]),-1):
        for j in range(i+2,2*N):#pay attention to i+2
            nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*Disruption*OD2[i,j]
            nWantOn2[i,j]=nWantOn[i,j]
        PsnAtStationStart[i]=sum(nWantOn2[i])
    temp=int(startStatus[k])-1
    i=int(startStatus[k])
    for j in range(i+2,2*N):
        nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*Disruption*OD2[i,j]
    if i<2*N-1:
        nOff=PsnInTrain[i-1,i]+PsnInTrain[i-1,i+1]
    else:
        nOff=PsnInTrain[i-1,i] #number of passengers getting off
    nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i]))
    nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i])
    PsnInTrain0[k]=[0*j for j in range(i+2)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+2,2*N)]
    nWantOn2[i]=[nWantOn[i,j]*(1-nCanOn_VS_nWantOn) for j in range(2*N)]
    PsnAtStationStart[i]=sum(nWantOn2[i])

if temp>=13:
    for i in range(temp,12,-1):
        for j in range(i+2,2*N):
            nWantOn[i,j]=PsnVol0[2*N-1-i,H1]*Disruption*OD2[i,j]
            nWantOn2[i,j]=nWantOn[i,j]
        PsnAtStationStart[i]=sum(nWantOn2[i])

#--------------------------------upstream existing services-----------------------
temp=N-2
for k in range(10,existingService):
    if startStatus[k]==13:
        continue;
    if startStatus[k]==12:
        continue;
    for i in range(temp,int(startStatus[k]),-1):
        for j in range(i+2,2*N):
            nWantOn[i,j]=PsnVol0[i,H1]*Disruption*OD2[i,j]
            nWantOn2[i,j]=nWantOn[i,j]
        PsnAtStationStart[i]=sum(nWantOn[i])
    temp=int(startStatus[k])-1
    i=int(startStatus[k])
    for j in range(i+2,N):
        nWantOn[i,j]=PsnVol0[i,H1]*Disruption*OD2[i,j]
    if i<N-1:
        nOff=PsnInTrain[i-1,i]+PsnInTrain[i-1,i+1]
    else:
        nOff=PsnInTrain[i-1,i]
    nCanOn=min(Capacity-sum(PsnInTrain[i-1])+nOff,sum(nWantOn[i]))
    nCanOn_VS_nWantOn=float(nCanOn)/sum(nWantOn[i])
    PsnInTrain0[k]=[0*j for j in range(i+2)]+[PsnInTrain[i-1,j]+round(nWantOn[i,j]*nCanOn_VS_nWantOn) for j in range(i+2,2*N)]
    nWantOn2[i]=[nWantOn[i,j]*(1-nCanOn_VS_nWantOn) for j in range(2*N)]
    PsnAtStationStart[i]=sum(nWantOn2[i])

if temp>=0:
    for i in range(temp,-1,-1):
        for j in range(i+2,N):
            nWantOn[i,j]=PsnVol0[i,H1]*Disruption*OD2[i,j]
            nWantOn2[i,j]=nWantOn[i,j]
        PsnAtStationStart[i]=sum(nWantOn2[i])