process_auto_counts <- function(auto_counts) {
    
    auto_counts$interval_60 <- as.POSIXct(auto_counts$interval_60, origin="1970-01-01", tz="UTC")
    auto_counts$hour <- hour(auto_counts$interval_60)
    auto_counts$day <- as.Date(auto_counts$interval_60)
    auto_counts$weekday <- wday(auto_counts$interval_60, week_start = 1)
    auto_counts$is_weekend <- auto_counts$weekday %in% c(6, 7)
    return(auto_counts)
}

auto_counts <- process_auto_counts(auto_counts)
process_cabi_data <- function(cabi_data) {
    cabi_data$hour <- ifelse(is.na(cabi_data$hour), -1, as.integer(cabi_data$hour))
    cabi_data$dotw <- ifelse(is.na(cabi_data$dotw), -1, as.integer(cabi_data$dotw))
    cabi_data$is_weekend <- cabi_data$dotw %in% c(6, 7)
    return(cabi_data)
}
cabi_data <- process_cabi_data(cabi_sum)

plot_hourly_patterns <- function(auto_counts, cabi_data) {
    hourly_auto <- auto_counts %>% group_by(hour) %>% summarize(Flow_Count = mean(Flow_Count, na.rm=TRUE))
    hourly_cabi <- cabi_data %>% group_by(hour) %>% summarize(cabi_count = mean(cabi_count, na.rm=TRUE))
    
    p1 <- ggplot(hourly_auto, aes(x = hour, y = Flow_Count, fill = "Auto Counts")) + 
        geom_bar(stat="identity", color = "#C4AEE8") + 
        scale_fill_manual(values = c("Auto Counts" = "#562DAB")) + 
        ggtitle("Average Hourly Auto Counts") + 
        xlab("Hour of Day") + 
        ylab("Average Count") + 
        theme_minimal()

    p2 <- ggplot(hourly_cabi, aes(x = hour, y = cabi_count, fill = "CABI Usage")) + 
        geom_bar(stat="identity", color = "#C4AEE8") + 
        scale_fill_manual(values = c("CABI Usage" = "#6A5ACD")) + 
        ggtitle("Average Hourly CABI Usage") + 
        xlab("Hour of Day") + 
        ylab("Average Count") + 
        theme_minimal()

weekday_auto <- auto_counts %>% group_by(hour, is_weekend) %>% summarize(Flow_Count = mean(Flow_Count, na.rm=TRUE))

weekday_cabi <- cabi_data %>% group_by(hour, is_weekend) %>% summarize(cabi_count = mean(cabi_count, na.rm=TRUE))

    p3 <- ggplot(weekday_auto, aes(x = hour, y = Flow_Count, color = factor(is_weekend))) + 
        geom_line(size = 1.2) + 
        scale_color_manual(values = c("FALSE" = "#C4AEE8", "TRUE" = "#6A5ACD")) + 
        ggtitle("Weekday vs Weekend - Auto Counts") + 
        xlab("Hour of Day") + 
        ylab("Average Count") + 
        theme_minimal()

    p4 <- ggplot(weekday_cabi, aes(x = hour, y = cabi_count, color = factor(is_weekend))) +
        geom_line(size = 1.2) + 
        scale_color_manual(values = c("FALSE" = "#C4AEE8", "TRUE" = "#6A5ACD")) + 
        ggtitle("Weekday vs Weekend - CABI Usage") + 
        xlab("Hour of Day") + 
        ylab("Average Count") + 
        theme_minimal()
    return(list(p1, p2, p3, p4))
}

analyze_peak_hours <- function(auto_counts, cabi_data) {
    morning_peak <- 7:9
    evening_peak <- 16:18

    auto_peaks <- data.frame(
        morning = mean(auto_counts$Flow_Count[auto_counts$hour %in% morning_peak], na.rm=TRUE),
        evening = mean(auto_counts$Flow_Count[auto_counts$hour %in% evening_peak], na.rm=TRUE),
        off_peak = mean(auto_counts$Flow_Count[!(auto_counts$hour %in% c(morning_peak, evening_peak))], na.rm=TRUE))

    cabi_peaks <- data.frame(
        morning = mean(cabi_data$cabi_count[cabi_data$hour %in% morning_peak], na.rm=TRUE),
        evening = mean(cabi_data$cabi_count[cabi_data$hour %in% evening_peak], na.rm=TRUE),
        off_peak = mean(cabi_data$cabi_count[!(cabi_data$hour %in% c(morning_peak, evening_peak))], na.rm=TRUE)
    )
    return(list(Auto_Counts = auto_peaks, CABI_Usage = cabi_peaks))
}

plots <- plot_hourly_patterns(auto_counts, cabi_data)
peak_stats <- analyze_peak_hours(auto_counts, cabi_data)

print(plots[[2]])

print(plots[[4]])

seg_nongem <- st_drop_geometry(segments)
landuse <- seg_nongem %>%
  select(Commercial, Industrial, Institutional, Mixed, Recreational, Religious, Residential, Vacant, Other) %>%
  summarise_all(sum) %>%
  pivot_longer(cols = everything(), names_to = "Land_Use_Type", values_to = "Count")

ggplot(landuse, aes(x = reorder(Land_Use_Type, Count), y = Count, fill = Land_Use_Type)) +   
  geom_col(width = 0.7) +  # Adjust bar width (default is 0.9, reduce for thinner bars)
  coord_flip() +   
scale_fill_manual(values = c(
  "Vacant"        = "#562DAB",  # Very dark purple
  "Residential"   = "#6A3BB8",
  "Institutional" = "#7A50C7",
  "Commercial"    = "#8E68CC",
  "Religious"     = "#9C7BD8",
  "Mixed"         = "#A989E0",
  "Other"         = "#C4AEE8",
  "Industrial"    = "#D6C4F0",
  "Recreational"  = "#E3D7F5"   # Lightest
), guide = "none") +  
  labs(title = "Distribution of Bike Lanes Across Adjacent Land Use Types",
       x = "Land Use Type", y = "Count") +   
  theme_minimal() +
  theme(
    legend.position = "right",   # Keep legend on the right
    panel.grid.major = element_blank())

ggplot() +
  geom_sf(data = basemap, fill = "grey90", color = "white") +
  geom_sf(
    data = segments,
    aes(color = median_incomeE),
    alpha = 0.8) +
    scale_color_gradientn(
    colors = rev(pp),
    name = "Income",
    na.value = "grey80") +
  labs(title = "Median Income") +
  theme_minimal()

ggplot() +
  geom_sf(data = basemap, fill = "grey90", color = "white") +
  geom_sf(
    data = segments,
    aes(color = commute_bicycleE),
    alpha = 0.8) +
    scale_color_gradientn(
    colors = rev(pp),
    name = "commute_bicycle") +
  labs(title = "Bicycle Commuters") +
  theme_minimal()

ggplot() +
  geom_sf(data = basemap, fill = "grey90", color = "white") +
  geom_sf(
    data = segments,
    aes(color = pct_white),
    alpha = 0.8) +
    scale_color_gradientn(
    colors = rev(pp),
    name = "pct_white") +
  labs(title = "Percentage of White") +
  theme_minimal()

ggplot() +
  geom_sf(data = basemap, fill = "grey90", color = "white") +
  geom_sf(
    data = segments,
    aes(color = edu_bachelors_plusE),
    alpha = 0.8) +
    scale_color_gradientn(
    colors = rev(pp),
    name = " edu_bachelor_plus") +
  labs(title = "Education Level") +
  theme_minimal()

ggplot() +
  geom_sf(data = basemap, fill = "grey97", color = "white") +
  geom_sf(
    data = segments,
    aes(color = betweenness_max),
    alpha = 0.8) +
    scale_color_gradientn(
    colors = rev(pp),
    name = "Betweenness") +
  labs(title = "Bike Lane Network: Binned Betweenness Centrality") +
  theme_minimal()

segments_manual <- segments %>%
  semi_join(seg_manual, by = "seg_id") %>%     
  left_join(seg_manual, by = "seg_id")

ggplot() +
  geom_sf(data = basemap, fill = "grey97", color = "white") +
  geom_sf(
    data = segments_manual,
    aes(color = avg_count),
    alpha = 0.8) +
    scale_color_gradientn(
    colors = rev(pp),
    name = "Average Bike Counts") +
  labs(title = "Manual Bike Counts") +
  theme_minimal()

ggplot(segments_manual, aes(x = avg_count)) +
  geom_histogram(bins = 30, fill = "#7A50C7", color = "white") +
  theme_minimal() +
  labs(title = "Distribution of Average Manual Counts",
       x = "Average Count", y = "Number of Segments")

cabi <- read.csv('202408-capitalbikeshare-tripdata.csv')

cabi$duration <- difftime( as.POSIXct( cabi$ended_at),as.POSIXct( cabi$started_at), units = "hours") 
library(sf)
library(tidyr)
library(dplyr)
cabi <- cabi[!is.na(cabi$end_lat),]
start_points <- st_as_sf(cabi, coords = c("start_lng", "start_lat"), crs = 4326)
end_points <- st_as_sf(cabi, coords = c("end_lng", "end_lat"), crs = 4326)
cabi$distance_km <- as.numeric(st_distance(start_points, end_points, by_element = TRUE))/1000
cabi <- cabi %>%  cbind(cabi %>% unite("OD", start_lng,start_lat,end_lng,end_lat, sep = "_") %>% select(OD))
trip_sum= cabi %>% mutate(speed= distance_km/as.numeric(duration)) %>%  group_by(started_at,start_lng,start_lat,end_lng,end_lat, OD) %>% summarise(mean_speed=mean(speed), trip_count= n())

library(ggplot2)
library(sf)
# install.packages('r5r')
# Sys.setenv(JAVA_HOME = "C:/Program Files/Eclipse Adoptium/jdk-21.0.6.7-hotspot/bin/java.exe") 
library(r5r)
options(java.parameters = "-Xmx8G")
r5r_core <- setup_r5(data_path = 'C:/Users/Xian Lu Lee/OneDrive - PennO365/Datasets/Practicum DCDOT/Cabi')
# library(rJava)
dc_bound <- st_read('Washington_DC_Boundary.geojson')

  street_network <- street_network_to_sf(r5r_core) 
origins_sf <-  trip_sum %>% st_as_sf(coords = c("start_lng", "start_lat"), crs = 4326)
ggplot() +
  
  geom_sf(street_network$vertices$geometry, mapping=aes(), color = "gray", size = 0.3)+
  geom_sf(origins_sf,mapping=aes()) +
  geom_sf(dc_bound,mapping=aes())+
  theme_minimal() +
  labs(title = "R5R Street Network and Origins", color = "Origin ID")


origins_sf <- st_as_sf(trip_sum, coords = c("start_lng", "start_lat"), crs = 4326)
destinations_sf <- st_as_sf(trip_sum, coords = c("end_lng", "end_lat"), crs = 4326)

# Ensure CRS matches
dc_bound <- st_transform(dc_bound, crs = st_crs(origins_sf))

# Filter origins and destinations inside the shapefile
origins_within <- st_intersects(origins_sf,dc_bound, sparse = FALSE)
destinations_within <- st_intersects(destinations_sf, dc_bound, sparse = FALSE)

# Keep only rows where both origin and destination are inside the shapefile
trip_filtered <- trip_sum[origins_within & destinations_within, ]  %>% ungroup() %>% mutate(id=row_number())

#snap points
snapped_points <- find_snap(r5r_core, trip_filtered %>% rename('lon'='start_lng', 'lat'='start_lat') %>% select (id, lat, lon))
trip_filtered$start_lng <- snapped_points$lon
trip_filtered$start_lat <- snapped_points$lat

snapped_points <- find_snap(r5r_core, trip_filtered %>% rename('lon'='end_lng', 'lat'='end_lat') %>% select (id, lat, lon))
trip_filtered$end_lng <- snapped_points$lon
trip_filtered$end_lat <- snapped_points$lat
trip_filtered <- trip_filtered %>% filter(trip_count>10)


all_itineraries <- list()

# Loop through each trip
for (i in 1:nrow(trip_filtered)) {
  # Extract trip details
  origin <- data.frame(
    id = i,
    lon = trip_filtered$start_lng[i],
    lat = trip_sum$start_lat[i]
  )
  
  destination <- data.frame(
    id = i,
    lon = trip_filtered$end_lng[i],
    lat = trip_filtered$end_lat[i]
  )
  
  bike_speed <- trip_filtered$mean_speed[i]+3
  
  # Ensure bike speed is positive; otherwise, set a default value
  if (bike_speed <= 0) {
    bike_speed <- 12  # Default to 3 m/s (~10.8 km/h) if speed is 0
  }
  
  # Generate itinerary
  itinerary <- detailed_itineraries(
    r5r_core,
    origins = origin,
    destinations = destination,
    mode = "BICYCLE",
    bike_speed = bike_speed
  )
  
  # Store results
  all_itineraries[[i]] <- itinerary
  print(i)
}

bike_trips_df <- do.call(rbind, all_itineraries)

trip_sum$started_at <- as.POSIXct(cabi_sf$started_at)
library(lubridate)
trip_sum<- trip_sum %>% mutate(
              hour= hour(started_at),
              dotw= wday(started_at))
trip_sum <- trip_sum %>%group_by(hour,dotw,OD) %>% summarise(trip_count=sum(trip_count))
cabi_sf <- trip_sum %>% right_join(trip_filtered %>% select(OD,id))

cabi_sf <- bike_trips_df %>% select(from_id,geometry) %>% right_join(cabi_sf %>% mutate(id=as.character(id)), by= c("from_id"="id"))

library(tidygraph)
library(sfnetworks)

seg_network <- st_cast(seg_network, "LINESTRING")
seg_network <- as_sfnetwork(seg_network)

seg_network <- seg_network %>%
  activate("nodes") %>%
  mutate(
    degree = centrality_degree(),         # Number of direct connections
    betweenness = centrality_betweenness(), # Control over information flow
    closeness = centrality_closeness()      # Average distance to all nodes
  )

library(igraph)

# Convert to an igraph object for connectivity analysis
graph_obj <- as.igraph(seg_network)

# Check if the network is fully connected
network_connected <- is_connected(graph_obj)
print(network_connected)

# Identify connected components
components_info <- components(graph_obj)
print(components_info)

set.seed(717)

library(spatialsample)
# library(tidymodels)
library(recipes)
library(rsample)
library(parsnip)
library(workflows)
library(yardstick)
library(dplyr)
library(lubridate)
library(tidyverse)
library(tune)
library(sf)
theme_set(theme_bw())

"%!in%" <- Negate("%in%")
g <- glimpse

nn_function <- function(measureFrom,measureTo,k) {
  library(FNN)
  nn <-   
    FNN::get.knnx(measureTo, measureFrom, k)$nn.dist
  output <-
    as.data.frame(nn) %>%
    rownames_to_column(var = "thisPoint") %>%
    gather(points, point_distance, V1:ncol(.)) %>%
    arrange(as.numeric(thisPoint)) %>%
    group_by(thisPoint) %>%
    summarize(pointDistance = mean(point_distance)) %>%
    arrange(as.numeric(thisPoint)) %>% 
    dplyr::select(-thisPoint) %>%
    pull()
  
  return(output)  
}


get_week_of_month <- function(date) {
  week(date) - week(floor_date(date, "month")) + 1
}
#######train test split######
setwd('C:/Users/Xian Lu Lee/OneDrive - PennO365/Datasets/GeoCloud/DC_DOT')

cabi_sum <- read.csv('Data/cabi_sum.csv') 
segments <- st_read("C:/Users/Xian Lu Lee/OneDrive - PennO365/Datasets/GeoCloud/DC_DOT/Data/segments_network.geojson"
)
cabi_max <- cabi_sum %>% group_by(seg_id) %>% summarise(cabi_count=max(cabi_count))

# cabi_sum <- cabi_sum[cabi_sum$hour %in% c(7,8,9,10, 15, 16, 17, 18),] 
load("C:/Users/Xian Lu Lee/OneDrive - PennO365/Datasets/GeoCloud/DC_DOT/model/data_joined(1).RData")
manual_counts <- manual_234 %>% st_drop_geometry()
  # read.csv("Data/manual_21_to_24.csv")

# Function to standardize time format: Add leading zero & ensure ":00" is present
extract_hour_ampm <- function(time_str) {
  time_str <- str_trim(time_str)  # Remove extra spaces
  
  # Extract AM/PM
  am_pm <- ifelse(grepl("AM", time_str), "AM", "PM")
  
  # Extract numeric hour
  hour <- as.numeric(str_extract(time_str, "^[0-9]+"))
  
  return(data.frame(Hour = hour, AM_PM = am_pm, stringsAsFactors = FALSE))
}


manual <- segments %>%right_join(manual_counts %>% select(-SUBBLOCKKEY,-SUBBLOCKID) %>% ungroup(), by='seg_id' )%>%  filter(!is.na(Total))%>% 
  mutate(
    Start_Hour_AMPM = lapply(Start_Time, extract_hour_ampm),
    End_Hour_AMPM = lapply(End_Time, extract_hour_ampm)
  ) %>%
  tidyr::unnest_wider(Start_Hour_AMPM, names_sep = "_") %>%
  tidyr::unnest_wider(End_Hour_AMPM, names_sep = "_") %>%
  rename(Start_Hour = Start_Hour_AMPM_Hour, Start_AMPM = Start_Hour_AMPM_AM_PM,
         End_Hour = End_Hour_AMPM_Hour, End_AMPM = End_Hour_AMPM_AM_PM)%>%
  mutate(
    Hours_Duration = ifelse(End_AMPM == Start_AMPM, End_Hour - Start_Hour,
                            (12 - Start_Hour) + End_Hour)  # Adjust if AM/PM changes
  ) %>% group_by(seg_id, Start_AMPM) %>% mutate(Total=sum(Total))%>% ##remove direction and get totals
  ##
  mutate(avg_count= Total/Hours_Duration) %>% 
  ungroup() 
# centrality <- read.csv("Data/network.csv")


main <- manual[ ,c(1,2:4, 10:18, 22:63, 95, 126,137, 141,146,148, 164, 168)] 
main %>% write.csv("manual_full.csv")
main <- main[!is.na(main$seg_id), ]
main <- main %>% mutate(Ward=as.factor(Ward))
main <- main %>% left_join(cabi_max)
# main <- main %>% left_join(centrality)
# main <- main %>% st_drop_geometry() %>% select(-geometry)
# main <- main %>% mutate(obs_weight = if_else(avg_count > 100, 10, if_else(avg_count > 25, 5, 1)))


# View(main[is.na(main$SUBBLOCKKEY), ])

### Initial Split for Training and Test
data_split <- initial_split(main , strata = "avg_count", prop = 0.80)
train_set <- training(data_split)
test_set  <- testing(data_split)

# This generates synthetic examples focused on rare/extreme values
train_set <- train_set %>%
  mutate(count_bin = cut(avg_count, breaks = 5)) 
train_set %>% ggplot(aes(x=count_bin))+geom_histogram(stat = 'count')
rare_bins <- train_set %>%
  group_by(count_bin) %>%
  summarise(count = n()) %>%
  filter(count < 50)  # 
# Select high-count rows
high_count_rows <- train_set %>% filter(avg_count > 100)

# Create synthetic rows by adding random noise to the high-count values
synthetic_rows <- high_count_rows %>%
  mutate(avg_count = avg_count + rnorm(n(), mean = 0, sd = 5)) %>% rbind(high_count_rows %>%
                                                                       mutate(avg_count = avg_count + rnorm(n(), mean = 0, sd = 5))) %>% 
   rbind(high_count_rows %>%
              mutate(avg_count = avg_count + rnorm(n(), mean = 0, sd = 5))) %>% 
  rbind(high_count_rows %>%
                                                                       mutate(avg_count = avg_count + rnorm(n(), mean = 0, sd = 5))) # Adjust sd to control noise level

# Combine original data with synthetic rows
train_set_synthetic <- bind_rows(train_set, synthetic_rows)
train_set_synthetic %>% ggplot(aes(x=count_bin))+geom_histogram(stat = 'count')
# Replace with newer spatialsample package approach
library(spatialsample)

# Spatial CV using newer approach

## LOGOCV on Neighborhood with group_vfold_cv()
cv_splits_geo <- group_vfold_cv(train_set,  
                                group = "GEOID")



# Get the first fold's split object
fold_1 <- cv_splits_geo$splits[[1]]

# Extract the analysis (training) set from the first fold
analysis_set_1 <- analysis(fold_1)

# Extract the assessment (validation) set from the first fold
assessment_set_1 <- assessment(fold_1)

### Create Recipes
test_set%>% select(where(~ is.factor(.) || is.character(.))) %>% colnames()
# Feature Creation
model_rec <- recipe(avg_count ~ ., data = train_set_synthetic
  %>% select(-seg_id,-count_bin,-BLOCKKEY, -has_parking,
  -has_buffer, -has_doubleyellow, 
  -GEOID #geoid,,
  # -cabi_count,
  # -Ward
  # -TOTALBIKELANES.x, -TOTALBIKELANEWIDTH.x,
  # -TRAIL_CLASS, -CENTERLINE,- LIGHTING
  ))%>% 
  update_role(Ward, new_role = "Neighborhood") %>%
  # update_role(obs_weight, new_role = "case_weight") %>% 
  step_other(Ward, threshold = 0.005) %>%
  step_mutate_at(all_predictors()&all_numeric(), fn = ~ ifelse(is.na(.), 1e-6, .)) %>%
  step_novel(all_nominal()) %>% 
  step_unknown(all_nominal()) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>% 
# step_log(avg_count, offset = 1, skip = TRUE) %>%
# step_nzv(all_predictors()&all_numeric(),-avg_count,freq_cut = 95/5, unique_cut = 10) %>%
  step_normalize(where(is.double),-avg_count)


# # step_zv(all_predictors()) %>%
# step_center(where(is.double) )%>%
#   step_scale(where(is.double) )

# step_ns(y,x, options = list(df = 4),keep_original_cols=T)

glimpse(model_rec %>% prep() %>% juice())

# transform new data
model_rec %>% prep() %>%  bake(new_data = main %>%select(-seg_id,-BLOCKKEY,-has_parking, -has_buffer, -has_doubleyellow,
                                                         
                                                         # -TOTALBIKELANES.x, -TOTALBIKELANEWIDTH.x,
                                                         # -TRAIL_CLASS, -CENTERLINE,- LIGHTING
                                                         ))

cor_matrix <- cor(main %>% select(where(is.numeric)),use='complete.obs')
cor_matrix['avg_count',]
########################################model testing###################################


## Model specifications
lm_plan <- 
  linear_reg() %>% 
  set_engine("lm", case.weights = TRUE)

glmnet_plan <-
  linear_reg() %>%
  set_args(penalty  = tune()) %>%
  set_args(mixture  = tune()) %>%
  set_engine("glmnet",case.weights = TRUE)%>%
  set_mode("regression")

rf_plan <- rand_forest() %>%
  set_args(mtry  = tune()) %>%
  set_args(min_n = tune()) %>%  # ✅ Now only tuning min_n
  set_args(n_trees= 2000) %>% 
  set_engine("ranger", importance = "permutation",case.weights = TRUE) %>%  
  set_mode("regression")


XGB_plan <- boost_tree() %>%
  set_args(trees = tune()) %>%  
  set_args(learn_rate = tune()) %>%  # Learning rate (shrinkage)
  set_args(tree_depth = tune()) %>%  # Maximum tree depth
  set_args(min_n = tune()) %>%       # Minimum node size
  # set_args(colsample_bytree = tune()) %>% # Feature subsampling
  set_args(early_stopping_rounds = 10) %>%       # Early stopping
  set_engine("xgboost",validation = 0.2, case.weights = TRUE,
             # objective="count:poisson",
             lambda=0.1, alpha=0.1) %>% 
  set_mode("regression")


# Hyperparameter grid for glmnet (penalization)
glmnet_grid <- expand.grid(
  penalty = 10^seq(-4, 2, length.out = 50),  # Wider range of penalties
  mixture = seq(0, 1, by = 0.25)
)

# Hyperparameter grid for random forest (tuning min_n only)
rf_grid <- expand.grid(
  mtry = c( 13, 30, 45, 80, 110), 
  min_n = c(5, 50, 100, 500)  # ✅ Only tuning min_n
)

# Hyperparameter grid for XGBoost (tuning n_estimators with early stopping)
xgb_grid <- expand.grid(
  trees = seq(10, 800, by = 50),  
  learn_rate = c(0.01, 0.1, 0.3),  
  tree_depth = c( 5, 6, 7),  
  min_n = c( 5, 50, 80) 
  # colsample_bytree = c(0.5, 0.7, 1)
)
#ntrees for xgboost does not need to be as big

# create workflow
lm_wf <-
  workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(lm_plan)
glmnet_wf <-
  workflow() %>%
  add_recipe(model_rec) %>%
  add_model(glmnet_plan)
rf_wf <-
  workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(rf_plan)
xgb_wf <-
  workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(XGB_plan)


# fit model to workflow and calculate metrics
control <- control_resamples(save_pred = TRUE, verbose = TRUE)
control <- control_grid(
  save_pred = TRUE,   # Ensures predictions are saved
  extract = function(x) x,  # Stores the full model object
  verbose = TRUE)
metrics <- metric_set(rmse, rsq, mape, smape)
lm_tuned <- lm_wf %>%
  tune::fit_resamples(.,
                      resamples = cv_splits_geo,
                      control   = control,
                      metrics   = metrics)

glmnet_tuned <- glmnet_wf %>%
  tune::tune_grid(.,
                  resamples = cv_splits_geo,
                  grid      = glmnet_grid,
                  control   = control,
                  metrics   = metrics)

rf_tuned <- rf_wf %>%
  tune::tune_grid(.,
                  resamples = cv_splits_geo,
                  grid      = rf_grid,
                  control   = control,
                  metrics   = metrics)

xgb_tuned <- xgb_wf %>%
  tune::tune_grid(.,
                  resamples = cv_splits_geo,
                  grid      = xgb_grid,
                  control   = control,
                  metrics   = metrics)


# ## metrics across grid
# collect_metrics(xgb_tuned)
# autoplot(xgb_tuned)
# collect_metrics(glmnet_tuned)



## 'Best' by some metric and margin
show_best(lm_tuned, metric = "rmse", n = 15)
show_best(glmnet_tuned, metric = "rmse", n = 15)
show_best(rf_tuned, metric = "rmse", n = 15)
show_best(xgb_tuned, metric = "rmse", n = 15)

lm_best_params     <- select_best(lm_tuned, metric = "rmse"    )
glmnet_best_params <- select_best(glmnet_tuned, metric = "rmse")
rf_best_params     <- select_best(rf_tuned, metric = "rmse"    )
xgb_best_params    <- select_best(xgb_tuned, metric = "rmse"   )


## Final workflow
lm_best_wf     <- finalize_workflow(lm_wf, lm_best_params)
glmnet_best_wf <- finalize_workflow(glmnet_wf, glmnet_best_params)
rf_best_wf     <- finalize_workflow(rf_wf, rf_best_params)
xgb_best_wf    <- finalize_workflow(xgb_wf, xgb_best_params)


# last_fit() emulates the process where, after determining the best model, the final fit on the entire training set is needed and is then evaluated on the test set.
lm_val_fit_geo <- lm_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)

glmnet_val_fit_geo <- glmnet_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)

rf_val_fit_geo <- rf_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)

xgb_val_fit_geo <- xgb_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)

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

# Pull best hyperparam preds from out-of-fold predictions
lm_best_OOF_preds <- collect_predictions(lm_tuned) 

glmnet_best_OOF_preds <- collect_predictions(glmnet_tuned) %>%
  filter(penalty  == glmnet_best_params$penalty[1] & mixture == glmnet_best_params$mixture[1])

rf_best_OOF_preds <- collect_predictions(rf_tuned) %>%
  filter(mtry  == rf_best_params$mtry[1] & min_n == rf_best_params$min_n[1])

xgb_best_OOF_preds <- collect_predictions(xgb_tuned) %>%
  filter(min_n == xgb_best_params$min_n[1]& .config == xgb_best_params$.config)

# collect validation set predictions from last_fit model
lm_val_pred_geo     <- collect_predictions(lm_val_fit_geo)
glmnet_val_pred_geo <- collect_predictions(glmnet_val_fit_geo)
rf_val_pred_geo     <- collect_predictions(rf_val_fit_geo)
xgb_val_pred_geo    <- collect_predictions(xgb_val_fit_geo)




# Aggregate OOF predictions (they do not overlap with Validation prediction set)
OOF_preds <- rbind(
                  data.frame(dplyr::select(lm_best_OOF_preds, .pred, avg_count), model = "lm") ,
                   data.frame(dplyr::select(glmnet_best_OOF_preds, .pred, avg_count), model = "glmnet"),
                   data.frame(dplyr::select(rf_best_OOF_preds, .pred, avg_count), model = "rf"),
                   data.frame(dplyr::select(xgb_best_OOF_preds, .pred, avg_count), model = "xgb")
                   # %>% mutate(.pred=expm1(.pred))
                  ) %>% 
  group_by(model) %>% 
  mutate(
         RMSE = yardstick::rmse_vec(avg_count, .pred),
         MAE  = yardstick::mae_vec(avg_count, .pred),
         MAPE = yardstick::mape_vec(avg_count, .pred)) %>% 
  ungroup() %>% 
  mutate(model = factor(model, levels=c(
    "lm",
    "glmnet",
    "rf",
    "xgb"
    )))



# average error for each model
ggplot(data = OOF_preds %>% 
         dplyr::select(model, RMSE) %>% 
         distinct() , 
       aes(x = model, y = RMSE, group = 1)) +
  geom_path(color = "red") +
  geom_label(aes(label = paste0(round(RMSE,1)))) +
  theme_bw()

# OOF predicted versus actual
ggplot(OOF_preds, aes(x = avg_count, y = .pred, group = model)) +
  geom_point(alpha = 0.3) +
  geom_abline(linetype = "dashed", color = "red") +
  geom_smooth(method = "lm", color = "blue") +
  # coord_equal() +
  facet_grid(~model) +
  ylim(0,300)+
  theme_bw()
full_fit    <-rf_best_wf %>% fit(train_set)

full_fit_lm     <- rf_best_wf %>% fit(main)


rf_pred_pm = predict(full_fit_lm, new_data = segments %>% left_join(manual_counts %>% select(seg_id,Ward, TOTALBIKELANES, TOTALBIKELANEWIDTH,
                                                                                             TRAIL_CLASS, CENTERLINE, LIGHTING), by='seg_id') %>% left_join(cabi_max) %>%
                       # rename(TOTALBIKELANES.x= TOTALBIKELANES, TOTALBIKELANEWIDTH.x= TOTALBIKELANEWIDTH) %>%

                       mutate(Start_AMPM= 'PM',
                                 Type.y=ifelse(is.na(road_class), "Trails", "Lanes"
                                               ),
                              Ward=as.factor(Ward)))
# rf_pred_pm$.pred<- rf_pred_pm$.pred %>% expm1()

rf_pred_am = predict(full_fit_lm, new_data = segments %>% left_join(manual_counts %>% select(seg_id,Ward, TOTALBIKELANES, TOTALBIKELANEWIDTH,
                                                                                                          TRAIL_CLASS, CENTERLINE, LIGHTING),by='seg_id') %>% left_join(cabi_max)%>%
                                    # rename(TOTALBIKELANES.x= TOTALBIKELANES, TOTALBIKELANEWIDTH.x= TOTALBIKELANEWIDTH) %>%

                                    mutate(Start_AMPM= 'AM',
                                           Type.y=ifelse(is.na(road_class), "Trails", "Lanes"
                                           ),
                                           Ward=as.factor(Ward)))
                                                        Type.y=ifelse(is.na(road_class), "Trails", "Lanes")))
# rf_pred_am$.pred<- rf_pred_am$.pred %>% expm1()


segments %>% left_join(manual_counts %>% select(seg_id,Ward, TOTALBIKELANES, TOTALBIKELANEWIDTH,
                                                        TRAIL_CLASS, CENTERLINE, LIGHTING),by='seg_id') %>% select(seg_id,GEOID,Ward) %>%
  cbind(rf_pred_am,rf_pred_pm) %>% distinct() %>%
  st_write('spatial_predictions6.csv', append = FALSE)

save(full_fit, test_set, rf_pred,OOF_preds, file = "model/Time_log_no_time_of_day_kde.RData")

load.fun <- function(x) { 
  x <- as.character(x) 
  if(isTRUE(x %in% .packages(all.available=TRUE))) { 
    eval(parse(text=paste("require(", x, ")", sep=""))) 
    print(paste(c(x, " : already installed; requiring"), collapse=''))
  } else { 
    #update.packages()
    print(paste(c(x, " : not installed; installing"), collapse=''))
    eval(parse(text=paste("install.packages('", x, "')", sep=""))) 
    print(paste(c(x, " : installed and requiring"), collapse=''))
    eval(parse(text=paste("require(", x, ")", sep=""))) 
  } 
} 

########### Required Packages ###########
packages = c(
  # Core ML
  "tidymodels", "tidyverse", 
  # Specific model packages
  "glmnet", "ranger", "xgboost", 
  # Newer spatial packages
  "sf", "stars", "terra", "spatialsample", "nngeo", "mapview",
  # Newer extensions
  "finetune", "themis", "applicable", "vetiver",
  # Existing packages
  "bayesplot", "lme4", "RcppEigen", "AmesHousing"
)

for(i in seq_along(packages)){
  packge <- as.character(packages[i])
  load.fun(packge)
}

set.seed(717)
theme_set(theme_bw())

"%!in%" <- Negate("%in%")
g <- glimpse

nn_function <- function(measureFrom,measureTo,k) {
  library(FNN)
  nn <-   
    FNN::get.knnx(measureTo, measureFrom, k)$nn.dist
  output <-
    as.data.frame(nn) %>%
    rownames_to_column(var = "thisPoint") %>%
    gather(points, point_distance, V1:ncol(.)) %>%
    arrange(as.numeric(thisPoint)) %>%
    group_by(thisPoint) %>%
    summarize(pointDistance = mean(point_distance)) %>%
    arrange(as.numeric(thisPoint)) %>% 
    dplyr::select(-thisPoint) %>%
    pull()
  
  return(output)  
}


get_week_of_month <- function(date) {
  week(date) - week(floor_date(date, "month")) + 1
}
#######train test split######
setwd('C:/Users/Xian Lu Lee/OneDrive - PennO365/Datasets/GeoCloud/DC_DOT/')
segments <-read.csv(
#   # "Data/Lanes_KumaHao.geojson"
  "/spatial_predictions6.csv"
  )
load("data_joined(1).RData")
# auto_counts <- st_read("Data/auto_2024_by_hr.csv")
auto_counts$interval_60 <- as.numeric(auto_counts$interval_60)
auto_counts$interval_60 <- as.POSIXct(auto_counts$interval_60, origin = "1970-01-01", tz = "UTC")
library(lubridate)
auto_counts <- auto_counts %>% 
  mutate(WkoMn = get_week_of_month(interval_60))

# full_weeks %>% select(Site_Name) %>% distinct()

cabi_sum <- read.csv('C:/Users/Xian Lu Lee/OneDrive - PennO365/Datasets/GeoCloud/DC_DOT/Data/cabi_sum.csv') %>% select(-X) %>% drop_na()
main <- segments  %>% left_join(auto_counts, by= "seg_id") %>% mutate(Flow_Count= as.numeric(Flow_Count)) %>%
  filter(!is.na(interval_60)) %>% 
  mutate(hour = as.numeric(hour(interval_60)),dotw = as.numeric(wday(interval_60, label=TRUE))) %>% st_drop_geometry()
main <- main %>% 
  left_join(cabi_sum, by= c('seg_id','hour', 'dotw')) 

main <- main  %>% 
  mutate(cabi_count = ifelse(is.na(cabi_count), 0, cabi_count))

# main <- main %>% drop_na()
main$hour= as.factor(main$hour)
main$.pred <- as.numeric( main$.pred)
main$.pred.1 <- as.numeric( main$.pred.1)
library(ggplot2)


main_summary <- main %>%
  count(Site_Name, dotw, WkoMn) %>% mutate(count=if_else(n>0, 1, 0))


# ggplot(weekdays, aes(y= Flow_Count, x= hour)) +geom_point()
# weekdays %>% group_by(hour) %>% summarise(mean=mean(Flow_Count)) %>% View()

ggplot(main_summary, aes(x = dotw, y = as.factor(WkoMn), fill = count)) +
  geom_tile(color = "white") +
  scale_fill_viridis_c(name = "Data Collected") +
  facet_wrap(~ Site_Name, scales = "free_y") +  # One heatmap per site
  labs(
    title = "Data Collection by Week of Month and Day of Week",
    x = "Day of the Week", y = "Week of the Month"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

main <- main %>% mutate(weekend = ifelse(dotw %in% c("Sun", "Sat"), "Weekend", "Weekday"),
                             time_of_day = case_when(hour(interval_60) < 7 | hour(interval_60) > 18 ~ "Overnight",
                                                     hour(interval_60) >= 7 & hour(interval_60) < 10 ~ "AM Rush",
                                                     hour(interval_60) >= 10 & hour(interval_60) < 15 ~ "Mid-Day",
                                                     hour(interval_60) >= 15 & hour(interval_60) <= 18 ~ "PM Rush"))

main <- main %>% group_by(Site_Name,interval_60) %>% mutate(Flow_Count= sum(Flow_Count)) %>% ungroup() %>% select(-Flow_Name,-Flow_Identifier) %>% distinct()###remove directionality because

bg <- st_read("Data/bg/bg2020.GEOJSON")
bg$GEOID <-bg$GEOCODE
manual_clean <- manual_234[complete.cases(st_coordinates(manual_234)), ] %>% st_transform(3857) 
# create kdem
kde_hp = sf.kde(manual_clean%>% select(Total), y=manual_clean$Total, bw=1000)
kde_df <- as.data.frame(kde_hp, xy = TRUE) %>%
  rename(density = 3)  
# map the result
ggplot(kde_df, aes(x = x, y = y, fill = density)) +
  geom_raster() +
  scale_fill_viridis_c(option = "plasma") +
  coord_equal() +
  theme_minimal() +
  labs(title = "Kernel Density Estimation", fill = "Density") 
# Convert sf points to SpatVector
pts_vect <- vect(main %>% st_as_sf())

# Extract KDE values
kde_vals <- extract(kde_hp, pts_vect)

# Combine with original sf object
main$manual_kde <- kde_vals[,2] 

# library(tidymodels)
library(rsample)
library(parsnip)
library(workflows)
library(yardstick)
main <- main %>% mutate(YearMonth = floor_date(interval_60, "month"))
main$GEOID <- as.factor(main$GEOID)
main$Ward<- as.factor(main$Ward)
main$dotw <- as.factor(main$dotw)
dat <- main %>% select(1:6, 8,9 ,12:18) %>% ungroup()%>%
  mutate(month=month(interval_60)) %>% 
    filter(month==8, weekend=="Weekday")
random_week_per_site <- dat %>%
  group_by(Site_Name) %>%
  distinct(WkoMn) %>%  # Get unique weeks per site
  filter(WkoMn %in% c(4:6)) %>% filter(Site_Name!='1st St NE') %>% 
  filter(Site_Name!= 'R Street NW')



# test_month <- sample(unique(dat$YearMonth), 1)
# 3. Create train and test sets (just key columns for matching)
test_set <- dat %>%
  inner_join(random_week_per_site, by = c("Site_Name", "WkoMn"))

train_set <- dat %>%anti_join(test_set)

# 4. Match and extract indices using row-wise logic
train_indices <- dat%>%
  mutate(row_id = row_number()) %>%
  semi_join(train_set, by = c("Site_Name", "interval_60")) %>%
  pull(row_id)

test_indices <- dat %>%
  mutate(row_id = row_number()) %>%
  semi_join(test_set, by = c("Site_Name", "interval_60")) %>%
  pull(row_id)

# 5. Create rsample split
library(rsample)
data_split <- make_splits(
  list(analysis = train_indices, assessment = test_indices),
  data = dat
)


## LOGOCV on Neighborhood with group_vfold_cv()
cv_splits_geo <- group_vfold_cv(train_set,  
                                group = "Ward")

## Accessing CV split data
glimpse(cv_splits_geo[1,1][[1]][[1]][[1]])

# Get the first fold's split object
fold_1 <- cv_splits_geo$splits[[1]]

# Extract the analysis (training) set from the first fold
analysis_set_1 <- analysis(fold_1)

# Extract the assessment (validation) set from the first fold
assessment_set_1 <- assessment(fold_1)

### Create Recipes
test_set%>% select(where(~ is.factor(.) || is.character(.))) %>% colnames()

# dat$Flow_Count <- log1p(dat$Flow_Count)
# train_set$Flow_Count <- log1p(train_set$Flow_Count)
# test_set$Flow_Count <- log1p(test_set$Flow_Count)
model_rec <- recipe(Flow_Count ~ ., data = train_set %>% 
                      select(-seg_id, -Site_Name,-GEOID, -interval_60, -month,
                             -WkoMn,-cabi_count)) %>%
  update_role(Ward, new_role = "Neighborhood") %>%
  step_other(Ward, threshold = 0.005) %>%
  step_normalize(where(is.double), -Flow_Count) %>% 
  step_unknown(all_nominal()) %>% 
  step_dummy(all_nominal(), one_hot = TRUE) %>% 
  # step_interact(~ hour:.pred) %>%   # Specify interaction term correctly
  step_log(Flow_Count, offset = 1, skip=T)
    
# step_nzv(all_numeric(),-Flow_Count,freq_cut = 95/5, unique_cut = 10) %>%

  # step_zv(all_predictors()) %>%
  # step_center(all_predictors(), -Flow_Count, -has_buffer, -has_parking, -has_doubleyellow) %>%
  # step_scale(all_predictors(), -Flow_Count, -has_buffer, -has_parking, -has_doubleyellow)
  
# step_ns(Latitude, Longitude, options = list(df = 4),keep_original_cols=TRUE)

glimpse(model_rec %>% prep() %>% juice())

# transform new data
model_rec %>% prep() %>%  bake(new_data = dat)



set.seed(42)
## Model specifications
lm_plan <- 
  linear_reg() %>% 
  set_engine("lm")

glmnet_plan <-
  linear_reg() %>%
  set_args(penalty  = tune()) %>%
  set_args(mixture  = tune()) %>%
  set_engine("glmnet")%>%
  set_mode("regression")


rf_plan <- rand_forest() %>%
  set_args(mtry  = tune()) %>%
  set_args(n_trees  = 2000) %>%
  set_args(min_n = tune()) %>%  # ✅ Now only tuning min_n
  set_engine("ranger", importance = "permutation") %>%  
  set_mode("regression")

XGB_plan <- boost_tree() %>%
  set_args(trees = tune()) %>%  
  set_args(learn_rate = tune()) %>%  # Learning rate (shrinkage)
  set_args(tree_depth = tune()) %>%  # Maximum tree depth
  set_args(min_n = tune()) %>%       # Minimum node size
  # set_args(colsample_bytree = tune()) %>% # Feature subsampling
  set_args(early_stopping_rounds = 10) %>%        # Early stopping
  set_engine("xgboost",validation = 0.2) %>% 
  set_mode("regression")

# poisson_mod <- 
#   poisson_reg() %>% 
#   set_engine("glm") %>% 
#   set_mode("regression")


# Hyperparameter grid for glmnet (penalization)
glmnet_grid <- expand.grid(
  penalty = 10^seq(-4, 2, length.out = 50),  # Wider range of penalties
  mixture = seq(0, 1, by = 0.25)
)

# Hyperparameter grid for random forest (tuning min_n only)
rf_grid <- expand.grid(
  mtry = c( 45, 80, 110), 
  min_n = c(5, 50, 100, 500)  # ✅ Only tuning min_n
)
# Hyperparameter grid for XGBoost (tuning n_estimators with early stopping)
xgb_grid <- expand.grid(
  trees = seq(50, 500, by = 50),  
  learn_rate = c(0.01, 0.1, 0.3),  
  tree_depth = c( 5, 7, 15,25),  
  min_n = c( 50, 200, 400) 
  # colsample_bytree = c(0.5, 0.7, 1)
)


library(poissonreg)
# create workflow
lm_wf <-
  workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(lm_plan)
# poisson_wf <- workflow() %>%  
  # add_recipe(model_rec) %>% 
  # add_model(poisson_mod)
  
glmnet_wf <-
  workflow() %>%
  add_recipe(model_rec) %>%
  add_model(glmnet_plan)
rf_wf <-
  workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(rf_plan)
xgb_wf <-
  workflow() %>% 
  add_recipe(model_rec) %>% 
  add_model(XGB_plan)





# fit model to workflow and calculate metrics
control <- control_resamples(save_pred = TRUE, verbose = TRUE)
metrics <- metric_set(rmse, rsq, mape, smape)
lm_tuned <- lm_wf %>%
  tune::fit_resamples(.,
                      resamples = cv_splits_geo,
                      control   = control,
                      metrics   = metrics)

glmnet_tuned <- glmnet_wf %>%
  tune::tune_grid(.,
                  resamples = cv_splits_geo,
                  grid      = glmnet_grid,
                  control   = control,
                  metrics   = metrics)

rf_tuned <- rf_wf %>%
  tune::tune_grid(.,
                  resamples = cv_splits_geo,
                  grid      = rf_grid,
                  control   = control,
                  metrics   = metrics)

xgb_tuned <- xgb_wf %>%
  tune::tune_grid(.,
                  resamples = cv_splits_geo,
                  grid      = xgb_grid,
                  control   = control,
                  metrics   = metrics)

## 'Best' by some metric and margin
show_best(lm_tuned, metric = "rmse", n = 15)
# show_best(glmnet_tuned, metric = "rmse", n = 15)
show_best(rf_tuned, metric = "rmse", n = 15)
show_best(xgb_tuned, metric = "rmse", n = 15)

lm_best_params     <- select_best(lm_tuned, metric = "rmse"    )
# glmnet_best_params <- select_best(glmnet_tuned, metric = "rmse")
rf_best_params     <- select_best(rf_tuned, metric = "rmse"    )
xgb_best_params    <- select_best(xgb_tuned, metric = "rmse"   )


## Final workflow
lm_best_wf     <- finalize_workflow(lm_wf, lm_best_params)
# glmnet_best_wf <- finalize_workflow(glmnet_wf, glmnet_best_params)
rf_best_wf     <- finalize_workflow(rf_wf, rf_best_params)
xgb_best_wf    <- finalize_workflow(xgb_wf, xgb_best_params)


# last_fit() emulates the process where, after determining the best model, the final fit on the entire training set is needed and is then evaluated on the test set.
lm_val_fit_geo <- lm_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)

# glmnet_val_fit_geo <- glmnet_best_wf %>% 
#   last_fit(split     = data_split,
#            control   = control,
#            metrics   = metrics)

rf_val_fit_geo <- rf_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)

xgb_val_fit_geo <- xgb_best_wf %>% 
  last_fit(split     = data_split,
           control   = control,
           metrics   = metrics)
# poisson_val_fit_geo <- poisson_wf %>% 
#   last_fit(split     = data_split,
#            control   = control,
#            metrics   = metrics)

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

# Pull best hyperparam preds from out-of-fold predictions
lm_best_OOF_preds <- collect_predictions(lm_tuned)

# glmnet_best_OOF_preds <- collect_predictions(glmnet_tuned) %>%
#   filter(penalty  == glmnet_best_params$penalty[1] & mixture == glmnet_best_params$mixture[1])

rf_best_OOF_preds <- collect_predictions(rf_tuned) %>%
  filter(mtry  == rf_best_params$mtry[1] & min_n == rf_best_params$min_n[1])

xgb_best_OOF_preds <- collect_predictions(xgb_tuned) %>%
  filter( min_n == xgb_best_params$min_n[1])

# collect validation set predictions from last_fit model
lm_val_pred_geo     <- collect_predictions(lm_val_fit_geo)
glmnet_val_pred_geo <- collect_predictions(glmnet_val_fit_geo)
rf_val_pred_geo     <- collect_predictions(rf_val_fit_geo)
xgb_val_pred_geo    <- collect_predictions(xgb_val_fit_geo)
# poisson_val_pred_geo    <- collect_predictions(poisson_val_fit_geo)



# Aggregate OOF predictions (they do not overlap with Validation prediction set)
OOF_preds <- rbind(
  data.frame(dplyr::select(lm_best_OOF_preds, .pred, Flow_Count)%>% mutate(.pred=expm1(.pred)), model = "lm") ,
   # data.frame(dplyr::select(poisson_val_pred_geo, .pred, Flow_Count), model = "poisson")
  # data.frame(dplyr::select(glmnet_best_OOF_preds, .pred, Flow_Count), model = "glmnet") ,
  data.frame(dplyr::select(rf_best_OOF_preds, .pred, Flow_Count)%>% mutate(.pred=expm1(.pred)), model = "rf"),
  data.frame(dplyr::select(xgb_best_OOF_preds, .pred, Flow_Count) %>% mutate(.pred=expm1(.pred)), model = "xgb")
  
) %>% 
  group_by(model) %>% 
  mutate(
    RMSE = yardstick::rmse_vec(Flow_Count, .pred),
    MAE  = yardstick::mae_vec(Flow_Count, .pred),
    MAPE = yardstick::mape_vec(Flow_Count, .pred)) %>% 
  ungroup() %>% 
  mutate(model = factor(model, levels=c(
    "lm",
    "glmnet",
    "rf",
    "xgb"
  )))
# average error for each model
ggplot(data = OOF_preds %>% 
         dplyr::select(model, RMSE) %>% 
         distinct() , 
       aes(x = model, y = RMSE, group = 1)) +
  geom_path(color = "red") +
  geom_label(aes(label = paste0(round(RMSE,1)))) +
  theme_bw()

# OOF predicted versus actual
ggplot(OOF_preds, aes(x = Flow_Count, y = .pred,
                      group = model
                      )) +
  geom_point(alpha = 0.3) +
  geom_abline(linetype = "dashed", color = "red") +
  geom_smooth(method = "lm", color = "blue") +
  # coord_equal() +
  facet_wrap(~model, nrow = 2) +
  theme_bw()

full_fit     <- rf_best_wf %>% workflows::fit(train_set)
rf_pred <- predict(full_fit, test_set %>% select(-Flow_Count))
rf_pred$.pred=expm1(rf_pred$.pred)
full_fit$fit$fit$fit