#### PROJECTED DEMAND FOR SCHOOL PLACES ####
#### PAN-LONDON MODEL                   ####

## July 2015, Wil Tonkiss
## Revised October 2015, Wil Tonkiss
## (c) Greater London Authority, 2015
## contact: demography@london.gov.uk

## Takes as inputs ward to ward pupil flows from the NPD and
## SYA ward-level population projections

## Calculates three scanarios for future demand
## Each scanario is run twice using 2 different assumptions about independent school attendence

rm(list=ls())
library(data.table)
library(dplyr)
library(tidyr)

        
#################################################################################################
        
#### DEFINE VARIABLES ####

WD <- "~/Models/Pan London SRP"                         

subDir <- substr(as.character(Sys.time()),1,16)
subDir <- gsub(":", "", subDir)
dir.create(file.path(WD, "Outputs", subDir), showWarnings = F)
WD.OUT <- file.path(WD, "Outputs", subDir)


SHLAA <- "ward_SHLAA_2014.csv"                          #ward level sya projection csv name (input projection)

secondary.npd.data  <- "Secondary_wardflows_2015.csv"   #processed NPD data for secondary pupils

primary.npd.data  <- "Primary_wardflows_2015.csv"       #processed NPD data for primary pupils


Roll.Year <- 2015         	                        #Roll is for Jan of year X. Taken in Aug X-1. Acdemic year (X-1)/X. Compare to MYE of X-1.
Output.Start.Year <- Roll.Year -1	                #Output projection from year
Output.End.Year <- Output.Start.Year + 10	        #Output projection to year



#######################################################################################


#### SOURCE FUNCTIONS ####

#This Main module calls a number of bespoke functions which are sourced here:

func.loc <- file.path(WD,"Model/R Scripts/Functions/")
func.list<- list.files(path=func.loc)

for(i in 1:length(func.list)){
        func <- paste(func.loc,func.list[i],sep="/")
        source(func)
}



### Read in data ####

setwd(file.path(WD,"Model/Inputs"))

#Read in 2014 rnd ward level Input.Projection capped projections
Input.Projection<-read.csv(SHLAA, header=T, nrow=1)
Input.Projection<-read.csv(SHLAA, header=T,
                colClasses=c(rep("character",4), rep("numeric",(ncol(Input.Projection)-4))))

Data.Col <- which(colnames(Input.Projection)=="X0")             #column with first sya data

# Years

Proj.First.Year <- min(Input.Projection$Year)                   #first year in input projection
Proj.Last.Year <- max(Input.Projection$Year)                    #last year in input projection
Proj.Rng <- Proj.Last.Year - Proj.First.Year + 1                #Determine the range of the input projection


#Read in ward to district lookup
ward.district.lookup <- read.csv("ward to district.csv", header=T, colClasses="character")

#A list of ward code for wards in the City of London
city.wards <- c(paste("E0500000",c(1:9),sep=""),paste("E050000",c(10:25),sep=""))



####################################################################################

####Function for Primary and Secondary####

Pan.London.Projections <- function(Cohort){

        if(Cohort == "Secondary"){
                LBound.Age <- 11
                UBound.Age <- 15
                Roll.File <- secondary.npd.data   #ward to ward pupil flows
        }
        
        if(Cohort == "Primary"){
                LBound.Age <- 4
                UBound.Age <- 10
                Roll.File <- primary.npd.data    #ward to ward pupil flows
        }
        
        

#### Read NPD data ####

#The inital input data files are read into memory

setwd(file.path(WD,"Model/Inputs"))

#Read in <ward of residence to ward of school> rolls for Jan 2014
#Remove any flows from Wales
NPD.Extract <- read.csv(Roll.File, header=T, colClasses=c("character","character",
                                                        "numeric","numeric",
                                                        "character","character")) %>%
        mutate(Country = substr(HomeWard,1,1))%>%
        filter(Country!="W")%>%
        select(-Country) %>%
        setnames(c("HomeWard","SchoolWard",
                "PupilFlow","PupilsOnRoll","SchoolLA","HomeLSOA"))


#3. A list of london wards
London.Wards <- filter(ward.district.lookup,substr(ward.district.lookup$District,1,3)=="E09")

#4. A list of London wards except those in City of London
London.Wards.Excl.City <- filter(London.Wards,District!="E09000001")

#Convert the above two ward selections to character vectors
London.Wards<-as.character(London.Wards$Ward)
London.Wards.Excl.City <-as.character(London.Wards.Excl.City$Ward)



#### CREATE DATAFRAMES ####

#Different functions require data in different conifigurations
#This section of code takes the input projection and roll data
#and creates the vectors and dataframes that are needed later

#1. A dataframe containing the total children on roll in each ward of residence

ward.roll.residents<- select(NPD.Extract,HomeWard,PupilsOnRoll) %>%
        unique()


#2. A cut of the input projection data containing only school-age children
#Ages included are detemined by whether the outputs are for primary
#or secondary school demand

#First combine the sya ages to make the desired cohort

Base.Popn <- select(Input.Projection,GSS.Code,Year) %>%
       mutate(Data = rowSums(select(Input.Projection,c((Data.Col + LBound.Age):(Data.Col+UBound.Age))))) %>%
       spread(Year, Data) %>%  #transpose data into 2D array with wards in rows and years in columns
       setnames(c("Ward",Proj.First.Year:Proj.Last.Year))


#### WARDS OUTSIDE LONDON ####

#Wards that aren't in London aren't included in the GLA projections
#Therefore projection data must be approximated using 2011 Census ward populations
#as a base and projecting forward using borough-level growth rates from the 2012 SNPP

#Extract a data table of wards not in London from the roll data
outside.london <- filter(NPD.Extract,!HomeWard %in% London.Wards) %>%
        select(HomeWard) %>%
        unique()


#If there are no wards from outside London in the rolls data there is no need to
#do anything further.

#Relevant dataframes and variables are passed to a function to
#calculate ward level projections for wards outside London
#which are then appended to the GLA ward projections 

if(nrow(outside.london)>0){
        
census.ward.pop <- read.csv("WARD POP CENSUS 2011.CSV",header=T,
                     colClasses=c("character","character","numeric")) %>%
        setnames(c("Ward","Age","Popn"))

snpp <- read.csv("2012 SNPP Population persons.csv",header=T,
                 colClasses=c(rep("character",5),
                              rep("numeric",27)))

outside.london <- Outside.London(ward.district.lookup, outside.london,
                                 LBound.Age, UBound.Age, Proj.Rng,
                                 census.ward.pop, snpp) %>%
        setnames(names(Base.Popn))
        

Base.Popn <- rbind(Base.Popn, outside.london)

}


#### REBASE PROJECTION FROM JUNE TO AUGUST ####

Base.Popn <- Rebase(Base.Popn,2)



#### REMOVE UNECESSARY YEARS FROM PROJECTION ####
        
#Now that data is availble all relevant wards the data
#for extra years can be removed

start.col <- Roll.Year - Proj.First.Year + 1   
col.rng <- Output.End.Year - Roll.Year                # Number of years to project

#Retain data for start year (year for which MYE and roll are comparable)
Indep.School.Estimates <- select(Base.Popn,1,start.col)                  
Base.Popn <- select(Base.Popn,Ward,start.col:(start.col+col.rng+1))      



#### CITY OF LONDON ####

#Convert ward to codes for wards in City to borough code
#for consistency with the projections input

change.ward.to.borough <- filter(NPD.Extract, HomeWard %in% city.wards) %>%
        select(-HomeWard) %>%
        mutate(HomeWard = "E09000001") %>%
        select(HomeWard,SchoolWard:HomeLSOA)

NPD.Extract <- filter(NPD.Extract, !HomeWard %in%  city.wards) %>%
        rbind(change.ward.to.borough)
        
rm(change.ward.to.borough)


#Exract ward totals for wards in City
city.total <- filter(ward.roll.residents, HomeWard %in% city.wards) %>%
                select(PupilsOnRoll) %>%
                sum()

#Sum all flows from wards in city to one borough-level flow
city <- data.frame(HomeWard="E09000001",PupilsOnRoll=city.total) %>%
        setnames(names(ward.roll.residents))

#Remove city wards from the master list of ward flows
ward.roll.residents <- filter(ward.roll.residents,!HomeWard %in% city.wards) %>%
        rbind(city)     #Add the borough level data to the master ward list

#Add the summed ward population back in
        
NPD.Extract$PupilsOnRoll <- ifelse(NPD.Extract$HomeWard=="E09000001", city.total, NPD.Extract$PupilsOnRoll)



#### CALCULATE NUMBER & PROPORTION AT INDEPENDENT SCHOOL ####

# constants are taken by comaparing Input.Projection (2013) with Roll (2014)

Indep.School.Estimates <- left_join(Indep.School.Estimates,ward.roll.residents,by=c("Ward"="HomeWard")) %>%
        setnames(c("Ward","Total.Pop","State.Pop")) %>%
        mutate(Indy.Pop = Total.Pop - State.Pop,
               Indy.pc = ((Total.Pop - State.Pop)/Total.Pop)*100,
               State.pc = 100-(((Total.Pop - State.Pop)/Total.Pop)*100))

Indep.School.Estimates$Ward <- as.character(Indep.School.Estimates$Ward)



#### DEAL WITH ROLL > ESTIMATE ####

#Where Roll population is greater than pop estimate add in
#the difference to each year of the estimate

#calculate top-up amount

Adjustment <- Indep.School.Estimates %>%
        mutate(adjustment = ifelse(Indep.School.Estimates$Indy.Pop<0,
                                   Indep.School.Estimates$Indy.Pop*-1,0)) %>%
        select(Ward, adjustment)


#Join top-ups to Input.Projection projection
nm <- names(Base.Popn)
Base.Popn <- left_join(Base.Popn,Adjustment,by="Ward")


#Add top-up constant to each year of projection
LastColumn <- ncol(Base.Popn)

Base.Popn <- cbind(select(Base.Popn,Ward),
                   lapply(Base.Popn[2:(LastColumn-1)], function(x) x+Base.Popn[LastColumn])) %>%     
                setnames(nm)

res.pop.for.output <- Base.Popn

rm(LastColumn,nm)


#Recreate the indy/state data based on the new projection numbers
#Shortcut is to set indy pop to 0 rather than doing calculation
       
Indep.School.Estimates <- mutate(Indep.School.Estimates,
                        temp = ifelse(Indep.School.Estimates$Indy.Pop<0, 0, Indep.School.Estimates$Indy.Pop),
                        temp1 = ifelse(Indep.School.Estimates$Indy.Pop<0, 0, Indep.School.Estimates$Indy.pc), 
                        temp2 = ifelse(Indep.School.Estimates$Indy.Pop<0, 100, Indep.School.Estimates$State.pc)) %>%
                        select(Ward,Total.Pop,State.Pop,temp,temp1,temp2) %>%
                        setnames(c("Ward","Total.Pop","State.Pop","Indy.Pop","Indy.pc","State.pc"))
  


##### DEMAND BY WARD OF RESIDENCE ####

# Determine school place demand if all future demand were met in ward of residence

col.labels <- c("Ward",Output.Start.Year:Output.End.Year)

Base.Popn <- left_join(Base.Popn,Indep.School.Estimates,by="Ward")

state.indep.demand.1 <- Demand.Residence(Base.Popn,type="abs")
state.indep.demand.2 <- Demand.Residence(Base.Popn,type="rate")

Residence.1 <- state.indep.demand.1[[1]]  %>%
        setnames(col.labels)
Independent.1 <-state.indep.demand.1[[2]]  %>%
        setnames(col.labels)

Residence.2 <- state.indep.demand.2[[1]] %>%
        setnames(col.labels)
Independent.2 <-state.indep.demand.2[[2]]  %>%
        setnames(col.labels)



#### DEMAND BY CURRENT MOBILITY PATTERNS ####

# Determine school place demand if all future demand were
# a function of current mobility distrution


#Determine home ward to school ward propensity rates

Mobility.Rates <- NPD.Extract %>%
        mutate(Prop=PupilFlow/PupilsOnRoll) %>%
        select(HomeWard,SchoolWard,Prop)


# Pass propensity rate and residence projection data to function

Mobility.1 <- Mobility.Demand(Mobility.Rates,Residence.1) %>%
        setnames(col.labels)
Mobility.2 <- Mobility.Demand(Mobility.Rates,Residence.2) %>%
        setnames(col.labels)



#### HYBRID ####

Hybrid.1 <- Hybrid(Residence.1,NPD.Extract) %>%
        setnames(col.labels)
Hybrid.2 <- Hybrid(Residence.2,NPD.Extract) %>%
        setnames(col.labels)



#Create a template for output data
#Ward, Ward Name, Borough, Borough Name
row.labels <- data.frame(Wd = c("E09000001", London.Wards.Excl.City)) %>%
        mutate(Ward = as.character(Wd))%>%
        select(Ward)%>%
        left_join(ward.district.lookup,by="Ward")





#### OUTPUTS ####


# Ward roll populations & borough roll projections & Independent School population estimates

#Column Labels for output files
col.labels <- c("Ward.Code","Ward.Name","Borough.Code","Borough.Name")

for(i in Output.Start.Year:Output.End.Year){
        col.labels <- c(col.labels,paste(i,"/",i-1999,sep=""))
}

output(Residence.1, "Residence", "N", col.labels, row.labels, Cohort, WD.OUT)
output(Residence.2, "Residence", "P", col.labels, row.labels, Cohort, WD.OUT)
output(Independent.1, "Independent", "N", col.labels, row.labels, Cohort, WD.OUT)
output(Independent.2, "Independent", "P", col.labels, row.labels, Cohort, WD.OUT)
output(Mobility.1, "Mobility", "N", col.labels, row.labels, Cohort, WD.OUT)
output(Mobility.2, "Mobility", "P", col.labels, row.labels, Cohort, WD.OUT)
output(Hybrid.1, "Hybrid", "N", col.labels, row.labels, Cohort, WD.OUT)
output(Hybrid.2, "Hybrid", "P", col.labels, row.labels, Cohort, WD.OUT)

output.residents.base(row.labels, res.pop.for.output, WD.OUT, Cohort)

}              #close Pan.London.Projections function



#### RUN MODEL FUNCTION ####

Pan.London.Projections("Primary")
Pan.London.Projections("Secondary")








