# Make segregation measures using year-specific tract definitions


# Bring in MSAs
msas = read.csv('raw_data/msas/list1_2020.csv')
msas$msa = str_pad(msas$CBSA.Code,width = 5,pad = '0')
msas$mettype = str_replace(msas$Metropolitan.Micropolitan.Statistical.Area,pattern = ' Statistical Area','')

msas$st = str_pad(msas$FIPS.State.Code,width = 2,pad = '0')
msas$cnty = str_pad(msas$FIPS.County.Code, width = 3, pad = '0')
msas = filter(msas, mettype != '', st != '72')
msas = select(msas, st,cnty,mettype, msa)

# Bring in Census data
# 1970
cen7 = read_csv(paste0('raw_data/nhgis_segregation/',nhgis_seg_name,'_ds98_1970_tract.csv') )
cen7 = cen7 %>% select(st = STATEA, cnty = COUNTYA, tract = TRACTA,pop_b = C0X002,pop_w = C0X001, pop_o = C0X003)
cen7$year = 1969

# 1980-2010
cen_ts = read_csv(paste0('raw_data/nhgis_segregation/',nhgis_seg_name,'_ts_nominal_tract.csv'))
cen_ts = mutate(cen_ts,
             pop_b_1979 = CV4AB1980 + CV4AG1980,
             pop_b_1989 = CV4AB1990 + CV4AG1990,
             pop_b_1999 = CV4AB2000 + CV4AG2000, 
             pop_o_1979 = CV4AC1980 + CV4AD1980 + CV4AF1980 + CV4AH1980 + CV4AI1980, 
             pop_o_1989 = CV4AC1990 + CV4AD1990 + CV4AF1990 + CV4AH1990 + CV4AI1990,
             pop_o_1999 = CV4AC2000 + CV4AD2000 + CV4AF2000 + CV4AH2000 + CV4AI2000 + CV4AE2000 + CV4AJ2000
             ) %>% select(st = STATEFP,cnty = COUNTYFP, tract = TRACTA,
                         pop_w_1979 = CV4AA1980,
                         pop_w_1989 = CV4AA1990,
                         pop_w_1999 = CV4AA2000,
                         all_of(c(contains('pop_b'), contains('pop_o')))
             )

cen_ts = pivot_longer(cen_ts,-c(st,cnty,tract))
cen_ts$year = str_extract(cen_ts$name,'[0-9]{4}')
cen_ts$race = substr(cen_ts$name,5,5)
cen_ts = cen_ts %>% pivot_wider(id_cols = c(st,cnty,tract,year),names_from = race, names_prefix = 'pop_',values_from = value, values_fill = 0)
cen_ts = mutate(cen_ts,
             pop_b = replace_na(pop_b,replace = 0),
             pop_o = replace_na(pop_o,replace = 0),
             pop_w = replace_na(pop_w,replace = 0)
)

cen_ts$year = as.numeric(cen_ts$year)
cen_ts %>% count(year)



# 2005-2009
cen9 = read_csv(paste0('raw_data/nhgis_segregation/',nhgis_seg_name,'_ds195_20095_tract.csv'))
cen9$pop_b = cen9$RLIE004 + cen9$RLIE014
cen9$pop_o = cen9$RLIE001 - cen9$pop_b - cen9$RLIE003
cen9 = select(cen9, 
              st = STATEA,cnty = COUNTYA,tract = TRACTA,
              pop_w = RLIE003,pop_b,pop_o)
cen9$year = 2007

# 2010 - 2014
cen14 = read_csv(paste0('raw_data/nhgis_segregation/',nhgis_seg_name,'_ds206_20145_tract.csv'))
cen14 = mutate(cen14, pop_b = ABBAE004 + ABBAE014)
cen14$pop_o = cen14$ABBAE001 - cen14$pop_b - cen14$ABBAE003

cen14 = select(cen14, 
              st = STATEA,cnty = COUNTYA,tract = TRACTA,
              pop_w = ABBAE003,pop_b,pop_o)
cen14$year = 2012

# 2015-2019
cen19 = read_csv(paste0('raw_data/nhgis_segregation/',nhgis_seg_name,'_ds244_20195_tract.csv'))
cen19 = mutate(cen19, pop_b = ALUKE004 + ALUKE014)
cen19$pop_o = cen19$ALUKE001 - cen19$pop_b - cen19$ALUKE003

cen19 = select(cen19, 
               st = STATEA,cnty = COUNTYA,tract = TRACTA,
               pop_w = ALUKE003,pop_b,pop_o)

cen19$year = 2017

cen = bind_rows(cen7,cen_ts,cen9,cen14,cen19)
stopifnot(!is.na(cen$pop_b))
stopifnot(!is.na(cen$pop_w))
stopifnot(!is.na(cen$pop_o))

cen %>% group_by(year)%>% summarise(
  b = sum(is.na(pop_b)),
  w = sum(is.na(pop_w))
)

cen %>% group_by(year)%>% summarise(
  b = sum(pop_b),
  w = sum(pop_w),
  o = sum(pop_o)
) # Definitely are missing some in the 1970 and 1980 Census, these are presumably parts of the country not tracted at that point
# Total pop was 22 million black in 1970, 26.5 million in 1980, vs 17 and 24 here
# Total was 178 million white in 1970 vs 125 million here, 188 vs 154 in 1980 (those are for all white include hisp though)

# Add on MSA
cen = left_join(cen,msas, by = c('st','cnty'))

cen = cen %>% mutate(msafull = case_when(mettype == 'Metropolitan' ~ msa,
                                          mettype == 'Micropolitan' ~ paste('nm',st,sep = '_'),
                                          is.na(mettype) ~ paste('nm',st,sep = '_'))) %>% select(-msa)

cen$msa = cen$msafull

  
cenlong = cen %>% select(year,msa,st,cnty,tract,pop_b,pop_w) %>% pivot_longer(
  cols = c(pop_w, pop_b),
  names_to = "race",
  values_to = "population"
) %>% mutate(fips = paste0(st,cnty,tract) ) 


cenlongtot = cen %>% select(year,msa,st,cnty,tract,pop_b,pop_w,pop_o) %>% pivot_longer(
  cols = c(pop_w, pop_b,pop_o),
  names_to = "race",
  values_to = "population"
) %>% mutate(fips = paste0(st,cnty,tract) ) 



dissimilarity_results = cenlong %>%
  # Group by CBSA
  group_by(msa,year) %>%
  # Calculate multiple segregation indices for each CBSA
  summarise(
    # Dissimilarity Index (evenness)
    dissimilarity_index = if(n_distinct(fips) <= 1) {
      0  # Return 0 if only one tract (no segregation possible)
    } else {
      dissimilarity(
        data = cur_data(),
        group = "race",
        unit = "fips",
        weight = "population"
      )$est  # Extract the estimate from the dissimilarity object
    },

    
    
    exposure_b_w = if(n_distinct(fips) <= 1 ) {
      #   # For single tract, isolation equals the proportion Black in that tract
      sum(population[race == "pop_w"], na.rm = TRUE) / sum(population, na.rm = TRUE)
    } else if(sum(population * (race == 'pop_b')) == 0) {
      0
    } else {
      exposure(data = cur_data(), 
               group = "race",
               unit = "fips",
               weight = "population") %>% 
        filter(of == 'pop_b',to == 'pop_w') %>% select(exposure) %>% unlist()
    },
    
    theil = if(n_distinct(fips) <= 1 ) {
      0
      } else if(sum(population * (race == 'pop_b')) == 0) {
      0
    } else {
      mutual_total(data = cur_data(), 
               group = "race",
               unit = "fips",
               weight = "population") %>% 
        filter(stat == 'H') %>% select(est) %>% unlist()
    },
    
    total_population = sum(population, na.rm = TRUE),
    total_black_pop = sum(population[race == "pop_b"], na.rm = TRUE),
    total_white_pop = sum(population[race == "pop_w"], na.rm = TRUE),
    num_tracts = n_distinct(fips),
    ) %>% ungroup()


iso_results = cenlongtot %>%
  # Group by CBSA
  group_by(msa,year) %>%
  # Calculate multiple segregation indices for each CBSA
  summarise(
iso_black =  if(n_distinct(fips) <= 1 ) {
  #   # For single tract, isolation equals the proportion Black in that tract
  sum(population[race == "pop_b"], na.rm = TRUE) / sum(population, na.rm = TRUE)
} else if(sum(population * (race == 'pop_b')) == 0) {
  0
} else {
  
  isolation(
    data = cur_data(),
    group = "race",
    unit = "fips",
    weight = "population"
  ) %>% filter(race == 'pop_b') %>% select(isolation) %>% unlist()
}
)

seg_results = full_join(dissimilarity_results, iso_results, by = c('msa','year'))
stopifnot(nrow(seg_results) == nrow(dissimilarity_results))



write_csv(seg_results,'intermediate_data/predictors/msa_segregation_nhgis.csv')

