---
title: "R (1/3) ABM Input Data - Synthetic population"
output: html_notebook
---
# LIBRARY
```{r}
# Basic Packages
library(sp) # Classes and Methods for Spatial Data.
library(rgdal)  # for vector work; sp package should always load with rgdal. ]
library(sf) # Simple Features for R, st_as_sf
library (raster)   # for metadata/attributes- vectors or rasters
library(dplyr)    # mutate, left_join
library(tidyverse) # tidy data 
library(data.table) # setnames
library(plyr)
library(igraph) # network analysis
library(tmap)
library(spatstat) # analysis of spatial point patterns,  ‘ppp’ (point pattern) object
library(maptools) # coerce SpatialPolygons into an ‘owin: observ window’ object
library(spdep) # local moran's I, nb2listw, 
library(rgeos)
library(magrittr)  # %>%
library(reshape2)  # dcast, cast category variables to 1/0
# STATISTIC
library(limSolve) # https://search.r-project.org/CRAN/refmans/limSolve/html/lsei.html 
library(GGally)       #correlation diagram  
# PLT
library(ggplot2)
library(corrplot)
library(htmlwidgets)
library(shinyjs)
library(contoureR)
library(geometry)
library(tmaptools)
library(RColorBrewer) # for color selection
library(hrbrthemes) # theme
library(ggbreak)     # break x/y axis     
library(ggpubr)
library(viridis)

# NETWORK 
library(igraph)
```

# Data Collection 
Study area: Chautauqua County 
Data: 2021 
    Tax Parcel data; 
    road network; schools; 
    2021 census data
          B01001 SEX BY AGE (2020 5yr)
            https://data.census.gov/cedsci/table?q=population&g=0500000US36013%241500000&tid=ACSDT5Y2020.B01001 
          B14002  Sex by school enrollment by level of school by type of school for the population 3 years and over
            https://data.census.gov/cedsci/table?q=education&g=0500000US36013%241500000&tid=ACSDT5Y2020.B14002
          B08007  Sex of workers by place of work--state and county level
            https://data.census.gov/cedsci/table?q=employment,%20age&g=0500000US36013%241500000&tid=ACSDT5Y2020.B08007
          B01002 MEDIAN AGE BY SEX (2020 5yr)
            https://data.census.gov/cedsci/table?q=population&g=0500000US36013%241500000&tid=ACSDT5Y2020.B01002
          B09019 HOUSEHOLD TYPE (INCLUDING LIVING ALONE) BY RELATIONSHIP (2020 5yr)
            https://data.census.gov/table?q=B09019&g=0500000US36013$1500000
          B09021  LIVING ARRANGEMENTS OF ADULTS 18 YEARS AND OVER BY AGE
            https://data.census.gov/cedsci/table?q=B09021&g=0500000US36013%241500000&tid=ACSDT5Y2020.B09021
          B11005 HOUSEHOLDS BY PRESENCE OF PEOPLE UNDER 18 YEARS BY HOUSEHOLD TYPE (2020 5yr)
            https://data.census.gov/cedsci/table?q=children&g=0500000US36013%241500000&tid=ACSDT5Y2020.B11005
          B09021  LIVING ARRANGEMENTS OF ADULTS 18 YEARS AND OVER BY AGE
            https://data.census.gov/cedsci/table?q=B09021&g=0500000US36013%241500000 
          B11004 FAMILY TYPE BY PRESENCE AND AGE OF RELATED CHILDREN UNDER 18 YEARS (2020 5yr)
            https://data.census.gov/cedsci/table?q=children&g=0500000US36013%241500000&tid=ACSDT5Y2020.B11004
          B09020  RELATIONSHIP BY HOUSEHOLD TYPE (INCLUDING LIVING ALONE) FOR THE POPULATION 65 YEARS AND OVER
          https://data.census.gov/cedsci/table?q=Household%20type%20child&g=0500000US36013%241500000&tid=ACSDT5Y2020.B09020
          B11001  HOUSEHOLD TYPE (INCLUDING LIVING ALONE)
          https://data.census.gov/cedsci/table?q=Household%20type&g=0500000US36013%241500000&tid=ACSDT5Y2020.B11001
          B11012  HOUSEHOLDS BY TYPE
          https://data.census.gov/cedsci/table?q=Household%20type&g=0500000US36013%241500000&tid=ACSDT5Y2020.B11012
          
          P5  Group Quarters Population By Major Group Quarters Type (2020 Decennial)
          https://data.census.gov/cedsci/table?q=Group%20Quarters%20Population&g=0500000US36013%241500000
          
          2020 Urban Area Clusters 
          https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2020&layergroup=Urban+Areas
          
          
          # UNUSED 
          B11012  HOUSEHOLDS BY TYPE
          https://data.census.gov/cedsci/table?q=Household%20type&g=0500000US36013%241500000&tid=ACSDT5Y2020.B11012
          B11016  HOUSEHOLD TYPE BY HOUSEHOLD SIZE
          https://data.census.gov/cedsci/table?q=Household%20type&g=0500000US36013%241500000&tid=ACSDT5Y2020.B11016
          
          P5  Group Quarters Population By Major Group Quarters Type (2020 Decennial)
          https://data.census.gov/cedsci/table?q=Group%20Quarters%20Population&g=0500000US36013%241500000

          
    U.S. Census Bureau’s Longitudinal Employer-Household Dynamics (LEHD) Origin-Destination Employment Statistics (LODES)
    0-18 School locations
    18-22 University

Chautauqua: 
    COUNTYFP: "013"
    COUNTYNS: "00974105"
    GEOID: "36013"

## Spatial Data
Projected CRS: NAD83 / UTM zone 18N
```{r}
# County
tt <- st_read("data/tl_2021_us_county/tl_2021_us_county.shp")
c_county <- tt[tt$COUNTYNS=="00974105",]

# Tract
tt <- st_read("data/tl_2021_36_tract/tl_2021_36_tract.shp")
c_tract <- tt[tt$COUNTYFP=="013",]

# Block Groups 
tt <- st_read("data/tl_2021_36_bg/tl_2021_36_bg.shp")
c_bgroup <- tt[tt$COUNTYFP=="013",]

# Parcel & Residential Parcel 
c_parcel <- st_read("data/Chautauqua_2021_Tax_Parcels_SHP_2203/Chautauqua_2021_Tax_Parcels_SHP_2203.shp")
c_parcel_r <- c_parcel[(c_parcel$PROP_CLASS > "199" & c_parcel$PROP_CLASS < "300"),] 
c_parcel_r <- c_parcel_r[!is.na(c_parcel_r$COUNTY),]

# Group Quarter Parcel: NAD83 / UTM zone 18N
c_parcel_gq <- c_parcel[c_parcel$PROP_CLASS %in% c("613","615","633","641","642","660","670"),] 

# Urban Areas 2020
c_urban <- st_read(paste(dirname(dirname(getwd())), "/zz_NYS_GIS_Data/tl_2020_us_uac10/tl_2020_us_uac10.shp", sep="")) %>% filter(NAME10 %like% ", NY")

# Project spatial data (e.g., county, tract, block groups) to "NAD83 / UTM zone 18N"
c_county <- st_transform(c_county, st_crs(c_parcel_r))
c_tract <- st_transform(c_tract, st_crs(c_parcel_r))
c_bgroup <- st_transform(c_bgroup, st_crs(c_parcel_r))
c_urban <- st_transform(c_urban, st_crs(c_parcel_r))

# Urban Areas only in Chautauqua County
c_urban <- c_urban[lengths(st_intersects(c_urban,c_county))>0,]
```

## Census (css) Data  
### Individual Level 
```{r}
# B01001 SEX BY AGE (2020 5yr)
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B01001-2022-05-19T145839.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
css_age_gender <- ttt

# B08007  Sex of workers by place of work--state and county level
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B08007-2022-05-23T034722.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
css_work_place <- ttt[c("GEOID","NAMELSAD","TRACT_NAMELSAD","m_ttl","m_Worked_in_county_of_residence",
                        "m_Worked_outside_county_of_residence","m_Worked_outside_state_of_residence",
                        "f_ttl","f_Worked_in_county_of_residence",
                        "f_Worked_outside_county_of_residence","f_Worked_outside_state_of_residence")]

# B09021  Living arrangements of adults 18 years and over by age
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B09021-2022-05-24T015931.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
css_over18_in_hh <- ttt
css_over18_in_hh$hhder_spouse_partner_ttl_over17 <- css_over18_in_hh$married_hhder_spouse_ttl_over17 +  css_over18_in_hh$unmarried_hhder_partner_ttl_over17
css_over18_in_hh$hhder_spouse_partner_18_34 <- css_over18_in_hh$married_hhder_spouse_18_34 + css_over18_in_hh$unmarried_hhder_partner_18_34
css_over18_in_hh$hhder_spouse_partner_35_64 <- css_over18_in_hh$married_hhder_spouse_35_64 + css_over18_in_hh$unmarried_hhder_partner_35_64
css_over18_in_hh$hhder_spouse_partner_over64 <- css_over18_in_hh$married_hhder_spouse_over64 + css_over18_in_hh$unmarried_hhder_partner_over64
css_over18_in_hh <- css_over18_in_hh[,!colnames(css_over18_in_hh) %in% c("married_hhder_spouse_ttl_over17","unmarried_hhder_partner_ttl_over17","married_hhder_spouse_18_34","unmarried_hhder_partner_18_34","married_hhder_spouse_35_64", "unmarried_hhder_partner_35_64","married_hhder_spouse_over64","unmarried_hhder_partner_over64")]

# P5  Group Quarters Population By Major Group Quarters Type (2020 Decennial)
lt <- c("Institutionalized.population","Other.institutional.facilities",
        "Noninstitutionalized.population","Military.quarters",
        "in_groupquarter","temp")
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_DECENNIALPL2020.P5-2022-06-06T011200.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
ttt <- ttt[, !names(ttt)%in% lt]
ttt[,6:length(ttt)] <- round(ttt[,6:length(ttt)]/ttt$Total,4)
ttt[ttt == "NaN"] <- 0
ttt$temp <- apply(ttt[,6:length(ttt)], 1, sum)
ttt[ttt$temp ==0,]$Other.noninstitutional.facilities <- 1

ttt <- left_join(ttt, css_hh_type[c("GEOID","in_groupquarter")])
ttt$Total <- ttt$in_groupquarter
ttt[,6:length(ttt)] <- round(ttt[,6:length(ttt)] * ttt$Total,0)
ttt <- ttt[, !names(ttt)%in% lt]
css_group_quarter_pop <- ttt

# B09020  RELATIONSHIP BY HOUSEHOLD TYPE (INCLUDING LIVING ALONE) FOR THE POPULATION 65 YEARS AND OVER
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B09020-2022-06-13T143729.csv", 
                colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
# combine parent & parent-in-law
ttt$over64_parent_or_inlaw <- ttt$over64_parent + ttt$over64_parent.in.law
# combine (family/non-family) non-relatives with relatives
ttt$over64_relatives_nonrelatives <- ttt$over64_other_relatives + ttt$over64_nonrelatives + ttt$over64_nonfamily_nonrelatives
# combine family householder with non-family not-alone householder 
ttt$over64_m_hhder_not_alone <- ttt$over64_m_hhder + ttt$over64_nonfamily_m_hhder_not_alone
ttt$over64_f_hhder_not_alone <- ttt$over64_f_hhder + ttt$over64_nonfamily_f_hhder_not_alone

lt <- c("over64_in_hh","over64_in_hh_family","over64_in_hh_nonfamily",
  "over64_parent","over64_parent.in.law",
  "over64_other_relatives","over64_nonrelatives","over64_nonfamily_nonrelatives",
  "over64_m_hhder","over64_nonfamily_m_hhder_not_alone",
  "over64_f_hhder","over64_nonfamily_f_hhder_not_alone")
ttt <- ttt[,!colnames(ttt) %in% lt]
css_over64_in_hh <- ttt

# ------------ Backup -----------
# # B14002  Sex by school enrollment by level of school by type of school for the population 3 years and over
# ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B14002-2022-05-23T033555.csv", colClasses = "character")
# ttt <- cln_census_geoid(ttt, idx_cbgroup)
# ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
# css_school_enroll <- ttt
```



### Household level
```{r}
# B09019 household type (including living alone) by relationship 
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B09019-2022-05-19T155332.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
css_hh_type <- ttt
css_hh_type$hhmber_Opposite_sex <- css_hh_type$hhmber_Opposite_sex_spouse + css_hh_type$hhmber_Opposite_sex_unmarried_partner
css_hh_type$hhmber_Same_sex <- css_hh_type$hhmber_Same_sex_spouse + css_hh_type$hhmber_Same_sex_unmarried_partner
css_hh_type$hhmber_Child <- css_hh_type$hhmber_Child + css_hh_type$hhmber_Foster_child
css_hh_type$hhmber_Other <- css_hh_type$hhmber_Other_nonrelatives + css_hh_type$hhmber_Other_relatives
css_hh_type <- css_hh_type[, !colnames(css_hh_type) %in% c("hhmber_Opposite_sex_spouse","hhmber_Opposite_sex_unmarried_partner",
                                            "hhmber_Same_sex_spouse", "hhmber_Same_sex_unmarried_partner","hhmber_Foster_child", 
                                            "hhmber_Other_nonrelatives", "hhmber_Other_relatives")]

# B11012  HOUSEHOLDS BY TYPE - married + cohabiting. This df helps to identify married & cohabiting hh. Ignore the udr18 information.  
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B11012-2022-06-20T211726.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
css_hh_type_2 <- ttt

# B11005 households by presence of people under 18 years by household type
ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B11005-2022-05-19T162036.csv", colClasses = "character")
ttt <- cln_census_geoid(ttt, idx_cbgroup)
ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
css_hh_with_udr18 <- ttt


# ------------ Backup -----------
# # B11004 family type by presence and age of related children under 18 years  
# ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B11004-2022-05-19T162839.csv", colClasses = "character")
# ttt <- cln_census_geoid(ttt, idx_cbgroup)
# ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
# css_hh_with_udr18_brkdown <- ttt

# # B11001  HOUSEHOLD TYPE (INCLUDING LIVING ALONE)
# ttt <- read.csv("data/ACSDT5Y2020_chau_block_group_cleaned/Chautanqua_bg_ACSDT5Y2020.B11001-2022-06-20T172331.csv", colClasses = "character")
# ttt <- cln_census_geoid(ttt, idx_cbgroup)
# ttt[,5:length(ttt)] <- cln_census_numeric(ttt[,5:length(ttt)])
# css_hh_type_2 <- ttt
```


## Reorganize household type 
```{r}
# 1: married family and cohabiting couple household (married_cohabiting) 
data <- css_hh_type_2[,c(1,5:21)]
## new columns 
data$hh_married_cohabiting <- data$hh_family_married + data$hh_cohabiting_couple
data$hh_married_cohabiting_udr18 <- data$hh_family_married_with_udr18 + data$hh_cohabiting_couple_with_udr18
data$hh_married_cohabiting_NO18 <- data$hh_family_married_NO_udr18 + data$hh_cohabiting_couple_NO_udr18

# join two household level data 
data <- left_join(data, css_hh_type[c("GEOID","hhder_ttl","hhder_m_ttl","hhder_f_ttl")])

# new columns 
data$hh_f_hhder_married_cohabiting <- data$hhder_f_ttl - data$hh_f_hhder_nospouse
data$hh_m_hhder_married_cohabiting <- data$hhder_m_ttl - data$hh_m_hhder_nospouse

# Validate
data$hh_family_married + data$hh_cohabiting_couple - 
  data$hh_f_hhder_married_cohabiting - data$hh_m_hhder_married_cohabiting

# Final household level dataset 
data <- data[c("GEOID","hh_ttl","hhder_m_ttl","hhder_f_ttl",
               "hh_m_hhder_married_cohabiting","hh_f_hhder_married_cohabiting",
               "hh_m_hhder_nospouse","hh_f_hhder_nospouse",
               "hh_married_cohabiting","hh_married_cohabiting_udr18","hh_married_cohabiting_NO18",
               
               "hh_m_hhder_nospouse_alone","hh_m_hhder_nospouse_with_udr18",
               "hh_m_hhder_nospouse_with_relatives","hh_m_hhder_nospouse_non_relatives",
               
               "hh_f_hhder_nospouse_alone","hh_f_hhder_nospouse_with_udr18",
               "hh_f_hhder_nospouse_with_relatives","hh_f_hhder_nospouse_non_relatives")]


# 2: Householder Alone >= 65
tt <- css_over64_in_hh[c("GEOID","over64_nonfamily_m_hhder_alone","over64_nonfamily_f_hhder_alone")]
setnames(tt, old=c("over64_nonfamily_m_hhder_alone","over64_nonfamily_f_hhder_alone"),
         new = c("hh_m_hhder_nospouse_alone_over64","hh_f_hhder_nospouse_alone_over64"))
data <- left_join(data, tt)

# 3: Classify household into 12 types 
## -------- Married or Cohabiting Couple ------
data$h1_mc_udr18 <- data$hh_married_cohabiting_udr18
data$h2_mc_NO18 <- data$hh_married_cohabiting_NO18

## -------- Male householder no spouse ------
data$h3_mns_alone_less65 <- data$hh_m_hhder_nospouse_alone - data$hh_m_hhder_nospouse_alone_over64
data$h4_mns_alone_over64 <- data$hh_m_hhder_nospouse_alone_over64
data$h5_mns_udr18 <- data$hh_m_hhder_nospouse_with_udr18
data$h6_mns_NO18_others <- data$hh_m_hhder_nospouse_with_relatives 
data$h7_mns_nonfamily <- data$hh_m_hhder_nospouse_non_relatives

## -------- Female householder no spouse ------
data$h8_fns_alone_less65 <- data$hh_f_hhder_nospouse_alone - data$hh_f_hhder_nospouse_alone_over64
data$h9_fns_alone_over64 <- data$hh_f_hhder_nospouse_alone_over64
data$h10_fns_udr18 <- data$hh_f_hhder_nospouse_with_udr18
data$h11_fns_NO18_others <- data$hh_f_hhder_nospouse_with_relatives 
data$h12_fns_nonfamily <- data$hh_f_hhder_nospouse_non_relatives

css_hh_type_fnl <- data

# 4： Group Quarters, pop >= 65
data <- css_over64_in_hh[c("GEOID","over64_in_gqquarter")]
data <- left_join(css_group_quarter_pop, data)
data$age_13_18_juvenilefacilities <- data$Juvenile.facilities
data$age_15_24_college <- data$College.University.student.housing
data$age_18_64_other <- data$Total - data$over64_in_gqquarter - data$age_13_18_juvenilefacilities - data$age_15_24_college
css_group_quarter_pop_2 <- data
```
# SET UP
```{r}
# Random Seed
rm_seed <- 315
set.seed(rm_seed)
# prepare index 
lt_geoid_bg <- idx_cbgroup$GEOID
lt_geoid_bg <- lt_geoid_bg[lt_geoid_bg!="360139900000"]
```

# Generate Individual Population
```{r}
# set random seed 
set.seed(rm_seed)

dfs_ind <- list()
for (i in 1:length(lt_geoid_bg)) {
  # create an empty dataframe for individuals in census bg i  
  tt <- css_age_gender[css_age_gender$GEOID == lt_geoid_bg[[i]],]$ttl 
  df <- setnames(data.frame(matrix(nrow = tt, ncol = 5)), new = c("ind_id","gender","age","hh_id","hh_role"))
  
  # unique identifier & gender & age
  df$ind_id <- paste(lt_geoid_bg[[i]], "_ind_",rownames(df), sep = "")
  df$gender <- rep(c("Male","Female"), times = c(css_age_gender[css_age_gender$GEOID == lt_geoid_bg[[i]],]$m_ttl,
                                                 css_age_gender[css_age_gender$GEOID == lt_geoid_bg[[i]],]$f_ttl))
  df$age <- ind_age_generator(lt_geoid_bg[[i]])
  
  # append dataframe to the list 
  dfs_ind[[length(dfs_ind) + 1]] <- df
}

names(dfs_ind) <- lt_geoid_bg
```

# Assign: Group Quarter
```{r}
# set random seed 
set.seed(rm_seed)

# Prepare data frame 
foo <- dfs_ind   # individual 
dfs_gquarter <- list()

for (i in 1:length(lt_geoid_bg)) {
  # Start Here
  idx <- lt_geoid_bg[[i]]
  m <- css_group_quarter_pop_2 %>% filter(GEOID == idx) %>%
    select(c("Total","age_13_18_juvenilefacilities","age_15_24_college","age_18_64_other","over64_in_gqquarter")) %>% 
    unlist() %>% as.matrix()
  m <- m + 0.01
  
  # Create an empty dataframe to store group quarter populations
  df <- data.frame(matrix(nrow = 0, ncol = 2)) %>% setnames(new = c("ind_id","type"))
  
  if(m[1] > 0){
    data <- foo[[idx]]
    
    # College Student 
    a <- data %>% filter(is.na(hh_role) & age %in% 15:24) %>% slice_sample( n = m[3]) %>% pull(ind_id)
    data[data$ind_id %in% a,"hh_role"] <- "in_gq"
    data[data$ind_id %in% a,"hh_id"] <- "gq_college_housing"
    
    # Juvenile_facility
    b <- data %>% filter(is.na(hh_role) & age %in% 13:18) %>% slice_sample( n = m[2]) %>% pull(ind_id)
    data[data$ind_id %in% b,"hh_role"] <- "in_gq"
    data[data$ind_id %in% b,"hh_id"] <- "gq_Juvenile_facility"
    
    # Other adult 
    c <- data %>% filter(is.na(hh_role) & age %in% 18:64) %>% slice_sample( n = m[4]) %>% pull(ind_id)
    data[data$ind_id %in% c,"hh_role"] <- "in_gq"
    data[data$ind_id %in% c,"hh_id"] <- "gq_others_adult"
    
    # Other for over 64
    d <- data %>% filter(is.na(hh_role) & age > 64) %>% slice_sample( n = m[5]) %>% pull(ind_id)
    data[data$ind_id %in% d,"hh_role"] <- "in_gq"
    data[data$ind_id %in% d,"hh_id"] <- "gq_others_adult_over64"
  }
  
  # Collect all individuals in gq 
  df <- data %>% filter(hh_role == "in_gq")
  dfs_gquarter[[length(dfs_gquarter)  + 1]] <- df
  
  foo[[idx]] <- data
  print(i)
}

names(dfs_gquarter) <- lt_geoid_bg
dfs_ind <- foo
```

# HOUSEHOLD 
## Create empty HH container for 12 types 
```{r}
# set random seed 
set.seed(rm_seed)

# Table: css_hh_type_fnl
foo <- dfs_ind   # individual 
dfs_hh <- list()

lt <- names(css_hh_type_fnl[,22:33])

# Start Loop to create household df containers
for (i in 1:length(lt_geoid_bg)) {
  # start here, get census block group index
  idx <- lt_geoid_bg[[i]]
  m <- css_hh_type_fnl %>% filter(GEOID == idx) %>% select(c(22:33)) %>% unlist()
  
  df <- data.frame(hh_id = NA,type = rep(lt, m), if_pop = 0)
  df$hh_id <- paste(idx, "_hh_",rownames(df), sep = "")
  dfs_hh[[1+length(dfs_hh)]] <- df
}

names(dfs_hh) <- lt_geoid_bg
```

## Fill Type 1: h1_mc_udr18
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# Initialize 
child_mom_brks <- c(0,17,41,100)
mom_dad_brks <- c(-100,-5,9,100)
labls <- c(0,1,0)

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # household to fill 
  lt_temp <- df %>% filter(type == "h1_mc_udr18") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age < 18) %>% nrow()
  n <- min(length(lt_temp), n)
  
  if(n > 0){
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]
      
      # While loop start 
      oo = 0
      while (oo == 0) {
        # Randomly select one child < 18
        a <- data %>% filter(is.na(hh_id) & age < 18) %>% slice_sample(n=1)
        b <- data %>% filter(is.na(hh_id) & gender == "Female" & age >=18) %>% slice_sample(n=1)
        c <- data %>% filter(is.na(hh_id) & gender == "Male" & age >=18) %>% slice_sample(n=1)
        
        o1 <- as.numeric(as.character(cut(b$age - a$age, breaks = child_mom_brks, labels = labls)))
        o2 <- as.numeric(as.character(cut(c$age - b$age, breaks = mom_dad_brks, labels = labls)))
        oo <- o1 * o2
      }
      
      # Save Result
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"child_udr18")
      data[data$ind_id == b$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"mom")
      data[data$ind_id == c$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"dad")
      
      print(j)
  }
    
  }
  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}

dfs_ind <- foo    


# # save the result 
# # imp_dfs_ind <- dfs_ind
# data <- export_imp_dfs(imp_dfs_ind)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step1_populate_type1_h1_mc_udr18.csv", row.names = F)
```

## Fill Type 5&10: single parent
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# Initialize 
child_mom_brks <- c(0,17,41,100)
child_dad_brks <- c(0,17,49,100)
labls <- c(0,1,0)

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  
  # type 5: h5_mns_udr18   household to fill 
  lt_temp <- df %>% filter(type == "h5_mns_udr18") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age < 18) %>% nrow()
  n <- min(n, length(lt_temp))
  
  if(n>0){
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]
      oo=0   # while loop stop condition oo=1
      n_loop = 1000
      while (oo==0 & n_loop > 0) {
        a <- data %>% filter(age < 18 & is.na(hh_id)) %>% slice_sample(n=1)
        b <- data %>% filter(gender=="Male" & age >=18 & is.na(hh_id)) %>% slice_sample(n=1)
        oo <- as.numeric(as.character(cut(b$age - a$age, breaks = child_dad_brks, label = labls)))
        n_loop = n_loop-1 }
      
      if(n_loop <=0){break}
      # Save Result
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"child_udr18")
      data[data$ind_id == b$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"single_dad")
      print(paste("single_dad   ",j, sep = "  "))
      
    }  }
  
  # type 10: h10_fns_udr18 
  lt_temp <- df %>% filter(type == "h10_fns_udr18") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age < 18) %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]            # household id
      oo=0
      n_loop = 1000
      while (oo==0 & n_loop > 0) {
        a <- data %>% filter(is.na(hh_role) & age < 18) %>% slice_sample(n=1)
        c <- data %>% filter(is.na(hh_role) & age >= 18 & gender=="Female") %>% slice_sample(n=1)
        oo <- as.numeric(as.character(cut(c$age - a$age, breaks = child_mom_brks, label = labls)))
        n_loop = n_loop-1
      }
      
      if(n_loop <=0){break}
      
      # Save Result
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"child_udr18")
      data[data$ind_id == c$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"single_mom")
      print(paste("single_mom   ",j, sep = "  "))
    }
  }

  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
  
}

# # save result 
# data <- export_imp_dfs(dfs_ind)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step2_populate_type5_10_single_parent.csv", row.names = F)

# dfs_ind <- foo
```

## Fill Type 4&9: alone >64
h4_mns_alone_over64
h9_fns_alone_over64
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# Populate 
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # ------------- MALE ALONE > 64
  # number of hh with older MALE singleton
  lt_temp <- df %>% filter(type=="h4_mns_alone_over64") %>% pull(hh_id)
  # count of individuals satisfying constraints 
  n <- data %>% filter(is.na(hh_role) & gender =="Male" & age > 64) %>% nrow()
  n <- min(length(lt_temp), n)
  
  if(n>0)(
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]            # household id
      a <- data %>% filter(is.na(hh_role) & gender =="Male" & age > 64) %>% slice_sample(n=1)
      # Save Result
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"male_oldersingleton")
    })
  
  # ------------- FEMALE ALONE > 64
  # number of hh with older FEMALE singleton
  lt_temp <- df %>% filter(type=="h9_fns_alone_over64") %>% pull(hh_id)
  # count of individuals satisfying constraints 
  n <- data %>% filter(is.na(hh_role) & gender =="Female" & age > 64) %>% nrow()
  n <- min(length(lt_temp), n)
  
  if(n>0)(
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]            # household id
      a <- data %>% filter(is.na(hh_role) & gender =="Female" & age > 64) %>% slice_sample(n=1)
      # Save Result
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"female_oldersingleton")
    })

  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}


# # # save result 
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step3_populate_type4_9_older_singleton.csv", row.names = F)

# dfs_ind <- foo

```

## Fill Type 2: h2_mc_NO18
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# Initialize
mom_dad_brks <- c(-100,-5,9,100)
labls <- c(0,1,0)

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # ------------- fill type 2: h2_mc_NO18
  # decide the minimum number of hh to fill
  lt_temp <- df %>% filter(type == "h2_mc_NO18") %>% pull(hh_id)
  n <- min(data %>% filter(is.na(hh_role) & age >= 18 & gender=="Male") %>% nrow(),
        data %>% filter(is.na(hh_role) & age >= 18 & gender=="Female") %>% nrow())
  n <- min(length(lt_temp), n)
  
  if(n>0)(
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]
      # Define stop condition 
      oo=0
      n_loop = 1000  # maximum runs
      
      while (oo ==0 & n_loop >0) {
        a <- data %>% filter(is.na(hh_role) & age >=18 & gender=="Male") %>% slice_sample(n=1)
        b <- data %>% filter(is.na(hh_role) & age >=18 & gender=="Female") %>% slice_sample(n=1)
        oo <- as.numeric(as.character(cut(a$age - b$age, breaks = mom_dad_brks, label = labls)))
        n_loop = n_loop -1
      }
      if(n_loop <=0){break}
      # Save Result
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"husband_NO_udr18")
      data[data$ind_id == b$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"wife_NO_udr18")
      print(paste("mc_NO_18",j, sep = "  "))
    }
  )

  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}

# # # save result 
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step4_populate_type2_h2_mc_NO18.csv", row.names = F)

# dfs_ind <- foo
```
```{r}
data <- read.csv("data/00_saved_dataset/01_step4_populate_type2_h2_mc_NO18.csv")
```

## Fill Type 3&8: Adult alone < 65
h3_mns_alone_less65
h8_fns_alone_less65
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # type 3: h3_mns_alone_less65
  lt_temp <- df %>% filter(type == "h3_mns_alone_less65") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >=18 & age <65 & gender =="Male") %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    a <- data %>% filter(is.na(hh_role) & age >=18 & age <65 & gender =="Male") %>% slice_sample(n=n) %>% pull("ind_id")
    data[data$ind_id %in% a,]$hh_id <- lt_temp[1:n]
    data[data$ind_id %in% a,]$hh_role <- "mns_alone_less65"
  }
  
  # type 8: h8_fns_alone_less65
  lt_temp <- df %>% filter(type == "h8_fns_alone_less65") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >=18 & age <65 & gender =="Female") %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    b <- data %>% filter(is.na(hh_role) & age >=18 & age <65 & gender =="Female") %>% slice_sample(n=n) %>% pull("ind_id")
    data[data$ind_id %in% b,]$hh_id <- lt_temp[1:n]
    data[data$ind_id %in% b,]$hh_role <- "fns_alone_less65"
  }

  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}

# save result
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step5_populate_type3_8_alone_less65.csv")
# 
# dfs_ind <- foo

```


## Fill Type 6&11: Adult not alone with other relatives
h6_mns_NO18_others
h11_fns_NO18_others
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# populate adult 
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  
  # Type 6: h6_mns_NO18_others
  lt_temp <- df %>% filter(type == "h6_mns_NO18_others") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >= 18 & gender=="Male") %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    a <- data %>% filter(is.na(hh_role) & age >= 18 & gender=="Male") %>% slice_sample(n=n) %>% pull(ind_id)
    data[data$ind_id %in% a,]$hh_id <- lt_temp[1:n]
    data[data$ind_id %in% a,]$hh_role <- "mns_no18_other_relative"
  }

  # Type 11: h11_fns_NO18_others
  lt_temp <- df %>% filter(type == "h11_fns_NO18_others") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >= 18 & gender=="Female") %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    b <- data %>% filter(is.na(hh_role) & age >= 18 & gender=="Female") %>% slice_sample(n=n) %>% pull(ind_id)
    data[data$ind_id %in% b,]$hh_id <- lt_temp[1:n]
    data[data$ind_id %in% b,]$hh_role <- "fns_no18_other_relative"
  }

  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}

# save result
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step6_populate_type6_11_adult_hhder_NOsp_NO18_otherrelatives.csv")
# dfs_ind <- foo

# imp_dfs_ind <- dfs_ind
# imp_dfs_hh <- dfs_hh
```

## Fill Child >= 18 
### Get child age limit 
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- imp_dfs_ind   # individual 
tt <- imp_dfs_hh  # household 

# Find the the lower and higher bound of child age
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # -------------

  # MALE householder -    potential father  
  a <- data %>% filter(hh_role %in% c("dad","husband_NO_udr18","single_dad","mns_no18_other_relative"))
  # FEMALE householder -  potential mother  
  b <- data %>% filter(hh_role %in% c("mom","wife_NO_udr18","single_mom","fns_no18_other_relative"))
  
  df <- left_join(df,a[c("hh_id","age")], by=c("hh_id"="hh_id"))
  df <- setnames(df, old="age",new="male_hhder_age")
  df <- left_join(df,b[c("hh_id","age")], by=c("hh_id"="hh_id"))
  df <- setnames(df, old="age",new="female_hhder_age")
  
  # Child age limit 
  df$child_lo_b <- 999
  df$child_hi_b <- -1
  df[!is.na(df$female_hhder_age),]$child_lo_b <- df[!is.na(df$female_hhder_age),]$female_hhder_age - 40
  df[!is.na(df$female_hhder_age),]$child_hi_b <- df[!is.na(df$female_hhder_age),]$female_hhder_age - 18
  
  
  df[(!is.na(df$male_hhder_age) & !is.na(df$female_hhder_age)),]$child_lo_b <- mapply(max, 
                                                                                      df[(!is.na(df$male_hhder_age) & !is.na(df$female_hhder_age)),]$male_hhder_age - 49, 
                                                                                      df[(!is.na(df$male_hhder_age) & !is.na(df$female_hhder_age)),]$child_lo_b)
  df[(!is.na(df$male_hhder_age) & !is.na(df$female_hhder_age)),]$child_hi_b <- mapply(min, 
                                                                                      df[(!is.na(df$male_hhder_age) & !is.na(df$female_hhder_age)),]$male_hhder_age - 18, 
                                                                                      df[(!is.na(df$male_hhder_age) & !is.na(df$female_hhder_age)),]$child_hi_b)
  
  df[(!is.na(df$male_hhder_age) & is.na(df$female_hhder_age)),]$child_lo_b <- df[(!is.na(df$male_hhder_age) & is.na(df$female_hhder_age)),]$male_hhder_age - 49
  df[(!is.na(df$male_hhder_age) & is.na(df$female_hhder_age)),]$child_hi_b <- df[(!is.na(df$male_hhder_age) & is.na(df$female_hhder_age)),]$male_hhder_age - 18

  tt[[idx]] <- df
  # print(paste(i,"done",sep = "----------------"))
}
```

### Assign child >= 18
```{r}
# set random seed 
set.seed(rm_seed)

# Use previous tables: foo (individual) & tt (household)
# Populate by census block group

lt_hhtype_temp <- c("h1_mc_udr18","h2_mc_NO18","h5_mns_udr18","h6_mns_NO18_others","h10_fns_udr18","h11_fns_NO18_others")
lt_age_range <- list(18:34,
                     35:64,
                     65:100)
labls <- c(0,1,0)

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # Census data 
  m1 <- css_over18_in_hh %>% filter(GEOID == idx) %>% select("child_hhmber_18_34","child_hhmber_35_64","child_hhmber_over64")
  
  for (zz in 1:3) {
    lt_age <- lt_age_range[[zz]]
    # in age range 
    m2 <- data %>% filter(is.na(hh_role) & age %in% lt_age) %>% nrow()
    n <- min(m1[1,zz], m2)
    
    if( n>0 ){
      for (j in 1:n) {
        oo = 0
        n_loop = 1000
        while (oo==0 & n_loop > 0) {
          # Household info (row)
          a <- df %>% filter(type %in% lt_hhtype_temp) %>% slice_sample(n=1)
          # Child info (row)
          b <- data %>% filter(is.na(hh_role) & age %in% lt_age) %>% slice_sample(n=1)
          # Check if in range 
          oo <- as.numeric(as.character(cut(b$age, breaks = c(-999, a$child_lo_b -1, a$child_hi_b, 999), label = labls)))
          n_loop = n_loop-1
        }
        
        if(n_loop <=0){break}
        # SAVE result to dataframe 
        data[data$ind_id== b$ind_id,c("hh_id","hh_role")] <- c(a$hh_id, "child_over17")
        # print(paste(names(m1)[[zz]],"child_over17",j, sep = "  "))
      }
    }
  }
  
  foo[[idx]] <- data
  # print(paste(i,"done",sep = "----------------"))
}

# save result
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step7_populate_ind_child_over17.csv")
dfs_ind <- foo   # individual 
dfs_hh <- tt
```

### Assign child <18
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# Populate by census block group
lt_hhtype_temp <- c("h1_mc_udr18","h5_mns_udr18","h10_fns_udr18")
labls <- c(0,1,0)

# Join individual info to household 
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # Household size 
  a <- table(data$hh_id) %>% as.data.frame() %>% setnames(old="Freq", new="hh_size")
  df <- left_join(df, a, by=c("hh_id" ="Var1"))
  # Household with child over 17
  a <- data %>% filter(hh_role == "child_over17")
  b <- table(a$hh_id) %>% as.data.frame() %>% setnames(old="Freq", new="child_over17_count")
  
  if(nrow(a) > 0){
    df <- left_join(df, b, by=c("hh_id" ="Var1"))
    a <- a[order(a$age, decreasing = T),]
    a <- a[!duplicated(a$hh_id),]
    
    a$grand_child_lo_b <- 999
    a$grand_child_hi_b <- -1
    
    a[a$gender =="Female",]$grand_child_lo_b <- a[a$gender =="Female",]$age - 40
    a[a$gender =="Female",]$grand_child_hi_b <- a[a$gender =="Female",]$age - 18
    a[a$gender =="Male",]$grand_child_lo_b <- a[a$gender =="Male",]$age - 49
    a[a$gender =="Male",]$grand_child_hi_b <- a[a$gender =="Male",]$age - 18
    
    a <- setnames(a, old="age", new="child_over17_age")
    df <- left_join(df, a[c("hh_id","child_over17_age","grand_child_lo_b","grand_child_hi_b")], by=c("hh_id" = "hh_id"))
    
    df[is.na(df)] <- 0
  } else{
    df$child_over17_count <- 0
  }

  tt[[idx]] <- df
  # print(paste(i,"done",sep = "----------------"))
}

# Assign child udr 18

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # -------------
  hh_max_size = 6
  lt_temp <- df %>% filter( type %in% lt_hhtype_temp & hh_size <= hh_max_size) %>% nrow()
  n <- data %>% filter(is.na(hh_role) & age < 18) %>% nrow()
  if(n>0 & lt_temp>0){
    for (j in 1:n) {
      oo=0   # while loop stop condition oo=1
      o1=0
      o2=0
      n_loop = 1000
      while(oo==0 & n_loop > 0){
        # Child 
        a <- data %>% filter(is.na(hh_role) & age < 18) %>% slice_sample(n=1)
        # Household
        b <- df %>% filter( type %in% lt_hhtype_temp & hh_size <= hh_max_size) %>% slice_sample(n=1)
        o1 <- as.numeric(as.character(cut(a$age, breaks = c(-999, b$child_lo_b -1, b$child_hi_b, 999), label = labls)))
        
        if(b$child_over17_count > 0){
          o2 <- as.numeric(as.character(cut(a$age, breaks = c(-999, b$grand_child_lo_b -1, b$grand_child_hi_b, 999), label = labls)))
        }
        
        oo=o1+o2
        n_loop = n_loop-1
      }
      
      if(n_loop <=0){break}
      # Save Result
      hh_idx <- b$hh_id
      df[df$hh_id == hh_idx,]$hh_size <- df[df$hh_id == hh_idx,]$hh_size + 1
      data[data$ind_id == a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"child_or_grand_udr18")
      print(paste("child_udr18   ",j, sep = "  "))
    }
  }

  foo[[idx]] <- data
  tt[[idx]] <- df
  print(paste(i,"done",sep = "----------------"))
}

# # # save result
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step8_populate_ind_udr18_rest.csv")
dfs_ind <- foo
dfs_hh <- tt  # household
```

## Fill Type 6&11 Other relatives
```{r}
# set random seed 
set.seed(rm_seed)

# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# Add other relatives to hh type 6&11: h6_mns_NO18_others, h11_fns_NO18_others   
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  
  # ----ONLY for check---------
  # a <- table(data$hh_id) %>% as.data.frame()
  # df <- left_join(df, a, by=c("hh_id"="Var1"))
  # df$z <- df$Freq - df$hh_size
  # print(sum(df$hh_size, na.rm = T) - nrow(data[!is.na(data$hh_id) & data$hh_role!="in_gq",]))
  
  # ---- Houshold with other relatives but only have one member --------- 
  lt_temp <- df %>% filter(type %in% c("h11_fns_NO18_others", "h6_mns_NO18_others") & hh_size ==1) %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >= 18) %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]
      # randomly select an individual 
      a <- data %>% filter(is.na(hh_role) & age >= 18) %>% slice_sample(n=1)
      data[data$ind_id==a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"other_relatives")
    }
    
    
  }
  
  # ----- Non family household: Male -------
  lt_temp <- df %>% filter(type == "h7_mns_nonfamily") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >= 18 & gender == "Male") %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]
      a <- data %>% filter(is.na(hh_role) & age >= 18 & gender == "Male") %>% slice_sample(n=1)
      data[data$ind_id==a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"mns_non_family")
      
    }
  }
  
  # ----- Non family household: Female -------
  lt_temp <- df %>% filter(type == "h12_fns_nonfamily") %>% pull(hh_id)
  n <- data %>% filter(is.na(hh_role) & age >= 18 & gender == "Female") %>% nrow()
  n <- min(n, length(lt_temp))
  if(n>0){
    for (j in 1:n) {
      hh_idx <- lt_temp[[j]]
      a <- data %>% filter(is.na(hh_role) & age >= 18 & gender == "Female") %>% slice_sample(n=1)
      data[data$ind_id==a$ind_id, c("hh_id","hh_role")] <- c(hh_idx,"fns_non_family")
      
    }
  }
  
  # Randomly assign 1 non-family members
  m <- data %>% filter(hh_role %in% c("fns_non_family","mns_non_family")) %>% nrow()
  n <- data %>% filter(is.na(hh_role)) %>% nrow()
  n <- min(m, n)
  
  if(n >0){
    lt_temp <- data %>% filter(hh_role %in% c("fns_non_family","mns_non_family")) %>% slice_sample(n=n) %>% pull(hh_id)
    a <- data %>% filter(is.na(hh_role)) %>% slice_sample(n=n) %>% pull(ind_id)
    data[data$ind_id %in% a,]$hh_id <- lt_temp
    data[data$ind_id %in% a,]$hh_role <- "non_family_member"
    
  }
  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}

# Randomly assign the rest individuals 

for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # ----- rest members --------
  a <- table(data$hh_id) %>% as.data.frame()
  df <- left_join(df, a, by=c("hh_id"="Var1"))
  df$hh_size <- df$Freq
  
  df[is.na(df)] <- 0
  df <- df[df$hh_size > 0,]
  
  # unassigned individuals 
  n <- data %>% filter(is.na(hh_role)) %>% nrow()
  if(n > 0){
    for (j in 1:n) {
      a <- data %>% filter(is.na(hh_role)) %>% slice_sample(n=1) %>% pull(ind_id)
      lt_temp <- df %>% filter(!type %in% c("h3_mns_alone_less65","h4_mns_alone_over64",
                                            "h8_fns_alone_less65","h9_fns_alone_over64") & 
                                 hh_size < 7) %>% slice_sample(n=1) %>% pull(hh_id)
      data[data$ind_id == a,c("hh_id","hh_role")] <- c(lt_temp, "other_members")
      df[df$hh_id == lt_temp,]$hh_size <- df[df$hh_id == lt_temp,]$hh_size + 1
    }
  }
  
  foo[[idx]] <- data
  print(paste(i,"done",sep = "----------------"))
}




# # # save result
# data <- export_imp_dfs(foo)
# write.csv(data,"data/01_synthetic_population_dataset/process_saved_dataset/XX_step9_populate_ind_type6_7_11_12_other_relatives_non_family.csv")

dfs_ind <- foo
dfs_hh <- tt  # household
```

## Update household size 
```{r}
# Table
foo <- dfs_ind   # individual 
tt <- dfs_hh  # household 

# 1 -- For LOOP
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # ------------- update household size + add geoid
  a <- table(data$hh_id) %>% as.data.frame()
  df <- left_join(df, a, by=c("hh_id"="Var1"))
  df$hh_size <- df$Freq
  df$GEOID <- idx
  df <- df[c("GEOID","hh_id","type","hh_size")]
  df[is.na(df)] <- 0
  df <- df[df$hh_size > 0,]
  
  # ------------- individual df + add geoid
  data$GEOID <- idx
  
  foo[[idx]] <- data
  tt[[idx]] <- df
  print(paste(i,"done",sep = "----------------"))
}

finl_dfs_ind <- foo
finl_dfs_hh <- tt
```

# SAVE & EXPORT 
```{r}
# Individual 
a <- export_imp_dfs(finl_dfs_ind)
# Household 
b <- export_imp_dfs(finl_dfs_hh)

write.csv(a,"data/01_synthetic_population_dataset/process_saved_dataset/XX_final_synthesized_individual.csv", row.names = F)
write.csv(b,"data/01_synthetic_population_dataset/process_saved_dataset/XX_final_synthesized_household.csv", row.names = F)
```





# Validation
## V - Household
```{r}
df_chisq <- data.frame(matrix(ncol = 2)) %>% setnames(new = c("GEOID","chisq_hh"))
foo <- finl_dfs_ind
tt <- finl_dfs_hh

# 1 -- For LOOP
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  # Individual Data
  data <- foo[[idx]]
  # Household Data 
  df <- tt[[idx]]
  # -------------
  # census data 
  a <- css_hh_type_fnl %>% filter(GEOID == idx) %>% select(c(22:33)) %>% 
    gather(key="type", value="count")
  b <- table(df$type) %>% as.data.frame() %>% setnames(old="Freq", new="simu_count")
  a <- left_join(a, b, by=c("type"="Var1"))
  a[is.na(a)] <- 0
  
  c <- sum((a$simu_count - a$count)**2 / a$simu_count, na.rm = T)
  df_chisq[1+nrow(df_chisq),c("GEOID","chisq_hh")] <- c(idx, c)
  
  print(paste(i,"done",sep = "----------------"))
}


data <- df_chisq[(df_chisq$chisq_hh!="Inf" & !is.na(df_chisq$chisq_hh)),]
data$chisq_hh <- as.numeric(data$chisq_hh)
ggplot(data, aes(x=chisq_hh)) + geom_density()

```


```{r}
df <- data.frame(matrix(ncol = 10)) %>% setnames(new = c("GEOID","s_h10_fns_udr18","s_h5_mns_udr18",
                                                        "s_h9_fns_alone_over64","s_h4_mns_alone_over64","s_h2_mc_NO18",
                                                        "s_h3_mns_alone_less65","s_h8_fns_alone_less65",
                                                        "s_h6_mns_NO18_others","s_h11_fns_NO18_others")) 

data <- foo
for (i in 1:length(lt_geoid_bg)) {
  idx <- lt_geoid_bg[[i]]
  df[i,] <- c(
    idx,
    data[[idx]] %>% filter(hh_role=="single_mom") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="single_dad") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="female_oldersingleton") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="male_oldersingleton") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="husband_NO_udr18") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="mns_alone_less65") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="fns_alone_less65") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="mns_no18_other_relative") %>% nrow(),
    data[[idx]] %>% filter(hh_role=="fns_no18_other_relative") %>% nrow()
  )
}

df[,2:ncol(df)] <- sapply(df[,2:ncol(df)], as.numeric) %>% as.data.frame()
df <- left_join(df, css_hh_type_fnl[,c(1,22:33)], by=c("GEOID"="GEOID"))

# CHEKC 
df$s_h10_fns_udr18 - df$h10_fns_udr18
df$s_h5_mns_udr18 - df$h5_mns_udr18
df$s_h9_fns_alone_over64 - df$h9_fns_alone_over64
df$s_h4_mns_alone_over64 - df$h4_mns_alone_over64
df$s_h2_mc_NO18 - df$h2_mc_NO18
df$s_h3_mns_alone_less65 - df$h3_mns_alone_less65
df$s_h8_fns_alone_less65 - df$h8_fns_alone_less65
df$s_h6_mns_NO18_others - df$h6_mns_NO18_others
df$s_h11_fns_NO18_others - df$h11_fns_NO18_others

```

