bs_combine.R

# script to load, analyze, and combine by scenario bs.csv files
# Peter Berrill
# May 2021
rm(list=ls()) # clear workspace i.e. remove saved variables
cat("\014") # clear console
library(dplyr)
library(ggplot2)
library(reshape2)
library(RColorBrewer)
setwd("~/Yale Courses/Research/Final Paper/HSM_github")
# read in large file containing 180k sample of the 2020 occupied housing stock, 
# generated by buildstock batch using the NREL ResStock housing characteristics data (https://github.com/NREL/resstock)
# bs2020<-read.csv('../resstock_projections/scen_bscsv/bs2020_180k.csv') # too large to store on github
# save(bs2020,file="Resstock_outputs/bs2020_180k.RData")
# files in the Resstock_outputs directory are too large to store on github, and this directory is included in .gitignore
# these large files can be requested from the author at peter.berrill@yale.edu
load("Resstock_outputs/bs2020_180k.RData") 
load("Resstock_outputs/bs_baseRFA.RData") # created by bs_adjust in the projection repository (https://github.com/NREL/resstock/tree/feature/projections/projection_scripts)
load("Intermediate_results/InitStock20.RData")
# similar to Resstock_outputs, HSM_results directory contains files too large to store on github and is therefore including in .gitignore
# Users can create the files in HSM_results themselves by running the hsm_cty script
load("HSM_results/US_smop_scenarios.RData") # us summary of stock model, created in the sm_result_vis script
load("Intermediate_results/ctycode.RData")

bsRFA<-bs_baseRFA[,1:113]

# check the distribution of housing types
ht<-table(bs2020$Geometry.Building.Type.RECS)/nrow(bs2020) # as expected, SF makes up about 68% of total occupied units
# lets show that it matches well
sum(ht[4:5]) # tot SF share from sample 68.2%
us_base[1,15]/us_base[1,14] # tot SF share from HSM, 68.2%

sum(ht[2:3]) # tot MF share from sample, 26.2%
us_base[1,16]/us_base[1,14] # tot MF share from HSM, 26.3%

ht[1] # MH share from sample, 5.5%
us_base[1,17]/us_base[1,14] # MH share from HSM, 5.5%

h20<-h20pc[-c(which(substr(h20pc$GeoID,1,2) %in% c("02","15"))),] # remove the counties in Alaska and Hawaii which are not part of the bs sample.
h20<-merge(h20,ctycode,by="GeoID")
# which counties are captured by the sample?
cc<-unique(bs2020$County)
# rsid<-unique(bs2020$rsid)
c<-which(!h20$RS_ID %in% cc)
# these ones are omitted, 24 low population counties with the 180k samples
h20$RS_ID[c]

ccfut<-unique(bs_baseRFA$County)
cfut<-which(!h20$RS_ID %in% ccfut)
# the following counties are omitted from future simulations, as there is insufficient new construction there for them to be included.
futomit<-h20$RS_ID[cfut]

cty_table<-table(bs2020$County)
# add description of simple type to bs file
bs2020$Type3<-"SF"
bs2020[bs2020$Geometry.Building.Type.RECS %in% c("Multi-Family with 2 - 4 Units","Multi-Family with 5+ Units"),]$Type3<-"MF"
bs2020[bs2020$Geometry.Building.Type.RECS == "Mobile Home",]$Type3<-"MH"
tcc_table<-as.data.frame(table(bs2020$County,bs2020$Vintage.ACS,bs2020$Type3))
names(tcc_table)<-c("County","Vintage ACS","Type3","count")
tcd_table<-as.data.frame(table(bs2020$Census.Division,bs2020$Vintage.ACS,bs2020$Type3))
names(tcd_table)<-c("Census.Division","Vintage ACS","Type3","count")

tcs_table<-as.data.frame(table(bs2020$State,bs2020$Vintage.ACS,bs2020$Type3))
names(tcs_table)<-c("State","Vintage ACS","Type3","count")

# do same for bsRFA
bsRFA$Type3<-"SF"
bsRFA[bsRFA$Geometry.Building.Type.RECS %in% c("Multi-Family with 2 - 4 Units","Multi-Family with 5+ Units"),]$Type3<-"MF"
bsRFA[bsRFA$Geometry.Building.Type.RECS == "Mobile Home",]$Type3<-"MH"
tcc_tableRFA<-as.data.frame(table(bsRFA$County,bsRFA$Type3))
names(tcc_tableRFA)<-c("County","Type3","count")
tcd_tableRFA<-as.data.frame(table(bsRFA$Census.Division,bsRFA$Type3))
names(tcd_tableRFA)<-c("Census.Division","Type3","count")
tcs_tableRFA<-as.data.frame(table(bsRFA$State,bsRFA$Type3))
names(tcs_tableRFA)<-c("State","Type3","count")

# add floor area estiamtes to each bs row, based on details in options_lookup (https://github.com/NREL/resstock) #######
bs2020$Floor.Area.m2<-0
bs2020[bs2020$Geometry.Floor.Area=="0-499"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(328/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="0-499"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(317/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="0-499"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(333/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="500-749"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(633/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="500-749"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(617/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="500-749"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(617/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="750-999"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(885/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="750-999"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(866/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="750-999"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(853/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="1000-1499"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(1220/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="1000-1499"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(1202/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="1000-1499"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(1138/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="1500-1999"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(1690/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="1500-1999"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(1675/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="1500-1999"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(1623/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="2000-2499"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(2176/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="2000-2499"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(2152/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="2000-2499"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(2115/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="2500-2999"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(2663/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="2500-2999"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(2631/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="2500-2999"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(2590/10.765,1)

bs2020[bs2020$Geometry.Floor.Area=="3000-3999"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(3301/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="3000-3999"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(3241/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="3000-3999"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(3138/10.765,1)
# 4000+. Using my own estimates here, consistent with my changes to options_lookup, but creating different value for MH
bs2020[bs2020$Geometry.Floor.Area=="4000+"&bs2020$Geometry.Building.Type.RECS %in% c("Single-Family Detached"),]$Floor.Area.m2<-round(7500/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="4000+"&bs2020$Geometry.Building.Type.RECS %in% c("Mobile Home"),]$Floor.Area.m2<-round(4200/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="4000+"&bs2020$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(7000/10.765,1)
bs2020[bs2020$Geometry.Floor.Area=="4000+"&bs2020$Type3=="MF",]$Floor.Area.m2<-round(7000/10.765,1)
# check distribution of house t-c in different regions ##########
h20$Occ_SF<-h20$Tot_HU_SF*rowSums(h20[,c(7:16)])
h20$Occ_MF<-h20$Tot_HU_MF*rowSums(h20[,c(27:36)])
h20$Occ_MH<-h20$Tot_HU_MH*rowSums(h20[,c(47:56)])
h20$Tot_Occ_Units<-h20$Occ_SF+h20$Occ_MF+h20$Occ_MH
locs<-unique(bs2020[,c("County","State","Census.Division")])
locs<-locs[order(locs$County),]
missing<-h20$RS_ID[c]
locsm<-data.frame("County"=missing,State=substr(missing,1,2))

StDiv<-unique(bs2020[,c("State","Census.Division")])
locsm<-merge(locsm,StDiv,by="State")
locs_all<-rbind(locs,locsm)
# add state and division to h20
h20<-merge(h20,locs,by.x = "RS_ID",by.y="County",all.x = TRUE) 

# RFA add floor area estiamtes to each bs row, based on details in options_lookup #######
bsRFA$Floor.Area.m2<-0
bsRFA[bsRFA$Geometry.Floor.Area=="0-499"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(328/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="0-499"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(317/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="0-499"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(333/10.765,1)

bsRFA[bsRFA$Geometry.Floor.Area=="500-749"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(633/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="500-749"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(617/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="500-749"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(617/10.765,1)

bsRFA[bsRFA$Geometry.Floor.Area=="750-999"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(885/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="750-999"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(866/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="750-999"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(853/10.765,1)

bsRFA[bsRFA$Geometry.Floor.Area=="1000-1499"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(1220/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="1000-1499"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(1202/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="1000-1499"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(1138/10.765,1)

bsRFA[bsRFA$Geometry.Floor.Area=="1500-1999"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(1690/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="1500-1999"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(1675/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="1500-1999"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(1623/10.765,1)

bsRFA[bsRFA$Geometry.Floor.Area=="2000-2499"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(2176/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="2000-2499"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(2152/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="2000-2499"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(2115/10.765,1)

bsRFA[bsRFA$Geometry.Floor.Area=="2500-2999"&bsRFA$Geometry.Building.Type.RECS %in% c("Single-Family Detached","Mobile Home"),]$Floor.Area.m2<-round(2663/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="2500-2999"&bsRFA$Geometry.Building.Type.RECS == "Single-Family Attached",]$Floor.Area.m2<-round(2631/10.765,1)
bsRFA[bsRFA$Geometry.Floor.Area=="2500-2999"&bsRFA$Type3=="MF",]$Floor.Area.m2<-round(2590/10.765,1)
# 
# calculate mean floor area by type, cohort, and geography #######
# type and cohort, nationally
tapply(bs2020$Floor.Area.m2, list(bs2020$Type3,bs2020$Vintage), mean)
meanFA_tc<-round(as.data.frame(tapply(bs2020$Floor.Area.m2, list(bs2020$Type3,bs2020$Vintage), mean)),2)
meanFA_tc$Type<-rownames(meanFA_tc)
mFA<-melt(meanFA_tc)
names(mFA)<-c("Type","Cohort","FloorArea")
mFA[mFA$Type=="MH" & mFA$Cohort=="<1940", ]$FloorArea<-80 # adjust floor area of <1940 MH
mFA[mFA$Type=="MH" & mFA$Cohort=="1940s", ]$FloorArea<-85 # adjust floor area of 1940s MH
windows()
ggplot(mFA,aes(x=Cohort,y=FloorArea,group=Type))+geom_line(aes(color=Type),size=1)+geom_point(aes(color=Type),size=2) + 
  labs(title = "Mean floor area of house types by cohort", y = "Mean Floor Area (m2)") + theme_bw() + ylim(c(10,260)) +
  theme(axis.text=element_text(size=11),axis.title=element_text(size=12,face = "bold"),plot.title = element_text(size = 12, face = "bold"))
# type and cohort, nationally, RFA
tapply(bsRFA$Floor.Area.m2, list(bsRFA$Type3,bsRFA$Vintage), mean)
meanFA_tc<-round(as.data.frame(tapply(bsRFA$Floor.Area.m2, list(bsRFA$Type3,bsRFA$Vintage), mean)),2)
meanFA_tc$Type<-rownames(meanFA_tc)
mFA_RFA<-melt(meanFA_tc)
names(mFA_RFA)<-c("Type","Cohort","FloorArea")
mFA$Scenario<-"1-4"
mFA_RFA$Scenario<-"5. Red. FA"

mFA[37:39,]<-mFA[34:36,]<-mFA[31:33,]<-mFA[28:30,]<-mFA[25:27,]
mFA$Cohort<-as.character(mFA$Cohort)
mFA[28:39,]$Cohort<-rep(c("2020s","2030s","2040s","2050s"),each=3)
mFA_all<-rbind(mFA,mFA_RFA)
mFA_all$Type_Scen<-paste(mFA_all$Type,mFA_all$Scenario,sep="_")
mFA_all$Characteristics<-"Base"
mFA_all[mFA_all$Scenario=="5. Red. FA",]$Characteristics<-"Reduced FA"
windows(width = 7.5,height = 5.3)
ggplot(mFA_all,aes(x=Cohort,y=FloorArea,group=Type_Scen))+geom_line(aes(color=Type,linetype=Characteristics),size=1)+geom_point(aes(color=Type),size=2) + scale_y_continuous(breaks = seq(0,250,50), limits = c(10,260),minor_breaks = NULL) +
  labs(title = "Mean floor area of house types by cohort", y = "Mean Floor Area (m2)") + theme_bw() + #ylim(c(10,260)) +
  theme(axis.text=element_text(size=10),axis.title=element_text(size=11,face = "bold"),plot.title = element_text(size = 13, face = "bold"), 
        legend.title=element_text(size=11), legend.text=element_text(size=10))
mFA_all[mFA_all$Scenario=="5. Red. FA",]$Scenario<-"5-6"
windows()
ggplot(mFA_all,aes(x=Cohort,y=FloorArea,group=Type_Scen))+geom_line(aes(color=Type,linetype=Scenario),size=1)+geom_point(aes(color=Type),size=2) + 
  labs(title = "Mean floor area of house types by cohort", y = "Mean Floor Area (m2)") + theme_bw() + ylim(c(10,260)) +
  theme(axis.text=element_text(size=11),axis.title=element_text(size=12,face = "bold"),plot.title = element_text(size = 13, face = "bold"), 
        legend.title=element_text(size=11), legend.text=element_text(size=10))
write.csv(mFA_all[,c('Cohort','Type','Characteristics','FloorArea')],file = '~/projects/Yale/resstock_projections/Figure_Results_Data/ED_Fig1.csv',row.names = FALSE)
# mFA by geograpy, base #####
# mean floor area by type, cohort, and region
meanFA_tcr<-round(as.data.frame(tapply(bs2020$Floor.Area.m2, list(bs2020$Census.Region,bs2020$Type3,bs2020$Vintage.ACS), mean)),2)
meanFA_tcr$`MH.<1940`<-80 # adjust floor area of <1940 MH
# at the regional level, if there is no data for a cohort, use the average for all non NA regions for that cohort for all regions. This just affects MH
for (k in 1:ncol(meanFA_tcr)) {if(any(is.na(meanFA_tcr[,k]))) {meanFA_tcr[,k]<- round(mean(meanFA_tcr[,k],na.rm = TRUE),2) }}

meanFA_tcr$Region<-rownames(meanFA_tcr)
rownames(meanFA_tcr)<-1:nrow(meanFA_tcr)
meanFA_tcr<-meanFA_tcr[,c(19,1:18)]
# type, cohort, and division
meanFA_tcd<-round(as.data.frame(tapply(bs2020$Floor.Area.m2, list(bs2020$Census.Division,bs2020$Type3,bs2020$Vintage.ACS), mean)),2)
meanFA_tcd$`MH.<1940`<-80
meanFA_tcd$Census.Division<-rownames(meanFA_tcd)
rownames(meanFA_tcd)<-1:nrow(meanFA_tcd)
DivReg<-unique(bs2020[,c("Census.Division","Census.Region")])
meanFA_tcd<-merge(meanFA_tcd,DivReg)
tcds<-tcd_table[tcd_table$count<5,] # which tcd combos have a small # of observations?
for (j in 1:nrow(tcds)) {
  meanFA_tcd[meanFA_tcd$Census.Division==tcds$Census.Division[j], paste(tcds$Type3[j],tcds$`Vintage ACS`[j],sep = ".")]<-NA
}

# fill in blanks based on regional average
for (k in 2:19) { 
  for (r in 1:9) {
   if (is.na(meanFA_tcd[r,k])) {
     meanFA_tcd[r,k]<-round(meanFA_tcr[meanFA_tcr$Region==meanFA_tcd[r,"Census.Region"],k],2)
   } 
  }
}
# type, cohort, and state
meanFA_tcs<-round(as.data.frame(tapply(bs2020$Floor.Area.m2, list(bs2020$State,bs2020$Type3,bs2020$Vintage.ACS), mean)),2)
meanFA_tcs$State<-rownames(meanFA_tcs)
rownames(meanFA_tcs)<-1:nrow(meanFA_tcs)
meanFA_tcs<-meanFA_tcs[,c(19,1:18)]
meanFA_tcs<-merge(meanFA_tcs,StDiv)
tcss<-tcs_table[tcs_table$count<5,] # which tcs combos have a small # of observations?
for (j in 1:nrow(tcss)) {
  meanFA_tcs[meanFA_tcs$State==tcss$State[j], paste(tcss$Type3[j],tcss$`Vintage ACS`[j],sep = ".")]<-NA
}
# fill in blanks based on division average
for (k in 2:19) { 
  for (r in 1:49) {
    if (is.na(meanFA_tcs[r,k])) {
      meanFA_tcs[r,k]<-round(meanFA_tcd[meanFA_tcd$Census.Division==meanFA_tcs[r,"Census.Division"],k],2)
    } 
  }
}

# type, cohort, and county
meanFA_tcc<-round(as.data.frame(tapply(bs2020$Floor.Area.m2, list(bs2020$County,bs2020$Type3,bs2020$Vintage.ACS), mean)),2)
meanFA_tcc$RS_ID<-rownames(meanFA_tcc)
rownames(meanFA_tcc)<-1:nrow(meanFA_tcc)
meanFA_tcc<-merge(meanFA_tcc,ctycode,all = TRUE)
meanFA_tcc<-meanFA_tcc[,c(20,1:19)]
meanFA_tcc<-meanFA_tcc[order(meanFA_tcc$GeoID),]
meanFA_tcc<-meanFA_tcc[-c(which(substr(meanFA_tcc$GeoID,1,2) %in% c("02","15"))),]
tccs<-tcc_table[tcc_table$count<5,] # which tcc combos have a small # of observations?
for (j in 1:nrow(tccs)) {
  meanFA_tcc[meanFA_tcc$RS_ID==tccs$County[j], paste(tccs$Type3[j],tccs$`Vintage ACS`[j],sep = ".")]<-NA
}

meanFA_tcc$State<-substr(meanFA_tcc$RS_ID,1,2)
# fill in blanks based on state average
for (k in 3:20) { 
  for (r in 1:3108) {
    if (is.na(meanFA_tcc[r,k])) {
      meanFA_tcc[r,k]<-round(meanFA_tcs[meanFA_tcs$State==meanFA_tcc[r,"State"],k-1],2)
    } 
  }
}

# mFA by geograpy, RFA #####

# mean floor area by type, cohort, and region
# meanFA_tcr_RFA<-round(as.data.frame(tapply(bsRFA$Floor.Area.m2, list(bsRFA$Census.Region,bsRFA$Type3,bsRFA$Vintage.ACS), mean)),2)
meanFA_tcr_RFA<-round(as.data.frame(tapply(bsRFA$Floor.Area.m2, list(bsRFA$Census.Region,bsRFA$Type3), mean)),2)

meanFA_tcr_RFA$Region<-rownames(meanFA_tcr_RFA)
rownames(meanFA_tcr_RFA)<-1:nrow(meanFA_tcr_RFA)
meanFA_tcr_RFA<-meanFA_tr_RFA[,c(4,1:3)]
# type, cohort, and division
meanFA_tcd_RFA<-round(as.data.frame(tapply(bsRFA$Floor.Area.m2, list(bsRFA$Census.Division,bsRFA$Type3), mean)),2)
meanFA_tcd_RFA$Census.Division<-rownames(meanFA_tcd_RFA)
rownames(meanFA_tcd_RFA)<-1:nrow(meanFA_tcd_RFA)
DivReg<-unique(bsRFA[,c("Census.Division","Census.Region")])
meanFA_tcd_RFA<-merge(meanFA_tcd_RFA,DivReg)
tcd_RFAs<-tcd_tableRFA[tcd_tableRFA$count<5,] # which td_RFA combos have a small # of observations?
# only needs done if td_RFAs was > 0
# for (j in 1:nrow(tcd_RFAs)) {
#   meanFA_tcd_RFA[meanFA_tcd_RFA$Census.Division==tcd_RFAs$Census.Division[j], paste(tcd_RFAs$Type3[j],tcd_RFAs$`Vintage ACS`[j],sep = ".")]<-NA
# }

# # fill in blanks based on regional average
# for (k in 2:19) { 
#   for (r in 1:9) {
#     if (is.na(meanFA_tcd_RFA[r,k])) {
#       meanFA_tcd_RFA[r,k]<-round(meanFA_tcr_RFA[meanFA_tcr_RFA$Region==meanFA_tcd_RFA[r,"Census.Region"],k],2)
#     } 
#   }
# }
# type, cohort, and state
meanFA_tcs_RFA<-round(as.data.frame(tapply(bsRFA$Floor.Area.m2, list(bsRFA$State,bsRFA$Type3), mean)),2)
meanFA_tcs_RFA$State<-rownames(meanFA_tcs_RFA)
rownames(meanFA_tcs_RFA)<-1:nrow(meanFA_tcs_RFA)
meanFA_tcs_RFA<-meanFA_tcs_RFA[,c(4,1:3)]
meanFA_tcs_RFA<-merge(meanFA_tcs_RFA,StDiv)
tcs_RFAs<-tcs_tableRFA[tcs_tableRFA$count<5,] # which tcs_RFA combos have a small # of observations?
for (j in 1:nrow(tcs_RFAs)) {
  meanFA_tcs_RFA[meanFA_tcs_RFA$State==tcs_RFAs$State[j], paste(tcs_RFAs$Type3[j],sep = ".")]<-NA
}
# fill in blanks based on division average
for (k in 2:4) { 
  for (r in 1:49) {
    if (is.na(meanFA_tcs_RFA[r,k])) {
      meanFA_tcs_RFA[r,k]<-round(meanFA_tcd_RFA[meanFA_tcd_RFA$Census.Division==meanFA_tcs_RFA[r,"Census.Division"],k],2)
    } 
  }
}

# type, cohort, and county
meanFA_tcc_RFA<-round(as.data.frame(tapply(bsRFA$Floor.Area.m2, list(bsRFA$County,bsRFA$Type3), mean)),2)
meanFA_tcc_RFA$RS_ID<-rownames(meanFA_tcc_RFA)
rownames(meanFA_tcc_RFA)<-1:nrow(meanFA_tcc_RFA)
meanFA_tcc_RFA<-merge(meanFA_tcc_RFA,ctycode,all = TRUE)
meanFA_tcc_RFA<-meanFA_tcc_RFA[,c(5,1:4)]
meanFA_tcc_RFA<-meanFA_tcc_RFA[order(meanFA_tcc_RFA$GeoID),]
meanFA_tcc_RFA<-meanFA_tcc_RFA[-c(which(substr(meanFA_tcc_RFA$GeoID,1,2) %in% c("02","15"))),] # remove AK, HI
tcc_RFAs<-tcc_tableRFA[tcc_tableRFA$count<5,] # which tcc_RFA combos have a small # of observations?
for (j in 1:nrow(tcc_RFAs)) {
  meanFA_tcc_RFA[meanFA_tcc_RFA$RS_ID==tcc_RFAs$County[j], paste(tcc_RFAs$Type3[j],sep = ".")]<-NA
}

meanFA_tcc_RFA$State<-substr(meanFA_tcc_RFA$RS_ID,1,2)
# fill in blanks based on state average
for (k in 3:5) { 
  for (r in 1:3108) {
    if (is.na(meanFA_tcc_RFA[r,k])) {
      meanFA_tcc_RFA[r,k]<-round(meanFA_tcs_RFA[meanFA_tcs_RFA$State==meanFA_tcc_RFA[r,"State"],k-1],2)
    } 
  }
}

# alter meanFA dfs so that the RFA version never has larger floor areas than the regular version. THIS IS CRUCIAL
# needs to be done for every future cohort
for (j in 1:3108) { 
  if (meanFA_tcc_RFA$MF[j]>meanFA_tcc$MF.2010s[j]) {
    meanFA_tcc_RFA$MF[j]<-meanFA_tcc$MF.2010s[j]
  }
  if (meanFA_tcc_RFA$SF[j]>meanFA_tcc$SF.2010s[j]) {
    meanFA_tcc_RFA$SF[j]<-meanFA_tcc$SF.2010s[j]
  }
  if (meanFA_tcc_RFA$MH[j]>meanFA_tcc$MH.2010s[j]) {
    meanFA_tcc_RFA$MH[j]<-meanFA_tcc$MH.2010s[j]
  }
}

# area and material  calculations ############
load("HSM_results/County_Scenario_SM_Results.RData") # 
SF_dem_factor<-0.35 # how much of SF homes leaving the stock are actually demolished right away?
MF_dem_factor<-0.2 # how much of MF homes leaving the stock are actually demolished right away?
MH_dem_factor<-0.45 # how much of MH homes leaving the stock are actually demolished right away? This was changed from 0.5 to 0.45
# add the resstock county ID's
smb<-smop_base
smop_base<-merge(smop_base,ctycode)
smop_hiDR<-merge(smop_hiDR,ctycode)
smop_hiMF<-merge(smop_hiMF,ctycode)
smop_hiDRMF<-merge(smop_hiDRMF,ctycode)
# remove AK and HI
smop_base<-smop_base[-c(which(substr(smop_base$GeoID,1,2) %in% c("02","15"))),]
smop_hiDR<-smop_hiDR[-c(which(substr(smop_hiDR$GeoID,1,2) %in% c("02","15"))),]
smop_hiMF<-smop_hiMF[-c(which(substr(smop_hiMF$GeoID,1,2) %in% c("02","15"))),]
smop_hiDRMF<-smop_hiDRMF[-c(which(substr(smop_hiDRMF$GeoID,1,2) %in% c("02","15"))),]

# define the RFA data frames, only for the first 3 stock scenarios
smop_RFA<-smop_base
smop_hiDR_RFA<-smop_hiDR
smop_hiMF_RFA<-smop_hiMF

all.equal(meanFA_tcc$GeoID,smop_base$GeoID) # if this is true, can just transplant the mean FA into the smop files
all.equal(meanFA_tcc$GeoID,smop_hiDR$GeoID) # if this is true, can just transplant the mean FA into the smop files
all.equal(meanFA_tcc$GeoID,smop_hiMF$GeoID) # if this is true, can just transplant the mean FA into the smop files
all.equal(meanFA_tcc$GeoID,smop_hiDRMF$GeoID) # if this is true, can just transplant the mean FA into the smop files
all.equal(meanFA_tcc_RFA$GeoID,smop_RFA$GeoID) # if this is true, can just transplant the mean FA into the smop files
q<-which(is.na(meanFA_tcc[,3:20]))
# if no na.s in q, add average floor areas per type and ACS vintage
smop_base[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]
smop_hiDR[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]
smop_hiMF[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]
smop_hiDRMF[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]
smop_RFA[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]
smop_hiDR_RFA[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]
smop_hiMF_RFA[,paste("FA",names(meanFA_tcc)[3:20],sep=".")]<-meanFA_tcc[,3:20]

# define different floor areas for new construction in the RFA scenarios
names(meanFA_tcc_RFA)[3:5]<-paste(names(meanFA_tcc_RFA)[3:5],'NC',sep='.')
smop_RFA[,paste("FA",names(meanFA_tcc_RFA)[3:5],sep=".")]<-meanFA_tcc_RFA[,3:5]
smop_hiDR_RFA[,paste("FA",names(meanFA_tcc_RFA)[3:5],sep=".")]<-meanFA_tcc_RFA[,3:5]
smop_hiMF_RFA[,paste("FA",names(meanFA_tcc_RFA)[3:5],sep=".")]<-meanFA_tcc_RFA[,3:5]

# now calculate m2/cap per house type in each county for each scenario ###########
# base
for (i in c(1:3108)) { # these loops are slow, take about 4 seconds per county
  print(i)
  smop_base[[3]][[i]]<-as.data.frame(smop_base[[3]][[i]])
  
  smop_base[[3]][[i]]$m2cap<-smop_base[[3]][[i]]$m2cap_MH<-smop_base[[3]][[i]]$m2cap_MF<- smop_base[[3]][[i]]$m2cap_SF<-
    smop_base[[3]][[i]]$Occ_m2<-smop_base[[3]][[i]]$Occ_m2_MH<-smop_base[[3]][[i]]$Occ_m2_MF<-smop_base[[3]][[i]]$Occ_m2_SF<-
    smop_base[[3]][[i]]$Tot_NC_m2<-smop_base[[3]][[i]]$NC_MH_m2<-smop_base[[3]][[i]]$NC_MF_m2<-smop_base[[3]][[i]]$NC_SF_m2<-
    smop_base[[3]][[i]]$Tot_Dem_m2<-smop_base[[3]][[i]]$Dem_MH_m2<-smop_base[[3]][[i]]$Dem_MF_m2<-smop_base[[3]][[i]]$Dem_SF_m2<-
    smop_base[[3]][[i]]$Tot_NC<-smop_base[[3]][[i]]$NC_MH<-smop_base[[3]][[i]]$NC_MF<-smop_base[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_base[[3]][[i]]$NC_SF[y]<-smop_base[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_base[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_base[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_base[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_base[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_base[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_base[[3]][[i]]$NC_MF[y]<-smop_base[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_base[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_base[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_base[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_base[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_base[[3]][[i]]$NC_MH[y]<-smop_base[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_base[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_base[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_base[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_base[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_base[[3]][[i]]$NC_SF[y]<-smop_base[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_base[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_base[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_base[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_base[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_base[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_base[[3]][[i]]$NC_MF[y]<-smop_base[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_base[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_base[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_base[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_base[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_base[[3]][[i]]$NC_MH[y]<-smop_base[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_base[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_base[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_base[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_base[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_base[[3]][[i]]$NC_SF[y]<-smop_base[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_base[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_base[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_base[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_base[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_base[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_base[[3]][[i]]$NC_MF[y]<-smop_base[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_base[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_base[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_base[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_base[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_base[[3]][[i]]$NC_MH[y]<-smop_base[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_base[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_base[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_base[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_base[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_base[[3]][[i]]$NC_SF[y]<-smop_base[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_base[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_base[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_base[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_base[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_base[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_base[[3]][[i]]$NC_MF[y]<-smop_base[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_base[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_base[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_base[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_base[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_base[[3]][[i]]$NC_MH[y]<-smop_base[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_base[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_base[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_base[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_base[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_base[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_base[[3]][[i]][,316:318][which(smop_base[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0 # make sure no negative inflows
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_base[[3]][[i]]$NC_SF[1:40]/smop_base[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_base[[3]][[i]]$NC_SF[1:40]/smop_base[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_base[[3]][[i]]$NC_SF[1:40]/smop_base[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_base[[3]][[i]]$NC_SF[1:40]/smop_base[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_base[[3]][[i]]$NC_SF[1:40]<-smop_base[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_base[[3]][[i]]$NC_SF[1:40])/sum(smop_base[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_base[[3]][[i]]$NC_MF[1:40]/smop_base[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_base[[3]][[i]]$NC_MF[1:40]/smop_base[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_base[[3]][[i]]$NC_MF[1:40]/smop_base[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_base[[3]][[i]]$NC_MF[1:40]/smop_base[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_base[[3]][[i]]$NC_MF[1:40]<-smop_base[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_base[[3]][[i]]$NC_MF[1:40])/sum(smop_base[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_base[[3]][[i]]$NC_MH[1:40]/smop_base[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_base[[3]][[i]]$NC_MH[1:40]/smop_base[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_base[[3]][[i]]$NC_MH[1:40]/smop_base[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_base[[3]][[i]]$NC_MH[1:40]/smop_base[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_base[[3]][[i]]$NC_MH[1:40]<-smop_base[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_base[[3]][[i]]$NC_MH[1:40])/sum(smop_base[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_base[[3]][[i]]$Tot_NC<-smop_base[[3]][[i]]$NC_SF+smop_base[[3]][[i]]$NC_MF+smop_base[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_base[[3]][[i]]$NC_SF_m2<-smop_base[[3]][[i]]$NC_SF*smop_base$FA.SF.2010s[i]
  smop_base[[3]][[i]]$NC_MF_m2<-smop_base[[3]][[i]]$NC_MF*smop_base$FA.MF.2010s[i]
  smop_base[[3]][[i]]$NC_MH_m2<-smop_base[[3]][[i]]$NC_MH*smop_base$FA.MH.2010s[i]
  smop_base[[3]][[i]]$Tot_NC_m2<-smop_base[[3]][[i]]$NC_SF_m2+smop_base[[3]][[i]]$NC_MF_m2+smop_base[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_base[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_base[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_base[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))+
    rowSums(smop_base[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_base[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_base[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_base[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_base[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))+
    rowSums(smop_base[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_base[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_base[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_base[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_base[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))+
    rowSums(smop_base[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_base[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_base[[3]][[i]]$Tot_Dem_m2<-smop_base[[3]][[i]]$Dem_SF_m2+smop_base[[3]][[i]]$Dem_MF_m2+smop_base[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type, assuming that 2010 m2/house are constant through 2050s
  smop_base[[3]][[i]]$Occ_m2_SF<-rowSums(smop_base[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_base[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_base[[3]][[i]]$Occ_m2_MF<-rowSums(smop_base[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_base[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_base[[3]][[i]]$Occ_m2_MH<-rowSums(smop_base[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_base[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_base[[3]][[i]]$Occ_m2<-smop_base[[3]][[i]]$Occ_m2_SF+smop_base[[3]][[i]]$Occ_m2_MF+smop_base[[3]][[i]]$Occ_m2_MH
  
  smop_base[[3]][[i]]$m2cap_SF<-smop_base[[3]][[i]]$Occ_m2_SF/smop_base[[3]][[i]]$Pop_SF
  smop_base[[3]][[i]]$m2cap_MF<-smop_base[[3]][[i]]$Occ_m2_MF/smop_base[[3]][[i]]$Pop_MF
  smop_base[[3]][[i]]$m2cap_MH<-smop_base[[3]][[i]]$Occ_m2_MH/smop_base[[3]][[i]]$Pop_MH
  smop_base[[3]][[i]]$m2cap<-smop_base[[3]][[i]]$Occ_m2/smop_base[[3]][[i]]$Population
}
# hiDR
for (i in c(1:3108)) {
  print(i)
  smop_hiDR[[3]][[i]]<-as.data.frame(smop_hiDR[[3]][[i]])
  
  smop_hiDR[[3]][[i]]$m2cap<-smop_hiDR[[3]][[i]]$m2cap_MH<-smop_hiDR[[3]][[i]]$m2cap_MF<- smop_hiDR[[3]][[i]]$m2cap_SF<-
    smop_hiDR[[3]][[i]]$Occ_m2<-smop_hiDR[[3]][[i]]$Occ_m2_MH<-smop_hiDR[[3]][[i]]$Occ_m2_MF<-smop_hiDR[[3]][[i]]$Occ_m2_SF<-
    smop_hiDR[[3]][[i]]$Tot_NC_m2<-smop_hiDR[[3]][[i]]$NC_MH_m2<-smop_hiDR[[3]][[i]]$NC_MF_m2<-smop_hiDR[[3]][[i]]$NC_SF_m2<-
    smop_hiDR[[3]][[i]]$Tot_Dem_m2<-smop_hiDR[[3]][[i]]$Dem_MH_m2<-smop_hiDR[[3]][[i]]$Dem_MF_m2<-smop_hiDR[[3]][[i]]$Dem_SF_m2<-
    smop_hiDR[[3]][[i]]$Tot_NC<-smop_hiDR[[3]][[i]]$NC_MH<-smop_hiDR[[3]][[i]]$NC_MF<-smop_hiDR[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_hiDR[[3]][[i]]$NC_SF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_hiDR[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_hiDR[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_hiDR[[3]][[i]]$NC_MF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_hiDR[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_hiDR[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_hiDR[[3]][[i]]$NC_MH[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_hiDR[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_hiDR[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_hiDR[[3]][[i]]$NC_SF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_hiDR[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_hiDR[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_hiDR[[3]][[i]]$NC_MF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_hiDR[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_hiDR[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_hiDR[[3]][[i]]$NC_MH[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_hiDR[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_hiDR[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_hiDR[[3]][[i]]$NC_SF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_hiDR[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_hiDR[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_hiDR[[3]][[i]]$NC_MF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_hiDR[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_hiDR[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_hiDR[[3]][[i]]$NC_MH[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_hiDR[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_hiDR[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_hiDR[[3]][[i]]$NC_SF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_hiDR[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_hiDR[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_hiDR[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_hiDR[[3]][[i]]$NC_MF[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_hiDR[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_hiDR[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_hiDR[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_hiDR[[3]][[i]]$NC_MH[y]<-smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_hiDR[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_hiDR[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_hiDR[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_hiDR[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_hiDR[[3]][[i]][,316:318][which(smop_hiDR[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_hiDR[[3]][[i]]$NC_SF[1:40]/smop_hiDR[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_hiDR[[3]][[i]]$NC_SF[1:40]/smop_hiDR[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_hiDR[[3]][[i]]$NC_SF[1:40]/smop_hiDR[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_hiDR[[3]][[i]]$NC_SF[1:40]/smop_hiDR[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_hiDR[[3]][[i]]$NC_SF[1:40]<-smop_hiDR[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_hiDR[[3]][[i]]$NC_SF[1:40])/sum(smop_hiDR[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_hiDR[[3]][[i]]$NC_MF[1:40]/smop_hiDR[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_hiDR[[3]][[i]]$NC_MF[1:40]/smop_hiDR[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_hiDR[[3]][[i]]$NC_MF[1:40]/smop_hiDR[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_hiDR[[3]][[i]]$NC_MF[1:40]/smop_hiDR[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_hiDR[[3]][[i]]$NC_MF[1:40]<-smop_hiDR[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_hiDR[[3]][[i]]$NC_MF[1:40])/sum(smop_hiDR[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_hiDR[[3]][[i]]$NC_MH[1:40]/smop_hiDR[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_hiDR[[3]][[i]]$NC_MH[1:40]/smop_hiDR[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_hiDR[[3]][[i]]$NC_MH[1:40]/smop_hiDR[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_hiDR[[3]][[i]]$NC_MH[1:40]/smop_hiDR[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_hiDR[[3]][[i]]$NC_MH[1:40]<-smop_hiDR[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_hiDR[[3]][[i]]$NC_MH[1:40])/sum(smop_hiDR[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_hiDR[[3]][[i]]$Tot_NC<-smop_hiDR[[3]][[i]]$NC_SF+smop_hiDR[[3]][[i]]$NC_MF+smop_hiDR[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_hiDR[[3]][[i]]$NC_SF_m2<-smop_hiDR[[3]][[i]]$NC_SF*smop_hiDR$FA.SF.2010s[i]
  smop_hiDR[[3]][[i]]$NC_MF_m2<-smop_hiDR[[3]][[i]]$NC_MF*smop_hiDR$FA.MF.2010s[i]
  smop_hiDR[[3]][[i]]$NC_MH_m2<-smop_hiDR[[3]][[i]]$NC_MH*smop_hiDR$FA.MH.2010s[i]
  smop_hiDR[[3]][[i]]$Tot_NC_m2<-smop_hiDR[[3]][[i]]$NC_SF_m2+smop_hiDR[[3]][[i]]$NC_MF_m2+smop_hiDR[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_hiDR[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_hiDR[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_hiDR[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))+
    rowSums(smop_hiDR[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_hiDR[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_hiDR[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_hiDR[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_hiDR[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))+
    rowSums(smop_hiDR[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_hiDR[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_hiDR[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_hiDR[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_hiDR[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))+
    rowSums(smop_hiDR[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_hiDR[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_hiDR[[3]][[i]]$Tot_Dem_m2<-smop_hiDR[[3]][[i]]$Dem_SF_m2+smop_hiDR[[3]][[i]]$Dem_MF_m2+smop_hiDR[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type, assuming that 2010 m2/house are constant through 2050s
  smop_hiDR[[3]][[i]]$Occ_m2_SF<-rowSums(smop_hiDR[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_hiDR[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_hiDR[[3]][[i]]$Occ_m2_MF<-rowSums(smop_hiDR[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_hiDR[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_hiDR[[3]][[i]]$Occ_m2_MH<-rowSums(smop_hiDR[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_hiDR[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_hiDR[[3]][[i]]$Occ_m2<-smop_hiDR[[3]][[i]]$Occ_m2_SF+smop_hiDR[[3]][[i]]$Occ_m2_MF+smop_hiDR[[3]][[i]]$Occ_m2_MH
  
  smop_hiDR[[3]][[i]]$m2cap_SF<-smop_hiDR[[3]][[i]]$Occ_m2_SF/smop_hiDR[[3]][[i]]$Pop_SF
  smop_hiDR[[3]][[i]]$m2cap_MF<-smop_hiDR[[3]][[i]]$Occ_m2_MF/smop_hiDR[[3]][[i]]$Pop_MF
  smop_hiDR[[3]][[i]]$m2cap_MH<-smop_hiDR[[3]][[i]]$Occ_m2_MH/smop_hiDR[[3]][[i]]$Pop_MH
  smop_hiDR[[3]][[i]]$m2cap<-smop_hiDR[[3]][[i]]$Occ_m2/smop_hiDR[[3]][[i]]$Population
}
# hiMF
for (i in c(1:3108)) {
  print(i)
  smop_hiMF[[3]][[i]]<-as.data.frame(smop_hiMF[[3]][[i]])
  
  smop_hiMF[[3]][[i]]$m2cap<-smop_hiMF[[3]][[i]]$m2cap_MH<-smop_hiMF[[3]][[i]]$m2cap_MF<- smop_hiMF[[3]][[i]]$m2cap_SF<-
    smop_hiMF[[3]][[i]]$Occ_m2<-smop_hiMF[[3]][[i]]$Occ_m2_MH<-smop_hiMF[[3]][[i]]$Occ_m2_MF<-smop_hiMF[[3]][[i]]$Occ_m2_SF<-
    smop_hiMF[[3]][[i]]$Tot_NC_m2<-smop_hiMF[[3]][[i]]$NC_MH_m2<-smop_hiMF[[3]][[i]]$NC_MF_m2<-smop_hiMF[[3]][[i]]$NC_SF_m2<-
    smop_hiMF[[3]][[i]]$Tot_Dem_m2<-smop_hiMF[[3]][[i]]$Dem_MH_m2<-smop_hiMF[[3]][[i]]$Dem_MF_m2<-smop_hiMF[[3]][[i]]$Dem_SF_m2<-
    smop_hiMF[[3]][[i]]$Tot_NC<-smop_hiMF[[3]][[i]]$NC_MH<-smop_hiMF[[3]][[i]]$NC_MF<-smop_hiMF[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_hiMF[[3]][[i]]$NC_SF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_hiMF[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_hiMF[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_hiMF[[3]][[i]]$NC_MF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_hiMF[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_hiMF[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_hiMF[[3]][[i]]$NC_MH[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_hiMF[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_hiMF[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_hiMF[[3]][[i]]$NC_SF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_hiMF[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_hiMF[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_hiMF[[3]][[i]]$NC_MF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_hiMF[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_hiMF[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_hiMF[[3]][[i]]$NC_MH[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_hiMF[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_hiMF[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_hiMF[[3]][[i]]$NC_SF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_hiMF[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_hiMF[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_hiMF[[3]][[i]]$NC_MF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_hiMF[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_hiMF[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_hiMF[[3]][[i]]$NC_MH[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_hiMF[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_hiMF[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_hiMF[[3]][[i]]$NC_SF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_hiMF[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_hiMF[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_hiMF[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_hiMF[[3]][[i]]$NC_MF[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_hiMF[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_hiMF[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_hiMF[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_hiMF[[3]][[i]]$NC_MH[y]<-smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_hiMF[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_hiMF[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_hiMF[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_hiMF[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_hiMF[[3]][[i]][,316:318][which(smop_hiMF[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_hiMF[[3]][[i]]$NC_SF[1:40]/smop_hiMF[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_hiMF[[3]][[i]]$NC_SF[1:40]/smop_hiMF[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_hiMF[[3]][[i]]$NC_SF[1:40]/smop_hiMF[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_hiMF[[3]][[i]]$NC_SF[1:40]/smop_hiMF[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_hiMF[[3]][[i]]$NC_SF[1:40]<-smop_hiMF[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_hiMF[[3]][[i]]$NC_SF[1:40])/sum(smop_hiMF[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_hiMF[[3]][[i]]$NC_MF[1:40]/smop_hiMF[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_hiMF[[3]][[i]]$NC_MF[1:40]/smop_hiMF[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_hiMF[[3]][[i]]$NC_MF[1:40]/smop_hiMF[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_hiMF[[3]][[i]]$NC_MF[1:40]/smop_hiMF[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_hiMF[[3]][[i]]$NC_MF[1:40]<-smop_hiMF[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_hiMF[[3]][[i]]$NC_MF[1:40])/sum(smop_hiMF[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_hiMF[[3]][[i]]$NC_MH[1:40]/smop_hiMF[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_hiMF[[3]][[i]]$NC_MH[1:40]/smop_hiMF[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_hiMF[[3]][[i]]$NC_MH[1:40]/smop_hiMF[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_hiMF[[3]][[i]]$NC_MH[1:40]/smop_hiMF[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_hiMF[[3]][[i]]$NC_MH[1:40]<-smop_hiMF[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_hiMF[[3]][[i]]$NC_MH[1:40])/sum(smop_hiMF[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_hiMF[[3]][[i]]$Tot_NC<-smop_hiMF[[3]][[i]]$NC_SF+smop_hiMF[[3]][[i]]$NC_MF+smop_hiMF[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_hiMF[[3]][[i]]$NC_SF_m2<-smop_hiMF[[3]][[i]]$NC_SF*smop_hiMF$FA.SF.2010s[i]
  smop_hiMF[[3]][[i]]$NC_MF_m2<-smop_hiMF[[3]][[i]]$NC_MF*smop_hiMF$FA.MF.2010s[i]
  smop_hiMF[[3]][[i]]$NC_MH_m2<-smop_hiMF[[3]][[i]]$NC_MH*smop_hiMF$FA.MH.2010s[i]
  smop_hiMF[[3]][[i]]$Tot_NC_m2<-smop_hiMF[[3]][[i]]$NC_SF_m2+smop_hiMF[[3]][[i]]$NC_MF_m2+smop_hiMF[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_hiMF[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_hiMF[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_hiMF[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))+
    rowSums(smop_hiMF[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_hiMF[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_hiMF[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_hiMF[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_hiMF[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))+
    rowSums(smop_hiMF[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_hiMF[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_hiMF[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_hiMF[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_hiMF[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))+
    rowSums(smop_hiMF[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_hiMF[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_hiMF[[3]][[i]]$Tot_Dem_m2<-smop_hiMF[[3]][[i]]$Dem_SF_m2+smop_hiMF[[3]][[i]]$Dem_MF_m2+smop_hiMF[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type, assuming that 2010 m2/house are constant through 2050s
  smop_hiMF[[3]][[i]]$Occ_m2_SF<-rowSums(smop_hiMF[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_hiMF[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_hiMF[[3]][[i]]$Occ_m2_MF<-rowSums(smop_hiMF[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_hiMF[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_hiMF[[3]][[i]]$Occ_m2_MH<-rowSums(smop_hiMF[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_hiMF[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_hiMF[[3]][[i]]$Occ_m2<-smop_hiMF[[3]][[i]]$Occ_m2_SF+smop_hiMF[[3]][[i]]$Occ_m2_MF+smop_hiMF[[3]][[i]]$Occ_m2_MH
  
  smop_hiMF[[3]][[i]]$m2cap_SF<-smop_hiMF[[3]][[i]]$Occ_m2_SF/smop_hiMF[[3]][[i]]$Pop_SF
  smop_hiMF[[3]][[i]]$m2cap_MF<-smop_hiMF[[3]][[i]]$Occ_m2_MF/smop_hiMF[[3]][[i]]$Pop_MF
  smop_hiMF[[3]][[i]]$m2cap_MH<-smop_hiMF[[3]][[i]]$Occ_m2_MH/smop_hiMF[[3]][[i]]$Pop_MH
  smop_hiMF[[3]][[i]]$m2cap<-smop_hiMF[[3]][[i]]$Occ_m2/smop_hiMF[[3]][[i]]$Population
}
# hi DRMF
for (i in c(1:3108)) {
  print(i)
  smop_hiDRMF[[3]][[i]]<-as.data.frame(smop_hiDRMF[[3]][[i]])
  
  smop_hiDRMF[[3]][[i]]$m2cap<-smop_hiDRMF[[3]][[i]]$m2cap_MH<-smop_hiDRMF[[3]][[i]]$m2cap_MF<- smop_hiDRMF[[3]][[i]]$m2cap_SF<-
    smop_hiDRMF[[3]][[i]]$Occ_m2<-smop_hiDRMF[[3]][[i]]$Occ_m2_MH<-smop_hiDRMF[[3]][[i]]$Occ_m2_MF<-smop_hiDRMF[[3]][[i]]$Occ_m2_SF<-
    smop_hiDRMF[[3]][[i]]$Tot_NC_m2<-smop_hiDRMF[[3]][[i]]$NC_MH_m2<-smop_hiDRMF[[3]][[i]]$NC_MF_m2<-smop_hiDRMF[[3]][[i]]$NC_SF_m2<-
    smop_hiDRMF[[3]][[i]]$Tot_Dem_m2<-smop_hiDRMF[[3]][[i]]$Dem_MH_m2<-smop_hiDRMF[[3]][[i]]$Dem_MF_m2<-smop_hiDRMF[[3]][[i]]$Dem_SF_m2<-
    smop_hiDRMF[[3]][[i]]$Tot_NC<-smop_hiDRMF[[3]][[i]]$NC_MH<-smop_hiDRMF[[3]][[i]]$NC_MF<-smop_hiDRMF[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_hiDRMF[[3]][[i]]$NC_SF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_hiDRMF[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_hiDRMF[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_hiDRMF[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MH[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_hiDRMF[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_hiDRMF[[3]][[i]]$NC_SF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_hiDRMF[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_hiDRMF[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_hiDRMF[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MH[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_hiDRMF[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_hiDRMF[[3]][[i]]$NC_SF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_hiDRMF[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_hiDRMF[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_hiDRMF[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_hiDRMF[[3]][[i]]$NC_MH[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_hiDRMF[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_hiDRMF[[3]][[i]]$NC_SF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_hiDRMF[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_hiDRMF[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MF[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_hiDRMF[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_hiDRMF[[3]][[i]]$NC_MH[y]<-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_hiDRMF[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_hiDRMF[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_hiDRMF[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_hiDRMF[[3]][[i]][,316:318][which(smop_hiDRMF[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_hiDRMF[[3]][[i]]$NC_SF[1:40]/smop_hiDRMF[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_hiDRMF[[3]][[i]]$NC_SF[1:40]/smop_hiDRMF[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_hiDRMF[[3]][[i]]$NC_SF[1:40]/smop_hiDRMF[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_hiDRMF[[3]][[i]]$NC_SF[1:40]/smop_hiDRMF[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_hiDRMF[[3]][[i]]$NC_SF[1:40]<-smop_hiDRMF[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_hiDRMF[[3]][[i]]$NC_SF[1:40])/sum(smop_hiDRMF[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_hiDRMF[[3]][[i]]$NC_MF[1:40]/smop_hiDRMF[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_hiDRMF[[3]][[i]]$NC_MF[1:40]/smop_hiDRMF[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_hiDRMF[[3]][[i]]$NC_MF[1:40]/smop_hiDRMF[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_hiDRMF[[3]][[i]]$NC_MF[1:40]/smop_hiDRMF[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_hiDRMF[[3]][[i]]$NC_MF[1:40]<-smop_hiDRMF[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_hiDRMF[[3]][[i]]$NC_MF[1:40])/sum(smop_hiDRMF[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_hiDRMF[[3]][[i]]$NC_MH[1:40]/smop_hiDRMF[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_hiDRMF[[3]][[i]]$NC_MH[1:40]/smop_hiDRMF[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_hiDRMF[[3]][[i]]$NC_MH[1:40]/smop_hiDRMF[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_hiDRMF[[3]][[i]]$NC_MH[1:40]/smop_hiDRMF[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_hiDRMF[[3]][[i]]$NC_MH[1:40]<-smop_hiDRMF[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_hiDRMF[[3]][[i]]$NC_MH[1:40])/sum(smop_hiDRMF[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_hiDRMF[[3]][[i]]$Tot_NC<-smop_hiDRMF[[3]][[i]]$NC_SF+smop_hiDRMF[[3]][[i]]$NC_MF+smop_hiDRMF[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_hiDRMF[[3]][[i]]$NC_SF_m2<-smop_hiDRMF[[3]][[i]]$NC_SF*smop_hiDRMF$FA.SF.2010s[i]
  smop_hiDRMF[[3]][[i]]$NC_MF_m2<-smop_hiDRMF[[3]][[i]]$NC_MF*smop_hiDRMF$FA.MF.2010s[i]
  smop_hiDRMF[[3]][[i]]$NC_MH_m2<-smop_hiDRMF[[3]][[i]]$NC_MH*smop_hiDRMF$FA.MH.2010s[i]
  smop_hiDRMF[[3]][[i]]$Tot_NC_m2<-smop_hiDRMF[[3]][[i]]$NC_SF_m2+smop_hiDRMF[[3]][[i]]$NC_MF_m2+smop_hiDRMF[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_hiDRMF[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_hiDRMF[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_hiDRMF[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))+
    rowSums(smop_hiDRMF[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_hiDRMF[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_hiDRMF[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_hiDRMF[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_hiDRMF[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))+
    rowSums(smop_hiDRMF[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_hiDRMF[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_hiDRMF[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_hiDRMF[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_hiDRMF[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))+
    rowSums(smop_hiDRMF[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_hiDRMF[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_hiDRMF[[3]][[i]]$Tot_Dem_m2<-smop_hiDRMF[[3]][[i]]$Dem_SF_m2+smop_hiDRMF[[3]][[i]]$Dem_MF_m2+smop_hiDRMF[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type, assuming that 2010 m2/house are constant through 2050s
  smop_hiDRMF[[3]][[i]]$Occ_m2_SF<-rowSums(smop_hiDRMF[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_hiDRMF[i,c(7,10,13,16,19,22,22,22,22,22)]),each=41),41,10))
  smop_hiDRMF[[3]][[i]]$Occ_m2_MF<-rowSums(smop_hiDRMF[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_hiDRMF[i,c(5,8,11,14,17,20,20,20,20,20)]),each=41),41,10))
  smop_hiDRMF[[3]][[i]]$Occ_m2_MH<-rowSums(smop_hiDRMF[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_hiDRMF[i,c(6,9,12,15,18,21,21,21,21,21)]),each=41),41,10))
  smop_hiDRMF[[3]][[i]]$Occ_m2<-smop_hiDRMF[[3]][[i]]$Occ_m2_SF+smop_hiDRMF[[3]][[i]]$Occ_m2_MF+smop_hiDRMF[[3]][[i]]$Occ_m2_MH
  
  smop_hiDRMF[[3]][[i]]$m2cap_SF<-smop_hiDRMF[[3]][[i]]$Occ_m2_SF/smop_hiDRMF[[3]][[i]]$Pop_SF
  smop_hiDRMF[[3]][[i]]$m2cap_MF<-smop_hiDRMF[[3]][[i]]$Occ_m2_MF/smop_hiDRMF[[3]][[i]]$Pop_MF
  smop_hiDRMF[[3]][[i]]$m2cap_MH<-smop_hiDRMF[[3]][[i]]$Occ_m2_MH/smop_hiDRMF[[3]][[i]]$Pop_MH
  smop_hiDRMF[[3]][[i]]$m2cap<-smop_hiDRMF[[3]][[i]]$Occ_m2/smop_hiDRMF[[3]][[i]]$Population
}
# RFA
for (i in c(1:3108)) { # these loops are slow, take about 4 seconds per county. Can redo the providence file here with runnung only i=2281
  print(i)
  smop_RFA[[3]][[i]]<-as.data.frame(smop_RFA[[3]][[i]])
  
  smop_RFA[[3]][[i]]$m2cap<-smop_RFA[[3]][[i]]$m2cap_MH<-smop_RFA[[3]][[i]]$m2cap_MF<- smop_RFA[[3]][[i]]$m2cap_SF<-
    smop_RFA[[3]][[i]]$Occ_m2<-smop_RFA[[3]][[i]]$Occ_m2_MH<-smop_RFA[[3]][[i]]$Occ_m2_MF<-smop_RFA[[3]][[i]]$Occ_m2_SF<-
    smop_RFA[[3]][[i]]$Tot_NC_m2<-smop_RFA[[3]][[i]]$NC_MH_m2<-smop_RFA[[3]][[i]]$NC_MF_m2<-smop_RFA[[3]][[i]]$NC_SF_m2<-
    smop_RFA[[3]][[i]]$Tot_Dem_m2<-smop_RFA[[3]][[i]]$Dem_MH_m2<-smop_RFA[[3]][[i]]$Dem_MF_m2<-smop_RFA[[3]][[i]]$Dem_SF_m2<-
    smop_RFA[[3]][[i]]$Tot_NC<-smop_RFA[[3]][[i]]$NC_MH<-smop_RFA[[3]][[i]]$NC_MF<-smop_RFA[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_RFA[[3]][[i]]$NC_SF[y]<-smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_RFA[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_RFA[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_RFA[[3]][[i]]$NC_MF[y]<-smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_RFA[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_RFA[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_RFA[[3]][[i]]$NC_MH[y]<-smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_RFA[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_RFA[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_RFA[[3]][[i]]$NC_SF[y]<-smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_RFA[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_RFA[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_RFA[[3]][[i]]$NC_MF[y]<-smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_RFA[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_RFA[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_RFA[[3]][[i]]$NC_MH[y]<-smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_RFA[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_RFA[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_RFA[[3]][[i]]$NC_SF[y]<-smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_RFA[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_RFA[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_RFA[[3]][[i]]$NC_MF[y]<-smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_RFA[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_RFA[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_RFA[[3]][[i]]$NC_MH[y]<-smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_RFA[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_RFA[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_RFA[[3]][[i]]$NC_SF[y]<-smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_RFA[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_RFA[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_RFA[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_RFA[[3]][[i]]$NC_MF[y]<-smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_RFA[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_RFA[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_RFA[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_RFA[[3]][[i]]$NC_MH[y]<-smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_RFA[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_RFA[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_RFA[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_RFA[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_RFA[[3]][[i]][,316:318][which(smop_RFA[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0 # make sure no negative inflows
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_RFA[[3]][[i]]$NC_SF[1:40]/smop_RFA[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_RFA[[3]][[i]]$NC_SF[1:40]/smop_RFA[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_RFA[[3]][[i]]$NC_SF[1:40]/smop_RFA[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_RFA[[3]][[i]]$NC_SF[1:40]/smop_RFA[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_RFA[[3]][[i]]$NC_SF[1:40]<-smop_RFA[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_RFA[[3]][[i]]$NC_SF[1:40])/sum(smop_RFA[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_RFA[[3]][[i]]$NC_MF[1:40]/smop_RFA[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_RFA[[3]][[i]]$NC_MF[1:40]/smop_RFA[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_RFA[[3]][[i]]$NC_MF[1:40]/smop_RFA[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_RFA[[3]][[i]]$NC_MF[1:40]/smop_RFA[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_RFA[[3]][[i]]$NC_MF[1:40]<-smop_RFA[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_RFA[[3]][[i]]$NC_MF[1:40])/sum(smop_RFA[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_RFA[[3]][[i]]$NC_MH[1:40]/smop_RFA[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_RFA[[3]][[i]]$NC_MH[1:40]/smop_RFA[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_RFA[[3]][[i]]$NC_MH[1:40]/smop_RFA[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_RFA[[3]][[i]]$NC_MH[1:40]/smop_RFA[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_RFA[[3]][[i]]$NC_MH[1:40]<-smop_RFA[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_RFA[[3]][[i]]$NC_MH[1:40])/sum(smop_RFA[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_RFA[[3]][[i]]$Tot_NC<-smop_RFA[[3]][[i]]$NC_SF+smop_RFA[[3]][[i]]$NC_MF+smop_RFA[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_RFA[[3]][[i]]$NC_SF_m2<-smop_RFA[[3]][[i]]$NC_SF*smop_RFA$FA.SF.NC[i]
  smop_RFA[[3]][[i]]$NC_MF_m2<-smop_RFA[[3]][[i]]$NC_MF*smop_RFA$FA.MF.NC[i]
  smop_RFA[[3]][[i]]$NC_MH_m2<-smop_RFA[[3]][[i]]$NC_MH*smop_RFA$FA.MH.NC[i]
  smop_RFA[[3]][[i]]$Tot_NC_m2<-smop_RFA[[3]][[i]]$NC_SF_m2+smop_RFA[[3]][[i]]$NC_MF_m2+smop_RFA[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_RFA[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_RFA[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))+
    rowSums(smop_RFA[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))
  smop_RFA[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_RFA[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))+
    rowSums(smop_RFA[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))
  smop_RFA[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_RFA[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))+
    rowSums(smop_RFA[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))
  smop_RFA[[3]][[i]]$Tot_Dem_m2<-smop_RFA[[3]][[i]]$Dem_SF_m2+smop_RFA[[3]][[i]]$Dem_MF_m2+smop_RFA[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type,now using update RFA assumption on floor area of new housing
  smop_RFA[[3]][[i]]$Occ_m2_SF<-rowSums(smop_RFA[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))
  smop_RFA[[3]][[i]]$Occ_m2_MF<-rowSums(smop_RFA[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))
  smop_RFA[[3]][[i]]$Occ_m2_MH<-rowSums(smop_RFA[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))
  smop_RFA[[3]][[i]]$Occ_m2<-smop_RFA[[3]][[i]]$Occ_m2_SF+smop_RFA[[3]][[i]]$Occ_m2_MF+smop_RFA[[3]][[i]]$Occ_m2_MH
  
  smop_RFA[[3]][[i]]$m2cap_SF<-smop_RFA[[3]][[i]]$Occ_m2_SF/smop_RFA[[3]][[i]]$Pop_SF
  smop_RFA[[3]][[i]]$m2cap_MF<-smop_RFA[[3]][[i]]$Occ_m2_MF/smop_RFA[[3]][[i]]$Pop_MF
  smop_RFA[[3]][[i]]$m2cap_MH<-smop_RFA[[3]][[i]]$Occ_m2_MH/smop_RFA[[3]][[i]]$Pop_MH
  smop_RFA[[3]][[i]]$m2cap<-smop_RFA[[3]][[i]]$Occ_m2/smop_RFA[[3]][[i]]$Population
}
# hiDR_RFA
for (i in c(1:3108)) { # these loops are slow, take about 4 seconds per county. Can redo the providence file here with runnung only i=2281
  print(i)
  smop_hiDR_RFA[[3]][[i]]<-as.data.frame(smop_hiDR_RFA[[3]][[i]])
  
  smop_hiDR_RFA[[3]][[i]]$m2cap<-smop_hiDR_RFA[[3]][[i]]$m2cap_MH<-smop_hiDR_RFA[[3]][[i]]$m2cap_MF<- smop_hiDR_RFA[[3]][[i]]$m2cap_SF<-
    smop_hiDR_RFA[[3]][[i]]$Occ_m2<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_MH<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_MF<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_SF<-
    smop_hiDR_RFA[[3]][[i]]$Tot_NC_m2<-smop_hiDR_RFA[[3]][[i]]$NC_MH_m2<-smop_hiDR_RFA[[3]][[i]]$NC_MF_m2<-smop_hiDR_RFA[[3]][[i]]$NC_SF_m2<-
    smop_hiDR_RFA[[3]][[i]]$Tot_Dem_m2<-smop_hiDR_RFA[[3]][[i]]$Dem_MH_m2<-smop_hiDR_RFA[[3]][[i]]$Dem_MF_m2<-smop_hiDR_RFA[[3]][[i]]$Dem_SF_m2<-
    smop_hiDR_RFA[[3]][[i]]$Tot_NC<-smop_hiDR_RFA[[3]][[i]]$NC_MH<-smop_hiDR_RFA[[3]][[i]]$NC_MF<-smop_hiDR_RFA[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_hiDR_RFA[[3]][[i]]$NC_SF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_hiDR_RFA[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MH[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_hiDR_RFA[[3]][[i]]$NC_SF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_hiDR_RFA[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MH[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_hiDR_RFA[[3]][[i]]$NC_SF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_hiDR_RFA[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_hiDR_RFA[[3]][[i]]$NC_MH[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_hiDR_RFA[[3]][[i]]$NC_SF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_hiDR_RFA[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MF[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_hiDR_RFA[[3]][[i]]$NC_MH[y]<-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_hiDR_RFA[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_hiDR_RFA[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_hiDR_RFA[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_hiDR_RFA[[3]][[i]][,316:318][which(smop_hiDR_RFA[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0 # make sure no negative inflows
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]<-smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40])/sum(smop_hiDR_RFA[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]<-smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40])/sum(smop_hiDR_RFA[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiDR_RFA[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]<-smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40])/sum(smop_hiDR_RFA[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_hiDR_RFA[[3]][[i]]$Tot_NC<-smop_hiDR_RFA[[3]][[i]]$NC_SF+smop_hiDR_RFA[[3]][[i]]$NC_MF+smop_hiDR_RFA[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_hiDR_RFA[[3]][[i]]$NC_SF_m2<-smop_hiDR_RFA[[3]][[i]]$NC_SF*smop_hiDR_RFA$FA.SF.NC[i]
  smop_hiDR_RFA[[3]][[i]]$NC_MF_m2<-smop_hiDR_RFA[[3]][[i]]$NC_MF*smop_hiDR_RFA$FA.MF.NC[i]
  smop_hiDR_RFA[[3]][[i]]$NC_MH_m2<-smop_hiDR_RFA[[3]][[i]]$NC_MH*smop_hiDR_RFA$FA.MH.NC[i]
  smop_hiDR_RFA[[3]][[i]]$Tot_NC_m2<-smop_hiDR_RFA[[3]][[i]]$NC_SF_m2+smop_hiDR_RFA[[3]][[i]]$NC_MF_m2+smop_hiDR_RFA[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_hiDR_RFA[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_hiDR_RFA[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))+
    rowSums(smop_hiDR_RFA[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))
  smop_hiDR_RFA[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_hiDR_RFA[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))+
    rowSums(smop_hiDR_RFA[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))
  smop_hiDR_RFA[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_hiDR_RFA[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))+
    rowSums(smop_hiDR_RFA[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))
  smop_hiDR_RFA[[3]][[i]]$Tot_Dem_m2<-smop_hiDR_RFA[[3]][[i]]$Dem_SF_m2+smop_hiDR_RFA[[3]][[i]]$Dem_MF_m2+smop_hiDR_RFA[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type,now using update RFA assumption on floor area of new housing
  smop_hiDR_RFA[[3]][[i]]$Occ_m2_SF<-rowSums(smop_hiDR_RFA[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))
  smop_hiDR_RFA[[3]][[i]]$Occ_m2_MF<-rowSums(smop_hiDR_RFA[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))
  smop_hiDR_RFA[[3]][[i]]$Occ_m2_MH<-rowSums(smop_hiDR_RFA[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_hiDR_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))
  smop_hiDR_RFA[[3]][[i]]$Occ_m2<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_SF+smop_hiDR_RFA[[3]][[i]]$Occ_m2_MF+smop_hiDR_RFA[[3]][[i]]$Occ_m2_MH
  
  smop_hiDR_RFA[[3]][[i]]$m2cap_SF<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_SF/smop_hiDR_RFA[[3]][[i]]$Pop_SF
  smop_hiDR_RFA[[3]][[i]]$m2cap_MF<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_MF/smop_hiDR_RFA[[3]][[i]]$Pop_MF
  smop_hiDR_RFA[[3]][[i]]$m2cap_MH<-smop_hiDR_RFA[[3]][[i]]$Occ_m2_MH/smop_hiDR_RFA[[3]][[i]]$Pop_MH
  smop_hiDR_RFA[[3]][[i]]$m2cap<-smop_hiDR_RFA[[3]][[i]]$Occ_m2/smop_hiDR_RFA[[3]][[i]]$Population
}
# hiMF_RFA
for (i in c(1:3108)) { # these loops are slow, take about 4 seconds per county. Can redo the providence file here with runnung only i=2281
  print(i)
  smop_hiMF_RFA[[3]][[i]]<-as.data.frame(smop_hiMF_RFA[[3]][[i]])
  
  smop_hiMF_RFA[[3]][[i]]$m2cap<-smop_hiMF_RFA[[3]][[i]]$m2cap_MH<-smop_hiMF_RFA[[3]][[i]]$m2cap_MF<- smop_hiMF_RFA[[3]][[i]]$m2cap_SF<-
    smop_hiMF_RFA[[3]][[i]]$Occ_m2<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_MH<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_MF<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_SF<-
    smop_hiMF_RFA[[3]][[i]]$Tot_NC_m2<-smop_hiMF_RFA[[3]][[i]]$NC_MH_m2<-smop_hiMF_RFA[[3]][[i]]$NC_MF_m2<-smop_hiMF_RFA[[3]][[i]]$NC_SF_m2<-
    smop_hiMF_RFA[[3]][[i]]$Tot_Dem_m2<-smop_hiMF_RFA[[3]][[i]]$Dem_MH_m2<-smop_hiMF_RFA[[3]][[i]]$Dem_MF_m2<-smop_hiMF_RFA[[3]][[i]]$Dem_SF_m2<-
    smop_hiMF_RFA[[3]][[i]]$Tot_NC<-smop_hiMF_RFA[[3]][[i]]$NC_MH<-smop_hiMF_RFA[[3]][[i]]$NC_MF<-smop_hiMF_RFA[[3]][[i]]$NC_SF<-0
  
  for (y in 1:10) { #2020-2029 New construction
    smop_hiMF_RFA[[3]][[i]]$NC_SF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2020_29[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2020_29[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_SF_Occ_2020_29[y]+smop_hiMF_RFA[[3]][[i]]$Dem_SF_Vac_2020_29[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2020_29[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2020_29[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MF_Occ_2020_29[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MF_Vac_2020_29[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MH[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2020_29[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2020_29[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MH_Occ_2020_29[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MH_Vac_2020_29[y]
  }
  
  for (y in 11:20) { #2030-2039 New construction
    smop_hiMF_RFA[[3]][[i]]$NC_SF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2030_39[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2030_39[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_SF_Occ_2030_39[y]+smop_hiMF_RFA[[3]][[i]]$Dem_SF_Vac_2030_39[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2030_39[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2030_39[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MF_Occ_2030_39[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MF_Vac_2030_39[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MH[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2030_39[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2030_39[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MH_Occ_2030_39[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MH_Vac_2030_39[y]
  }
  
  for (y in 21:30) { #2040-2049 New construction
    smop_hiMF_RFA[[3]][[i]]$NC_SF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2040_49[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2040_49[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_SF_Occ_2040_49[y]+smop_hiMF_RFA[[3]][[i]]$Dem_SF_Vac_2040_49[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2040_49[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2040_49[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MF_Occ_2040_49[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MF_Vac_2040_49[y]
    
    
    smop_hiMF_RFA[[3]][[i]]$NC_MH[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2040_49[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2040_49[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MH_Occ_2040_49[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MH_Vac_2040_49[y]
  }
  for (y in 31:40) { #2050-2059 New construction
    smop_hiMF_RFA[[3]][[i]]$NC_SF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Occ_2050_60[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_SF_Vac_2050_60[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_SF_Occ_2050_60[y]+smop_hiMF_RFA[[3]][[i]]$Dem_SF_Vac_2050_60[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MF[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Occ_2050_60[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MF_Vac_2050_60[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MF_Occ_2050_60[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MF_Vac_2050_60[y]
    
    smop_hiMF_RFA[[3]][[i]]$NC_MH[y]<-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y+1]+smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y+1]-
      smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Occ_2050_60[y]-smop_hiMF_RFA[[3]][[i]]$Tot_HU_MH_Vac_2050_60[y]+
      smop_hiMF_RFA[[3]][[i]]$Dem_MH_Occ_2050_60[y]+smop_hiMF_RFA[[3]][[i]]$Dem_MH_Vac_2050_60[y]
  }
  smop_hiMF_RFA[[3]][[i]][,316:318][which(smop_hiMF_RFA[[3]][[i]][,316:318]<0,arr.ind = TRUE)]<-0 # make sure no negative inflows
  # NC_SFpc[is.infinite(NC_SFpc)]<-NaN
  # correct for the extra loading of new construction towards the end of decades, by making the new construction be the same percentage of total construction in every year of each decade
  # if (any(!is.na(smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_SF[1:40]))){ 
  if (any(!is.na(smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_SF[1:40])&!is.infinite(smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_SF[1:40]))) {
    NC_SFpc<-smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_SF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_SFpc[is.infinite(NC_SFpc)]<-NaN
    NC_SFpc[which(is.nan(NC_SFpc))]<-mean(NC_SFpc,na.rm = TRUE)
    smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]<-smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc*sum(smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40])/sum(smop_hiMF_RFA[[3]][[i]]$NC_SF[1:40]*mean(NC_SFpc)/NC_SFpc) }
  
  # if (any(!is.na(smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MF[1:40]))){ 
  if (any(!is.na(smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MF[1:40])&!is.infinite(smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MF[1:40]))) {
    NC_MFpc<-smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MF[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MFpc[is.infinite(NC_MFpc)]<-NaN
    NC_MFpc[which(is.nan(NC_MFpc))]<-mean(NC_MFpc,na.rm = TRUE)
    smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]<-smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc*sum(smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40])/sum(smop_hiMF_RFA[[3]][[i]]$NC_MF[1:40]*mean(NC_MFpc)/NC_MFpc) }
  
  # if (any(!is.na(smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MH[1:40]))){ 
  if (any(!is.na(smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MH[1:40])&!is.infinite(smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MH[1:40]))) {
    NC_MHpc<-smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]/smop_hiMF_RFA[[3]][[i]]$Con_MH[1:40]
    # convert Infs to Nans, so that they can be removed in the mean function with na.rm
    NC_MHpc[is.infinite(NC_MHpc)]<-NaN
    NC_MHpc[which(is.nan(NC_MHpc))]<-mean(NC_MHpc,na.rm = TRUE)
    smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]<-smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc*sum(smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40])/sum(smop_hiMF_RFA[[3]][[i]]$NC_MH[1:40]*mean(NC_MHpc)/NC_MHpc) }
  # sum up total inflows to stock
  smop_hiMF_RFA[[3]][[i]]$Tot_NC<-smop_hiMF_RFA[[3]][[i]]$NC_SF+smop_hiMF_RFA[[3]][[i]]$NC_MF+smop_hiMF_RFA[[3]][[i]]$NC_MH
  
  # New construction floor area inflows
  smop_hiMF_RFA[[3]][[i]]$NC_SF_m2<-smop_hiMF_RFA[[3]][[i]]$NC_SF*smop_hiMF_RFA$FA.SF.NC[i]
  smop_hiMF_RFA[[3]][[i]]$NC_MF_m2<-smop_hiMF_RFA[[3]][[i]]$NC_MF*smop_hiDR_RFA$FA.MF.NC[i]
  smop_hiMF_RFA[[3]][[i]]$NC_MH_m2<-smop_hiMF_RFA[[3]][[i]]$NC_MH*smop_hiDR_RFA$FA.MH.NC[i]
  smop_hiMF_RFA[[3]][[i]]$Tot_NC_m2<-smop_hiMF_RFA[[3]][[i]]$NC_SF_m2+smop_hiMF_RFA[[3]][[i]]$NC_MF_m2+smop_hiMF_RFA[[3]][[i]]$NC_MH_m2
  # Demolition floor area outflows, using dem factors to convert stock losses to actual demolitions. Not currently calculating columns for total demolitions by type. but can sum up the groups of columns if needed to do so
  smop_hiMF_RFA[[3]][[i]]$Dem_SF_m2<-SF_dem_factor*rowSums(smop_hiMF_RFA[[3]][[i]][,248:257]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))+
    rowSums(smop_hiMF_RFA[[3]][[i]][,258:267]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))
  smop_hiMF_RFA[[3]][[i]]$Dem_MF_m2<-MF_dem_factor*rowSums(smop_hiMF_RFA[[3]][[i]][,268:277]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))+
    rowSums(smop_hiMF_RFA[[3]][[i]][,278:287]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))
  smop_hiMF_RFA[[3]][[i]]$Dem_MH_m2<-MH_dem_factor*rowSums(smop_hiMF_RFA[[3]][[i]][,288:297]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))+
    rowSums(smop_hiMF_RFA[[3]][[i]][,298:307]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))
  smop_hiMF_RFA[[3]][[i]]$Tot_Dem_m2<-smop_hiMF_RFA[[3]][[i]]$Dem_SF_m2+smop_hiMF_RFA[[3]][[i]]$Dem_MF_m2+smop_hiMF_RFA[[3]][[i]]$Dem_MH_m2
  
  # total occupied floor area stock by type,now using update RFA assumption on floor area of new housing
  smop_hiMF_RFA[[3]][[i]]$Occ_m2_SF<-rowSums(smop_hiMF_RFA[[3]][[i]][,110:119]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(7,10,13,16,19,22,25,25,25,25)]),each=41),41,10))
  smop_hiMF_RFA[[3]][[i]]$Occ_m2_MF<-rowSums(smop_hiMF_RFA[[3]][[i]][,130:139]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(5,8,11,14,17,20,23,23,23,23)]),each=41),41,10))
  smop_hiMF_RFA[[3]][[i]]$Occ_m2_MH<-rowSums(smop_hiMF_RFA[[3]][[i]][,150:159]*matrix(rep(as.numeric(smop_hiMF_RFA[i,c(6,9,12,15,18,21,24,24,24,24)]),each=41),41,10))
  smop_hiMF_RFA[[3]][[i]]$Occ_m2<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_SF+smop_hiMF_RFA[[3]][[i]]$Occ_m2_MF+smop_hiMF_RFA[[3]][[i]]$Occ_m2_MH
  
  smop_hiMF_RFA[[3]][[i]]$m2cap_SF<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_SF/smop_hiMF_RFA[[3]][[i]]$Pop_SF
  smop_hiMF_RFA[[3]][[i]]$m2cap_MF<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_MF/smop_hiMF_RFA[[3]][[i]]$Pop_MF
  smop_hiMF_RFA[[3]][[i]]$m2cap_MH<-smop_hiMF_RFA[[3]][[i]]$Occ_m2_MH/smop_hiMF_RFA[[3]][[i]]$Pop_MH
  smop_hiMF_RFA[[3]][[i]]$m2cap<-smop_hiMF_RFA[[3]][[i]]$Occ_m2/smop_hiMF_RFA[[3]][[i]]$Population
}

smop_base_FA<-smop_base
smop_hiDR_FA<-smop_hiDR
smop_hiMF_FA<-smop_hiMF
smop_hiDRMF_FA<-smop_hiDRMF
smop_RFA_FA<-smop_RFA
smop_hiDR_RFA_FA<-smop_hiDR_RFA
smop_hiMF_RFA_FA<-smop_hiMF_RFA
save(smop_base_FA,smop_hiDR_FA,smop_hiMF_FA,smop_hiDRMF_FA,smop_RFA_FA,file="HSM_results/County_FloorArea.RData") 
save(smop_base_FA,smop_hiDR_FA,smop_hiMF_FA,smop_RFA_FA,smop_hiDR_RFA_FA,smop_hiMF_RFA_FA,file="../resstock_projections/ExtData/County_FloorArea.RData")