train <- py$train
test <- py$test

# Keep IDs for submission
ids <- test %>% pull(Id)
test <- test %>% select(-Id)

# Keep outcome variable
y <- train$SalePrice

print(glue("The training dataset has {ncol(train)} features and {nrow(train)} observations."))

map_chr(train, class) %>% table() %>% {glue("There are {.} features of type {names(.)}")}

glimpse(train)

train <- train %>% mutate_if(is.list, unlist)

train <- train %>% na_if("NaN")

test <- test %>% mutate_if(is.list, unlist)

test <- test %>% na_if("NaN")

any(map_lgl(train, is.list))

map_dbl(train, ~ sum(is.na(.x))) %>% 
  table() %>% 
  {glue("There {ifelse(. == 1, 'is', 'are')} {.} {ifelse(. == 1, 'feature', 'features')} with {names(.)} missing values")}

train %>% 
  summarise_all(~sum(is.na(.))) %>% 
  pivot_longer(cols = everything()) %>% 
  filter(value != 0) %>% 
  arrange(desc(value)) %>% 
  {glue("{.[['name']]} has {.[['value']]} missing observations")}

train %>% 
  summarise_all(~round(sum(is.na(.))/length(.), 4) * 100) %>% 
  pivot_longer(cols = everything()) %>% 
  filter(value != 0) %>% 
  arrange(desc(value)) %>% 
  {glue("{.[['name']]} has {.[['value']]}% missing observations")}

p1 <- train %>%
  ggplot(aes(x = SalePrice)) +
  geom_histogram(aes(x = SalePrice, stat(density)),
                 bins = 100,
                 fill = "cornflowerblue",
                 alpha = 0.7) +
  geom_density(color = "midnightblue") +
  scale_x_continuous(breaks= seq(0, 800000, by=100000), labels = scales::comma) +
  labs(x = "Sale Price", y = "", title = "Density of Sale Price") +
  theme(axis.text.x = element_text(size = 10, angle = 45, hjust = 1))

p1

p2 <- train %>%
  ggplot() +
  geom_histogram(aes(x = log(SalePrice), stat(density)),
                 bins = 100,
                 fill = "cornflowerblue",
                 alpha = 0.7) +
  geom_density(aes(x = log(SalePrice)), color = "midnightblue") +
  labs(x = "Natural log of Sale Price", y = "", title = "Density of Logged Sale Price")

p1 + p2

p1 <- train %>%
  ggplot(aes(x = GrLivArea, y = SalePrice)) +
  geom_point(color = "cornflowerblue", alpha = 0.75) +
  scale_y_continuous(breaks= seq(0, 800000, by=200000), labels = scales::comma) +
  labs(x = "Above ground living area", y = "Sale Price", title = "Sale Price \nby Above Ground Living Area")

p2 <- train %>%
  ggplot(aes(x = GrLivArea, y = log10(SalePrice))) +
  geom_point(color = "cornflowerblue", alpha = 0.75) +
  labs(x = "Above ground living area", y = "Log of Sale Price", title = "Logged Sale Price \nby Above Ground Living Area")

p1 + p2  

p1 <- train %>%
  ggplot(aes(x = GrLivArea, y = SalePrice)) +
  geom_point(aes(color = factor(OverallQual)), alpha = 0.75) +
  scale_color_viridis_d(name = "Overall Quality", 
                        breaks = c("2", "4", "6", "8", "10")) +
  scale_y_continuous(breaks= seq(0, 800000, by=200000), labels = scales::comma) +
  labs(x = "Above ground living area", y = "Sale Price", title = "Sale Price \nby Above Ground Living Area") +
  theme(legend.position = c(0.2, 0.8),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 6),
        legend.background = element_blank(),
        legend.key.size = unit(5, "mm"),
        #legend.box.background = element_rect(color = "black", fill = "transparent"),
        legend.margin = ggplot2::margin(0.1, 0.1, 0.1, 0.1, "cm"))

p2 <- 
  train %>%
  mutate(MSZoning = factor(MSZoning, levels = c("RL", "RM", "C (all)", "FV", "RH"))) %>% 
  ggplot(aes(x = GrLivArea, y = SalePrice)) +
  geom_point(aes(color = MSZoning), alpha = 0.75) +
  scale_color_viridis_d(name = "Zoning", end = 0.7) +
  scale_y_continuous(breaks= seq(0, 800000, by=200000), labels = scales::comma) +
  labs(x = "Above ground living area", y = "Sale Price", title = "Sale Price \nby Above Ground Living Area") +
  theme(legend.position = c(0.2, 0.8),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 6),
        legend.background = element_blank(),
        legend.key.size = unit(5, "mm"),
        #legend.box.background = element_rect(color = "black", fill = "transparent"),
        legend.margin = ggplot2::margin(0.1, 0.1, 0.1, 0.1, "cm"))

p1 + p2

train <-
  train %>%
  filter(
    !(GrLivArea > 4000 & SalePrice < 200000)
  )

all_dat <- bind_rows(train, test) %>% 
  select(-Id, -SalePrice)

none_vars <- c("BsmtCond",
               "BsmtQual",
               "BsmtExposure",
               "BsmtFinType1",
               "BsmtFinType2",
               "MasVnrType",
               "GarageType",
               "GarageFinish",
               "GarageQual",
               "GarageCond",
               "Alley",
               "Fence",
               "FireplaceQu",
               "MiscFeature",
               "PoolQC")

all_dat <- 
  all_dat %>% 
  mutate_at(vars(all_of(none_vars)), ~ replace_na(., "None"))

zero_vars <- c("MasVnrArea",
               "BsmtFinSF1",
               "BsmtFinSF2",
               "BsmtFullBath",
               "BsmtHalfBath",
               "BsmtUnfSF",
               "TotalBsmtSF",
               "GarageArea",
               "GarageCars",
               "GarageYrBlt")

all_dat <- 
  all_dat %>% 
  mutate_at(vars(all_of(zero_vars)), ~ replace_na(., 0))

random_vars <- c("Electrical", "Exterior1st", "Exterior2nd", "Functional", "KitchenQual", "MSZoning", "SaleType", "Utilities")

all_dat <- 
  all_dat %>% 
  group_by(Neighborhood, MSSubClass) %>% 
  mutate_at(vars(all_of(random_vars)), ~ ifelse(is.na(.), which.max(table(.)) %>% names(), .)) %>% 
  group_by(Neighborhood) %>% 
  mutate(LotFrontage = ifelse(is.na(LotFrontage), median(LotFrontage, na.rm = TRUE), LotFrontage)) %>% 
  ungroup()

glue("There are {sum(is.na(all_dat))} features with missing observations")

all_dat <-
  all_dat %>% 
  mutate(TotalArea = `1stFlrSF` + `2ndFlrSF` + TotalBsmtSF,
         AvgRmSF = GrLivArea / TotRmsAbvGrd,
         TotalBaths = FullBath + (HalfBath * 0.5) + BsmtFullBath + (BsmtHalfBath * 0.5),
         OutsideArea = OpenPorchSF + `3SsnPorch` + EnclosedPorch + ScreenPorch + WoodDeckSF
         ) %>% 
  mutate(TotalArea2 = TotalArea^2,
         OverallQual2 = OverallQual^2)

qual_vars = c("BsmtCond", "BsmtQual", "ExterCond", 
             "ExterQual", "FireplaceQu", "GarageCond", 
             "GarageQual", "HeatingQC", "KitchenQual", "PoolQC")
all_dat <-
  all_dat %>% 
  mutate_at(vars(all_of(qual_vars)), ~ factor(., 
                                              levels = c("None",
                                                         "Po",
                                                         "Fa",
                                                         "TA",
                                                         "Gd",
                                                         "Ex"), 
                                              ordered = TRUE)) %>% 
  mutate(Functional = factor(Functional, 
                             levels = c("Sev",
                                        "Maj1",
                                        "Maj2",
                                        "Mod",
                                        "Min1",
                                        "Min2",
                                        "Typ"),
                             ordered = TRUE),
         GarageFinish = factor(GarageFinish,
                               levels = c("None",
                                          "Unf",
                                          "RFn",
                                          "Fin"),
                               ordered = TRUE),
         BsmtFinType1 = factor(BsmtFinType1,
                               levels = c("None",
                                          "Unf",
                                          "LwQ",
                                          "Rec",
                                          "BLQ",
                                          "ALQ",
                                          "GLQ"),
                               ordered = TRUE),
         BsmtFinType2 = factor(BsmtFinType2,
                               levels = c("None",
                                          "Unf",
                                          "LwQ",
                                          "Rec",
                                          "BLQ",
                                          "ALQ",
                                          "GLQ"),
                               ordered = TRUE))

# To make these match the features in Python, I'm going to convert them to numeric
all_dat <-
  all_dat %>% 
  mutate_at(vars(all_of(qual_vars),
                 contains("BsmtFin"),
                 Functional,
                 GarageFinish),
            ~ as.numeric(.))

nominal_vars = c("Alley", "BldgType", "Condition1", "Condition2", "Electrical", 
                "Exterior1st", "Exterior2nd", "Fence", "Foundation", "GarageType", 
                "Heating", "HouseStyle", "LandContour", "LotConfig", "MSSubClass",
                "MSZoning", "MasVnrType", "MiscFeature", "MoSold", "Neighborhood",
                "RoofMatl", "RoofStyle", "SaleCondition", "SaleType", "YrSold",
                "CentralAir", "BsmtExposure", "LandSlope", "LotShape", "PavedDrive",
                "YearBuilt", "YearRemodAdd", "Street", "Utilities")

all_dat <-
  all_dat %>% 
  mutate_at(vars(all_of(nominal_vars)), as.character)

dummies <- model.matrix(~., data = all_dat %>% select_if(is.character)) %>% as_tibble() %>% select(-contains("Inter"))

all_dat <-
  all_dat %>% 
  select_if(negate(is.character)) %>% 
  bind_cols(dummies)

skewed_vars <- 
  map_lgl(all_dat %>% select_if(is.numeric), ~ (moments::skewness(.x) > 0.75)) %>% which() %>% names()

lmbdas <- map_dbl(skewed_vars, ~ car::powerTransform(all_dat[[.x]] + 1)$lambda)

trans_vars <- map2_dfc(skewed_vars, lmbdas, ~ car::bcPower(all_dat[[.x]] + 1, lambda = .y))

colnames(trans_vars) <- skewed_vars

all_dat <-
  all_dat %>% 
  select(-all_of(skewed_vars)) %>% 
  bind_cols(trans_vars)

glue("Lambda value for LowQualFinSF in R is: {car::powerTransform(train[['LowQualFinSF']] + 1)$lambda}")

glue("There are {sum(is.na(all_dat))} features with missing observations")

sparse <- all_dat %>% 
  summarize_if(is.numeric, ~ sum(. == 0)/nrow(all_dat)) %>% 
  pivot_longer(cols = everything()) %>% 
  filter(value > 0.9994) %>% 
  pull(name)

all_dat <- all_dat %>% 
  select(-one_of(sparse))

glue("The are {ncol(all_dat)} features and {nrow(all_dat)} observations.")

all_dat_scaled <-
  all_dat %>% 
  mutate_all(~ scale(.)[,1], center = TRUE)

train_dat <- all_dat %>% slice(1:nrow(train)) %>% janitor::clean_names()
train_dat_scaled <- all_dat_scaled %>% slice(1:nrow(train)) %>% janitor::clean_names()

test_dat <- all_dat %>% slice(nrow(train)+1:nrow(.)) %>% janitor::clean_names()
test_dat_scaled <- all_dat_scaled %>% slice(nrow(train)+1:nrow(.)) %>% janitor::clean_names()

# Add the log of the target back to the training data for modeling
train_dat[["SalePrice"]] <- log(train[["SalePrice"]])
train_dat_scaled[["SalePrice"]] <- log(train[["SalePrice"]])

#sz <- ceiling(nrow(X)*.7)
#
#set.seed(1491)
#train_idx <- sort(sample(nrow(X), sz))
#X_train <- X[train_idx, ]
#X_test <- X[-train_idx, ]
#
#y_train <- y[train_idx]
#y_test <- y[-train_idx]

set.seed(1491)
train_split <- initial_split(train_dat, prob = 0.7)
set.seed(1491)
train_scaled_split <- initial_split(train_dat_scaled, prob = 0.7)
train_cv <- training(train_split) %>% vfold_cv(v = 5)
train_scaled_cv <- training(train_scaled_split) %>% vfold_cv(v = 5)

glue("The training data is {nrow(training(train_split))} rows by {ncol(training(train_split))} columns")

glue("The testing data is {nrow(testing(train_split))} rows by {ncol(testing(train_split))} columns")

lasso_model <-
    linear_reg(
        penalty = tune(),
        mixture = 1) %>%
    set_engine("glmnet",
               standardize = TRUE)

alphas <- grid_regular(penalty(range = c(-5, 0)), levels = 25)
alphas_py <- alphas %>% pull(penalty)

lasso_cv <-
    tune_grid(
        formula = SalePrice ~ .,
        model = lasso_model,
        resamples = train_cv,
        grid = alphas,
        metrics = metric_set(rmse),
        control = control_grid(verbose = FALSE)
    )

best_lasso <-
    lasso_cv %>%
    select_best("rmse", maximize = FALSE)

print(best_lasso)

best_lasso_model <-
    lasso_model %>%
    finalize_model(parameters = best_lasso)

lasso_fit <-
    best_lasso_model %>%
    fit(SalePrice ~ ., training(train_split))

lasso_predictions <- predict(lasso_fit, testing(train_split))

testing(train_split) %>%
  select(SalePrice) %>%
  bind_cols(lasso_predictions) %>%
  rmse(SalePrice, .pred) %>% 
  pull(.estimate) %>% 
  {glue("For the LASSO model in R, the RMSE is {round(., 4)}")}

ridge_model <-
    linear_reg(
        penalty = tune(),
        mixture = 0) %>%
    set_engine("glmnet",
               standardize = TRUE)

alphas <- grid_regular(penalty(range = c(-5, 0)), levels = 25)
alphas_py <- alphas %>% pull(penalty)

ridge_cv <-
    tune_grid(
        formula = SalePrice ~ .,
        model = ridge_model,
        resamples = train_cv,
        grid = alphas,
        metrics = metric_set(rmse),
        control = control_grid(verbose = FALSE)
    )

best_ridge <-
    ridge_cv %>%
    select_best("rmse", maximize = FALSE)

print(best_ridge)

best_ridge_model <-
    ridge_model %>%
    finalize_model(parameters = best_ridge)

ridge_fit <-
    best_ridge_model %>%
    fit(SalePrice ~ ., training(train_split))

ridge_predictions <- predict(ridge_fit, testing(train_split))

testing(train_split) %>%
  select(SalePrice) %>%
  bind_cols(ridge_predictions) %>%
  rmse(SalePrice, .pred) %>% 
  pull(.estimate) %>% 
  {glue("For the ridge model in R, the RMSE is {round(., 4)}")}

en_model <-
    linear_reg(
        penalty = tune(),
        mixture = tune()) %>%
    set_engine("glmnet",
               standardize = TRUE)

en_grid <- grid_regular(penalty(range = c(-5, 0)), mixture(c(0.01, 0.95)), levels = 10)
alphas_py <- en_grid %>% pull(penalty) %>% unique() %>% round(5)
l1 <- en_grid %>% pull(mixture) %>% unique()  

en_cv <-
    tune_grid(
        formula = SalePrice ~ .,
        model = en_model,
        resamples = train_cv,
        grid = en_grid,
        metrics = metric_set(rmse),
        control = control_grid(verbose = FALSE)
    )

best_en <-
    en_cv %>%
    select_best("rmse", maximize = FALSE)

print(best_en)

best_en_model <-
    en_model %>%
    finalize_model(parameters = best_en)

en_fit <-
    best_en_model %>%
    fit(SalePrice ~ ., training(train_split))

en_predictions <- predict(en_fit, testing(train_split))

testing(train_split) %>%
  select(SalePrice) %>%
  bind_cols(en_predictions) %>%
  rmse(SalePrice, .pred) %>% 
  pull(.estimate) %>% 
  {glue("For the Elastic Net model in R, the RMSE is {round(., 4)}")}

mars_model <-
  mars(
    mode = "regression",
    prod_degree = tune(),
    prune_method = tune()
  ) %>%
  set_engine("earth")

mars_grid <-
    grid_regular(
        prod_degree(),
        prune_method(c('backward', 'none', 'exhaustive', 'forward', 'seqrep')),
        levels = 5
    )

mars_tune <-
  tune_grid(
    formula   = SalePrice ~ .,
    model     = mars_model,
    resamples = train_cv,
    grid      = mars_grid,
    metrics   = metric_set(rmse),
    control   = control_grid(verbose = FALSE)
  )

best_mars <-
  mars_tune %>%
  select_best("rmse", maximize = FALSE)

print(best_mars)

mars_fit <-
  mars_model %>%
  finalize_model(parameters = best_mars) %>%
  fit(SalePrice ~ ., training(train_split))

mars_predictions <-
  testing(train_split) %>%
  select(SalePrice) %>%
  bind_cols(
    predict(mars_fit, testing(train_split))
  )

testing(train_split) %>%
  select(SalePrice) %>%
  bind_cols(mars_predictions) %>%
  rmse(SalePrice, .pred) %>% 
  pull(.estimate) %>% 
  {glue("For the MARS model in R, the RMSE is {round(., 4)}")}

# R Parsnip to Python translator:
# sample_size = subsample
# loss_reduction = gamma
# trees = n_estimators
# tree_depth = max_depth
# min_n = min_child_weight
# mtry = colsample_bytree

#xgb_model <-
#  boost_tree(
#    mode = "regression",
#    trees = tune(),
#    learn_rate = tune(),
#    min_n = tune(),
#    tree_depth = tune()
#  ) %>%
#  set_engine(
#    "xgboost",
#    booster = "gbtree"
#  )
#
#xgb_grid <-
#    grid_regular(
#        trees(c(100, 1000)),
#        learn_rate(c(-4, 0)),
#        min_n(c(2, 10)),
#        tree_depth(c(2, 8)),
#        levels = 3
#    )
#
#learn_rate_py <- xgb_grid %>% pull(learn_rate) %>% unique()
#min_n_py <- xgb_grid %>% pull(min_n) %>% unique()
#tree_depth_py <- xgb_grid %>% pull(tree_depth) %>% unique()
#
#xgb_tune <-
#  tune_grid(
#    formula   = SalePrice ~ .,
#    model     = xgb_model,
#    resamples = train_cv,
#    grid      = xgb_grid,
#    metrics   = metric_set(rmse),
#    control   = control_grid(verbose = TRUE)
#  )
#
#best_xgb <-
#  xgb_tune %>%
#  select_best("rmse", maximize = FALSE)
#
#print(best_xgb)
#
#xgb_fit <-
#  xgb_model %>%
#  finalize_model(parameters = best_xgb) %>%
#  fit(SalePrice ~ ., training(train_split))

xgb_model <-
  boost_tree(
    mode = "regression",
    trees = 2000,
    sample_size = 0.5,
    mtry = 1,
    loss_reduction = 0.001,
    learn_rate = 0.01,
    min_n = 10,
    tree_depth = 5
  ) %>%
  set_engine(
    "xgboost",
    booster = "gbtree"
  )

xgb_fit <-
  xgb_model %>%
  fit(SalePrice ~ ., training(train_scaled_split))

xgb_predictions <-
  testing(train_split) %>%
  select(SalePrice) %>%
  bind_cols(
    predict(xgb_fit, testing(train_scaled_split))
  )

testing(train_scaled_split) %>%
  select(SalePrice) %>%
  bind_cols(xgb_predictions) %>%
  rmse(SalePrice, .pred) %>% 
  pull(.estimate) %>% 
  {glue("For the XGBoost model in R, the RMSE is {round(., 4)}")}