PROJ2_code.Rmd

---
title: "Project 2 Cloud Data Code Appendix"
author: "Yicheng Shen (Student ID: 2806571) & Yunhong Bao (Student ID: 2427527)"
date: "December 6, 2022"
header-includes: 
      - \usepackage{amsmath}
output: 
    pdf_document    
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = T, eval = F, cache = TRUE, warning = F, message = F)
library(mosaic)
library(GGally)
library(caret)
library(ggfortify)
library(gridExtra)
library(MASS)
library(e1071)
library(DiagrammeR)
library(class)
library(tree)
library(randomForest)
library(gbm)
library(xgboost)
ggplot2::theme_set(ggplot2::theme_bw())
knitr::opts_chunk$set(out.width = "100%", fig.align = 'center')

source("CVmaster.R")
```


## EDA 

*scale in regression*

```{r read data and give name}
image_1 <- read.table("image_data/imagem1.txt")
image_2 <- read.table("image_data/imagem2.txt")
image_3 <- read.table("image_data/imagem3.txt")
var_name <- c("Y_coord", "X_coord", "Cloud", "NDAI", "SD", "CORR", "DF","CF","BF","AF","AN")
colnames(image_1) = colnames(image_2) = colnames(image_3) <- var_name
image_1$Cloud <- factor(image_1$Cloud)
image_2$Cloud <- factor(image_2$Cloud)
image_3$Cloud <- factor(image_3$Cloud)
image_1$logSD <- log(image_1$SD)
image_2$logSD <- log(image_2$SD)
image_3$logSD <- log(image_3$SD)
```

```{r summary and dim}
summary(image_1)
summary(image_2)
summary(image_3)
dim(image_1)
dim(image_2)
dim(image_3)
```

```{r check na}
any(apply(image_1, 2, is.na))
any(apply(image_2, 2, is.na))
any(apply(image_3, 2, is.na))
```


```{r feature boxplot}
# image_1 %>% 
# ggplot() + geom_boxplot(aes(x=Cloud, y=SD)) 

all_image <- rbind(image_1 %>% mutate(image = "1"), image_2 %>% mutate(image = "2"), image_3 %>% mutate(image = "3"))
 
all_image%>% 
  ggplot() + geom_boxplot(aes(x=Cloud, y=NDAI, color = image)) 

all_image %>% 
  ggplot() + geom_boxplot(aes(x=Cloud, y=SD, color = image)) 

all_image %>% 
  ggplot() + geom_boxplot(aes(x=Cloud, y=CORR, color = image)) 
```


```{r pairwise, message = F, warning = F}
image_1 %>%
  dplyr::select(Cloud, NDAI, SD, CORR) %>% ggpairs()

image_2 %>%
  dplyr::select(Cloud, NDAI, SD, CORR) %>% ggpairs()

image_3 %>%
  dplyr::select(Cloud, NDAI, SD, CORR) %>% ggpairs()
```


*Do a PCA on the angle *
```{r angle corr}
image_1 %>% dplyr::select(DF:AN) %>% 
  ggpairs()
image_2 %>% dplyr::select(DF:AN) %>% 
  ggpairs()
image_3 %>% dplyr::select(DF:AN) %>% 
  ggpairs()
```


```{r cloud map}
# library(sf)
# ggplot(data = world) +
#     geom_sf() +
#     geom_point(data = image_1, aes(x = X-coord, y = Y-coord), size = 4, 
#         shape = 23, fill = "darkred") +
#     coord_sf(xlim = c(-88, -78), ylim = c(0, 33), expand = FALSE)

Palette_1 <- c("gray","black", "white")

map_1 <- image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Cloud), size = 0.5) +
  theme_map() + theme(legend.position = "none") +
  labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_1)

map_2 <- image_2 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Cloud), size =  0.5) +
  theme_map() + theme(legend.position = "none") +
  labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_1) 

map_3 <- image_3 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Cloud), size = 0.5) +
  theme_map() + theme(legend.position = "none") +
  labs(title = "Image 3", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_1)
```

```{r arranged map, fig.width = 12, fig.height = 4}
grid.arrange(arrangeGrob(map_1 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         map_2 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         map_3 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         nrow = 1))
```


*Better map*

```{r new map}
Palette_2 <- c("gray", "white")

map1 <- image_1 %>%
  filter(Cloud != 0) %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Cloud), size = 0.5) +
  theme_dark() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_2)

map2 <- image_2 %>%
  filter(Cloud != 0) %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Cloud), size = 0.5) +
  theme_dark() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_2)

map3 <- image_3 %>%
  filter(Cloud != 0) %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Cloud), size = 0.5) +
  theme_dark() + labs(title = "Image 3", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_2)
# map1
# map2
# map3
```

```{r arrange map, fig.width=12,fig.height=4}
map_legend <- lemon::g_legend(map1 +  guides(colour = guide_legend(nrow = 1)))

grid.arrange(arrangeGrob(map1 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         map2 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         map3 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         nrow = 1),
             map_legend, nrow = 2, heights = c(10, 1))
```


```{r image 1 covariate map, fig.width=12,fig.height=4}
NDAI_plot <- image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = NDAI), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate")

SD_plot <- image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = logSD), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate")

CORR_plot <- image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = CORR), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 


grid.arrange(NDAI_plot + ggeasy::easy_center_title(),
             SD_plot + ggeasy::easy_center_title(),
             CORR_plot + ggeasy::easy_center_title(),
             nrow = 1)
```


```{r image 2 covariate map, fig.width=12,fig.height=4}
NDAI_plot <- image_2 %>% 
  # filter (NDAI >= 1.2 )%>% filter (NDAI <= 2.5 )%>% 
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = NDAI), size = 0.5) +
  theme_map() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") 

SD_plot <- image_2 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = logSD), size = 0.5) +
  theme_map() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") 

CORR_plot <- image_2 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = CORR), size = 0.5) +
  theme_map() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") 

grid.arrange(NDAI_plot + ggeasy::easy_center_title(),
             SD_plot + ggeasy::easy_center_title(),
             CORR_plot + ggeasy::easy_center_title(),
             nrow = 1)
```


```{r image 3 covariate map, fig.width=12,fig.height=4}
NDAI_plot <- image_3%>% 
  # filter (NDAI >= 1.2 )%>% filter (NDAI <= 2.5 )%>% 
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = NDAI), size = 0.5) +
  theme_map() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") 

SD_plot <- image_3 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = logSD), size = 0.5) +
  theme_map() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") 

CORR_plot <- image_3 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = CORR), size = 0.5) +
  theme_map() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate")

grid.arrange(NDAI_plot + ggeasy::easy_center_title(),
             SD_plot + ggeasy::easy_center_title(),
             CORR_plot + ggeasy::easy_center_title(),
             nrow = 1)
```


```{r angle map}
image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = DF), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 

image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = CF), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 

image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = BF), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 

image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = AF), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 

image_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = AN), size = 0.5) +
  theme_map() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 
```

# Remove cloud = 0 from here.

```{r filter zero, eval = T, echo = T}
Image_1 <- image_1 %>% filter(Cloud != "0") %>% mutate(Cloud = factor(Cloud))
Image_2 <- image_2 %>% filter(Cloud != "0") %>% mutate(Cloud = factor(Cloud))
Image_3 <- image_3 %>% filter(Cloud != "0") %>% mutate(Cloud = factor(Cloud))

dim(Image_1)
dim(Image_2)
dim(Image_3)
```


```{r hist covariate, fig.width=12, fig.height=3}
hist_1 <- Image_1 %>% 
  ggplot(aes(x = NDAI)) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_2 <- Image_2 %>% 
  ggplot(aes(x = NDAI)) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_3 <- Image_3 %>% 
  ggplot(aes(x = NDAI)) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

grid.arrange(hist_1 + ggeasy::easy_center_title(),
             hist_2 + ggeasy::easy_center_title(),
             hist_3 + ggeasy::easy_center_title(),
             nrow = 1)

hist_1 <- Image_1 %>% 
  ggplot(aes(x = SD)) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_2 <- Image_2 %>% 
  ggplot(aes(x = SD)) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_3 <- Image_3 %>% 
  ggplot(aes(x = SD)) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

grid.arrange(hist_1 + ggeasy::easy_center_title(),
             hist_2 + ggeasy::easy_center_title(),
             hist_3 + ggeasy::easy_center_title(),
             nrow = 1)


hist_1 <- Image_1 %>% 
  ggplot(aes(x = log(SD))) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_2 <- Image_2 %>% 
  ggplot(aes(x = log(SD))) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_3 <- Image_3 %>% 
  ggplot(aes(x = log(SD))) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

grid.arrange(hist_1 + ggeasy::easy_center_title(),
             hist_2 + ggeasy::easy_center_title(),
             hist_3 + ggeasy::easy_center_title(),
             nrow = 1)



hist_1 <- Image_1 %>% 
  ggplot(aes(x = CORR )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_2 <- Image_2 %>% 
  ggplot(aes(x = CORR )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_3 <- Image_3 %>% 
  ggplot(aes(x = CORR )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

grid.arrange(hist_1 + ggeasy::easy_center_title(),
             hist_2 + ggeasy::easy_center_title(),
             hist_3 + ggeasy::easy_center_title(),
             nrow = 1)


All_Image <- rbind(Image_1, Image_2, Image_3)

hist_1 <- All_Image %>% 
  ggplot(aes(x = NDAI )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_2 <- All_Image %>% 
  ggplot(aes(x = CORR )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_3 <- All_Image %>% 
  ggplot(aes(x = SD )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_4 <- All_Image %>% 
  ggplot(aes(x = log(SD) )) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

grid.arrange(hist_1 + ggeasy::easy_center_title(),
             hist_2 + ggeasy::easy_center_title(),
             hist_3 + ggeasy::easy_center_title(),
             hist_4 + ggeasy::easy_center_title(),
             nrow = 2)
```



```{r angle hist, fig.width=12, fig.height=3.5}
hist_1 <- Image_1 %>% 
  ggplot(aes(x = (AF))) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_2 <- Image_2 %>% 
  ggplot(aes(x = (AF))) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

hist_3 <- Image_3 %>% 
  ggplot(aes(x = (AF))) +
  geom_density(aes(fill = Cloud, color = Cloud), alpha = 0.7) 

grid.arrange(hist_1 + ggeasy::easy_center_title(),
             hist_2 + ggeasy::easy_center_title(),
             hist_3 + ggeasy::easy_center_title(),
             nrow = 1)
```



# 2.Data Preparation 


# 2.a


#### Make a table of availble data: 


```{r index groups, eval = F, echo = F}
# set.seed(8848)
# sample_index_1 <- sample(nrow(Image_1), nrow(Image_1)*0.8)
# Image_1_Train <- Image_1[sample_index_1,]
# train_index_1 <- sample(nrow(Image_1_Train), nrow(Image_1_Train)*0.8)
# Image_1_val <- Image_1_Train[-train_index_1,]
# Image_1_train <- Image_1_Train[train_index_1,]
# Image_1_Test <- Image_1[-sample_index_1,]
# 
# 
# sample_index_2 <- sample(nrow(Image_2), nrow(Image_2)*0.8)
# Image_2_Train <- Image_2[sample_index_2,]
# train_index_2 <- sample(nrow(Image_2_Train), nrow(Image_2_Train)*0.8)
# Image_2_val <- Image_2_Train[-train_index_2,]
# Image_2_train <- Image_2_Train[train_index_2,]
# Image_2_Test <- Image_2[-sample_index_2,]
# 
# 
# sample_index_3 <- sample(nrow(Image_3), nrow(Image_3)*0.8)
# Image_3_Train <- Image_3[sample_index_3,]
# train_index_3 <- sample(nrow(Image_3_Train), nrow(Image_3_Train)*0.8)
# Image_3_val <- Image_3_Train[-train_index_3,]
# Image_3_train <- Image_3_Train[train_index_3,]
# Image_3_Test <- Image_3[-sample_index_3,]
```


```{r 5 groups, fig.width=12, fig.height=4}
# summary(Image_1_train[,1:2])
# quantile(Image_1_train$X_coord, c(0.5))
# quantile(Image_1_train$Y_coord, c(0.2, 0.4, 0.6, 0.8))

seperate_fold <- function(train_data) {
  
  fold_1 <- train_data %>% filter( Y_coord  < 383*0.2 )
  fold_2 <- train_data %>% filter( 383*0.2 <= Y_coord, Y_coord  < 383*0.4 )
  fold_3 <- train_data %>% filter( 383*0.4 <= Y_coord, Y_coord  < 383*0.6 )
  fold_4 <- train_data %>% filter( 383*0.6 <= Y_coord, Y_coord  < 383*0.8 )
  fold_5 <- train_data %>% filter( 383*0.8 <= Y_coord )
  
return( rbind(fold_1 %>% mutate(Type = "Train"),
      fold_2 %>% mutate(Type = "Train"),
      fold_3 %>% mutate(Type = "Train"),
      fold_4 %>% mutate(Type = "Test"),
      fold_5 %>% mutate(Type = "Validation") ) ) 
}

Cut_1 <- seperate_fold(Image_1)
Cut_2 <- seperate_fold(Image_2)
Cut_3 <- seperate_fold(Image_3)


Cut_1_Map <- Cut_1 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Type), size = 0.5) +
  # geom_hline(yintercept = 383*0.2 ) + 
  # geom_hline(yintercept = 383*0.4) + 
  labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 

Cut_2_Map <- Cut_2 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Type), size = 0.5) +
  # geom_hline(yintercept = 383*0.2) + 
  # geom_hline(yintercept = 383*0.4) + 
  labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate")

Cut_3_Map <- Cut_3 %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Type), size = 0.5) +
  # geom_hline(yintercept = 383*0.2 ) + 
  # geom_hline(yintercept = 383*0.4) + 
  labs(title = "Image 3", x = "X Coordinate", y = "Y Coordinate") + theme(legend.position="none")


grid.arrange(arrangeGrob(Cut_1_Map + ggeasy::easy_center_title() + theme(legend.position="none"),
                         Cut_2_Map + ggeasy::easy_center_title() + theme(legend.position="none"),
                         Cut_3_Map + ggeasy::easy_center_title() + theme(legend.position="none"),
                         nrow = 1),
             lemon::g_legend(Cut_1_Map +  guides(colour = guide_legend(nrow = 1))), nrow = 2, heights = c(10, 1))

Image_all_Cut <- rbind(Cut_1, Cut_2, Cut_3)

prop.table(table(Image_all_Cut$Type))

Image_all_Cut %>% 
  group_by(Type) %>% 
  dplyr::summarise(`Perct Cloud` = mean(Cloud == "1"))
```


```{r 5 KM groups, fig.width=12, fig.height=4}
set.seed(8848)

KM_1 <- kmeans(Image_1[,1:2], centers = 5, nstart = 25)
KM_2 <- kmeans(Image_2[,1:2], centers = 5, nstart = 25)
KM_3 <- kmeans(Image_3[,1:2], centers = 5, nstart = 25)

Image_all_KM <- rbind(Image_1 %>% mutate(image = "1", cluster = KM_1$cluster),
                      Image_2 %>% mutate(image = "2", cluster = KM_2$cluster),
                      Image_3 %>% mutate(image = "3", cluster = KM_3$cluster))

Image_all_KM <- Image_all_KM %>% 
  mutate(Type = ifelse(cluster == "2", "Test", "Train"))  %>% 
  mutate(Type = ifelse(cluster == "5", "Validation", Type)) 

KM_1_Map <- Image_all_KM %>% filter(image == "1") %>% 
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Type), size = 0.5) +
  theme_bw() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate") 
# +  scale_color_manual(values = c("red", "blue", "green", "green", "green"))  

KM_2_Map <- Image_all_KM %>% filter(image == "2") %>% 
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Type), size = 0.5) +
  theme_bw() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") 

KM_3_Map <-Image_all_KM %>% filter(image == "3") %>% 
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = Type), size = 0.5) +
  theme_bw() + labs(title = "Image 3", x = "X Coordinate", y = "Y Coordinate") 

grid.arrange(arrangeGrob(KM_1_Map + ggeasy::easy_center_title() + theme(legend.position="none"),
                         KM_2_Map + ggeasy::easy_center_title() + theme(legend.position="none"),
                         KM_3_Map + ggeasy::easy_center_title() + theme(legend.position="none"),
                         nrow = 1),
             lemon::g_legend(KM_1_Map +  guides(colour = guide_legend(nrow = 1))), nrow = 2, heights = c(10, 1))


prop.table(table(Image_all_KM$Type))

Image_all_KM %>% 
  group_by(Type) %>% 
  dplyr::summarise(`Perct Cloud` = mean(Cloud == "1"))
```

# 2.b

```{r accuaracy of trival predictor}
cut_v <- Image_all_Cut %>% 
  filter(Type == "Validation")%>% 
  mutate(Predict = "-1") %>% 
  summarise(`Data Type` = "Validation", Accuracy = mean(Predict == Cloud))

cut_t <- Image_all_Cut %>% 
  filter(Type == "Test")%>% 
  mutate(Predict = "-1") %>% 
  summarise(`Data Type` = "Test", Accuracy = mean(Predict == Cloud))


km_v <- Image_all_KM %>% 
  filter(Type == "Validation")%>% 
  mutate(Predict = "-1") %>% 
  summarise(`Data Type` = "Validation", Accuracy = mean(Predict == Cloud))

km_t <- Image_all_KM %>% 
  filter(Type == "Test")%>% 
  mutate(Predict = "-1") %>% 
  summarise(`Data Type` = "Test", Accuracy = mean(Predict == Cloud))

knitr::kable(rbind(cut_v, cut_t, km_v, km_t),
    caption = "Proportions of Cloudy and Clear Surfaces by Expert Labeling",
    booktabs = T ) %>%
  kable_styling(latex_options = "HOLD_position", position = "center", font_size = 8)
```
 
 
# 2.c

```{r}
Image_all_Cut <- Image_all_Cut %>% dplyr::select(-SD)
Image_all_KM <- Image_all_KM %>% dplyr::select(-SD)
```


```{r}
# save(Image_all_Cut, file = "Image_all_Cut.RData")
# save(Image_all_KM, file = "Image_all_KM.RData")

load("Image_all_Cut.RData")
load("Image_all_KM.RData")
```


```{r}
pca_result <- prcomp(Image_all_Cut %>% dplyr::select(NDAI:logSD), cetner = T, scale. = T) 
as.matrix(pca_result$rotation)

eigenvalues <- (pca_result$sdev)^2
eigs_cum <- cumsum(eigenvalues)/sum(eigenvalues)
ggplot() + geom_point(aes(x = 1:length(eigenvalues), y=eigs_cum)) +
geom_line(aes(x = 1:length(eigenvalues), y=eigs_cum)) + theme_bw() +
  labs(x = "First PCs", y = "Fraction of total variance explained")       

cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$NDAI)
cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$logSD)
cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$CORR)

cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$DF)
cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$CF)
cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$BF)
cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$AF)
cor( as.numeric(Image_all_Cut$Cloud), Image_all_Cut$AN)
```


```{r PCA visual}
trn_pca <- Image_1_train %>% 
  filter(Cloud != "0") %>% 
  dplyr::select(-Cloud)  %>% 
  FactoMineR::PCA(graph = F)

proj <- trn_pca$ind$coord[, 1:2]  %>% 
  data.frame()  %>% 
  tibble()  %>% 
  mutate(Cloud_True = (Image_1_train %>% 
  filter(Cloud != "0") )$Cloud)

p1 <- proj %>% 
  ggplot(aes(x = Dim.1, y = Dim.2)) +
  geom_point(aes(color = Cloud_True)) +
  labs(x = "PC1", y = "PC2", title = "PCA Projections") 

p2 <- plot(trn_pca, choix = "var")

p1
```

Best are NDAI, CORR and logSD 


# 2d CV Master

Refer to the R file. 


# 3. Model Fitting & Validation

# 3.a

```{r train and test}
load("Image_all_Cut.RData")
load("Image_all_KM.RData")
Image_Train_Cut <- Image_all_Cut %>%  filter(Type != "Test")
Image_Test_Cut <- Image_all_Cut %>%  filter(Type == "Test")
Image_Train_KM <- Image_all_KM %>%  filter(Type != "Test")
Image_Test_KM <- Image_all_KM %>%  filter(Type == "Test")
```


### call CVmaster from separate R file

```{r, warning=F}
set.seed(8848)

# Cut
cv_result = CVmaster(
  classifier = "Logistic",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, c(4:11)],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

glm_test <- glm(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
        data = Image_Train_Cut,family = 'binomial')
probt <- predict(glm_test, Image_Train_Cut, type = "response")
t = pROC::coords(pROC::roc(Image_Train_Cut$Cloud, probt), "best", transpose = FALSE)

prob <- predict(glm_test, Image_Test_Cut, type = "response")
prediction <- ifelse(prob > t[[1]], 1,-1)
test_accuracy <- mean(prediction == Image_Test_Cut$Cloud)

result_logit_Cut <- cbind(classifier = "Logistic", data = "Horizontal Cut", cv_result, `Test Accuracy` = test_accuracy)

# KM
cv_result <- CVmaster(
  classifier = "Logistic",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

glm_test <- glm(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
        data = Image_Train_KM,family = 'binomial')
probt <- predict(glm_test, Image_Train_KM, type = "response")
t = pROC::coords(pROC::roc(Image_Train_KM$Cloud, probt), "best", transpose = FALSE)

prob <- predict(glm_test, Image_Test_KM, type = "response")
prediction <- ifelse(prob > t[[1]], 1,-1)

result_logit_KM <- cbind(classifier = "Logistic", data = "K-means", cv_result, `Test Accuracy` = test_accuracy) 


rbind(result_logit_Cut, result_logit_KM)  %>% 
  as_tibble() %>% 
  mutate_at(c("fold 1","fold 2","fold 3","fold 4","fold 5","fold 6","fold 7","fold 8","fold 9","fold 10","CV Average", "Test Accuracy"), as.numeric) %>% 
  mutate(across(where(is.numeric), ~ round(., 4)))
```


```{r, warning=F}
# Cut
set.seed(8848)
cv_result <- CVmaster(
  classifier = "LDA",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

lda_test <- lda(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
    data = Image_Train_Cut)
prediction <- predict(lda_test, Image_Test_Cut)$class
test_accuracy <- mean(prediction == Image_Test_Cut$Cloud)

result_lda_Cut <- cbind(classifier = "LDA", data = "Horizontal Cut", cv_result, `Test Accuracy` = test_accuracy)

# KM
cv_result <- CVmaster(
  classifier = "LDA",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

lda_test <- lda(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
    data = Image_Train_KM)
prediction = predict(lda_test, Image_Test_KM)$class
test_accuracy <- mean(prediction == Image_Test_KM$Cloud)

result_lda_KM <- cbind(classifier = "LDA", data = "K-means", cv_result, `Test Accuracy` = test_accuracy)

rbind(result_lda_Cut, result_lda_KM)  %>% 
  as_tibble() %>% 
  mutate_at(c("fold 1","fold 2","fold 3","fold 4","fold 5","fold 6","fold 7","fold 8","fold 9","fold 10","CV Average", "Test Accuracy"), as.numeric) %>% 
  mutate(across(where(is.numeric), ~ round(., 4)))
```



```{r, warning=F}
set.seed(8848)

# Cut
cv_result <- CVmaster(
  classifier = "QDA",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

qda_test <- qda(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
    data = Image_Train_Cut)
prediction <- predict(qda_test, Image_Test_Cut)$class
test_accuracy <- mean(prediction == Image_Test_Cut$Cloud)

result_qda_Cut <- cbind(classifier = "QDA", data = "Horizontal Cut", cv_result, `Test Accuracy` = test_accuracy)

# KM
cv_result <- CVmaster(
  classifier = "QDA",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

qda_test <- qda(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
    data = Image_Train_KM)
prediction <- predict(qda_test, Image_Test_KM)$class
test_accuracy <- mean(prediction == Image_Test_KM$Cloud)

result_qda_KM <- cbind(classifier = "QDA", data = "K-means", cv_result, `Test Accuracy` = test_accuracy)


rbind(result_qda_Cut, result_qda_KM)  %>% 
  as_tibble() %>% 
  mutate_at(c("fold 1","fold 2","fold 3","fold 4","fold 5","fold 6","fold 7","fold 8","fold 9","fold 10","CV Average", "Test Accuracy"), as.numeric) %>% 
  mutate(across(where(is.numeric), ~ round(., 4)))
```


```{r, warning=F}
set.seed(8848)

# Cut
cv_result <- CVmaster(
  classifier = "NB",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

NB_test <- naiveBayes(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
    data = Image_Train_Cut)
prediction <- predict(NB_test, Image_Test_Cut)
test_accuracy <- mean(prediction == Image_Test_Cut$Cloud)

result_NB_Cut <- cbind(classifier = "Naive Bayes", data = "Horizontal Cut", cv_result, `Test Accuracy` = test_accuracy)

# KM
cv_result <- CVmaster(
  classifier = "NB",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

NB_test <- naiveBayes(Cloud ~ NDAI + CORR + logSD + DF + CF + BF + AF + AN,
    data = Image_Train_KM)
prediction <- predict(NB_test, Image_Test_KM)
test_accuracy <- mean(prediction == Image_Test_KM$Cloud)


result_NB_KM <- cbind(classifier = "Naive Bayes", data = "K-means", cv_result, `Test Accuracy` = test_accuracy)


rbind(result_NB_Cut, result_NB_KM)  %>% 
  as_tibble() %>% 
  mutate_at(c("fold 1","fold 2","fold 3","fold 4","fold 5","fold 6","fold 7","fold 8","fold 9","fold 10","CV Average", "Test Accuracy"), as.numeric) %>% 
  mutate(across(where(is.numeric), ~ round(., 4)))
```


```{r, warning=F}
set.seed(8848)
CVmaster(
  classifier = "KNN",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)


CVmaster(
  classifier = "KNN",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)
```


```{r, warning=F}
set.seed(8848)
CVmaster(
  classifier = "SVM",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)


CVmaster(
  classifier = "SVM",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)
```


# ROC 

```{r}
source("ROC.R")
```

```{r}
glm_result = glm(Cloud ~ NDAI + logSD + CORR + DF + CF + BF + AF + AN , data=Image_Train_Cut, 
family = 'binomial')

preds <- predict(glm_result,  type = "response")

pROC::coords(roc(Image_Train_Cut$Cloud, preds), "best", best.method=c("closest.topleft"), transpose = FALSE)

preds <- predict(glm_result, Image_Test_Cut,  type = "response")

# pROC::coords(roc(Image_Test_Cut$Cloud, preds), "best", best.method=c("closest.topleft"), transpose = FALSE)

roc(Image_Test_Cut$Cloud, preds, print.auc = T, plot = T)


prediction <- ifelse(preds > 0.3575289, 1, -1)
test_accuracy <- mean(prediction == Image_Test_Cut$Cloud)


# KM

glm_result = glm(Cloud ~ NDAI + logSD + CORR + DF + CF + BF + AF + AN , data=Image_Train_KM, 
family = 'binomial')

preds <- predict(glm_result,  type = "response")

pROC::coords(roc(Image_Train_KM$Cloud, preds), "best", best.method=c("closest.topleft"), transpose = FALSE)

# pROC::coords(roc(Image_Test_Cut$Cloud, preds), "best", best.method=c("closest.topleft"), transpose = FALSE)

preds <- predict(glm_result, Image_Test_KM,  type = "response")

# pROC::coords(roc(Image_Test_KM$Cloud, preds), "best", best.method=c("closest.topleft"), transpose = FALSE)

roc(Image_Test_KM$Cloud, preds, print.auc = T, plot = T)

prediction <- ifelse(preds > 0.3700108, 1, -1)
test_accuracy <- mean(prediction == Image_Test_KM$Cloud)
```

```{r, out.width="50%"}
glm_result = glm(Cloud ~ NDAI + logSD + CORR + DF + CF + BF + AF + AN , data=Image_Train_Cut, 
family = 'binomial')

preds <- predict(glm_result, type = "response")

t = pROC::coords(roc(Image_Train_Cut$Cloud, preds), "best", transpose = FALSE)

preds <- predict(glm_result, Image_Test_Cut, type = "response")

roc(Image_Test_Cut$Cloud, preds, print.auc = T, plot = T, print.thres = t[[1]], main = "Logistic Regression")


glm_result = glm(Cloud ~ NDAI + logSD + CORR + DF + CF + BF + AF + AN , data=Image_Train_KM, 
family = 'binomial')

preds <- predict(glm_result, type = "response")

t = pROC::coords(roc(Image_Train_KM$Cloud, preds), "best", transpose = FALSE)

preds <- predict(glm_result, Image_Test_KM, type = "response")

roc(Image_Test_KM$Cloud, preds, print.auc = T, plot = T, print.thres = t[[1]], main = "Logistic Regression")
```


```{r}
library(pROC)

roc_logreg <- ROC(
  classifier = "Logistic",
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  test_data = Image_Test_Cut
)

plot(roc_logreg)


roc_LDA <- ROC(
  classifier = "LDA",
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3]
)

plot(roc_LDA)

roc_QDA <- ROC(
  classifier = "QDA",
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3]
)

plot(roc_QDA)


roc_NB <- ROC(
  classifier = "NB",
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3]
)

plot(roc_NB)
```


```{r}
roc_knn <- list()

for (i in 2:10) {
roc_knn[[i]] <- ROC(
  classifier = "KNN",
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = i
)
}
```


# Boosted tree

```{r}
Image_Train_Cut$Cloud01 <- factor(ifelse(Image_Train_Cut$Cloud == "-1", "0", Image_Train_Cut$Cloud))
levels(Image_Train_Cut$Cloud01) <- c("0", "1")
Image_Test_Cut$Cloud01 <- factor(ifelse(Image_Test_Cut$Cloud == "-1", "0", Image_Test_Cut$Cloud))
levels(Image_Test_Cut$Cloud01) <- c("0", "1")

Image_Train_KM$Cloud01 <- factor(ifelse(Image_Train_KM$Cloud == "-1", "0", Image_Train_KM$Cloud))
levels(Image_Train_KM$Cloud01) <- c("0", "1")
Image_Test_KM$Cloud01 <- factor(ifelse(Image_Test_KM$Cloud == "-1", "0", Image_Test_KM$Cloud))
levels(Image_Test_KM$Cloud01) <- c("0", "1")
```

```{r, warning=F, cache = T}
source("CVmaster.R")

# cut 
cv_result <- CVmaster(
  classifier = "Boosting Tree",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, 4:11],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

boost_model <- xgboost(
  data = as.matrix(Image_Train_Cut[, 4:11]),
  label = as.matrix(Image_Train_Cut$Cloud01),
  max.depth = 4,
  eta = 0.05,
  nthread = parallel::detectCores(),
  nrounds = 1000,
  objective = "binary:logistic"
)

pred <- predict(boost_model, as.matrix(Image_Train_Cut[, 4:11]))
pROC::coords(roc(Image_Train_Cut$Cloud01, pred), "best", transpose = FALSE)

pred <- predict(boost_model, as.matrix(Image_Test_Cut[, 4:11]))

# Cut ROC
roc(Image_Test_Cut$Cloud01, pred, print.auc = T, plot = T, print.thres = 0.4955721)
# pROC::coords(roc(Image_Test_Cut$Cloud01, pred), "best", best.method=c("closest.topleft"), transpose = FALSE)


prediction <- as.numeric(pred > 0.5)
test_accuracy <- mean(prediction == Image_Test_Cut$Cloud01)


# KM
cv_result <- CVmaster(
  classifier = "Boosting Tree",
  train_coord = Image_Train_KM[, 1:2],
  train_feature = Image_Train_KM[, 4:11],
  train_label = Image_Train_KM[, 3],
  K = 10,
  loss_fun = "Accuracy"
)

boost_model <- xgboost(
  data = as.matrix(Image_Train_KM[, 4:11]),
  label = as.matrix(Image_Train_KM$Cloud01),
  max.depth = 4,
  eta = 0.01,
  nthread = parallel::detectCores(),
  nrounds = 3,
  objective = "binary:logistic"
)

pred <- predict(boost_model, as.matrix(Image_Train_KM[, 4:11]))
pROC::coords(roc(Image_Train_KM$Cloud01, pred), "best", transpose = FALSE)


pred <- predict(boost_model, as.matrix(Image_Test_KM[, 4:11]))

# KM ROC
roc_test_boosting <- roc(Image_Test_KM$Cloud01, pred)
roc(Image_Test_KM$Cloud01, pred, print.auc = T, plot = T, print.thres = 0.4980169)
# pROC::coords(roc(Image_Test_KM$Cloud01, pred), "best", best.method=c("closest.topleft"), transpose = FALSE)


prediction <- as.numeric(pred > 0.4980169)
test_accuracy <- mean(prediction == Image_Test_KM$Cloud01)

result_BT_KM <- cbind(classifier = "Boosting Trees", data = "K-means", cv_result, `Test Accuracy` = test_accuracy)
```

```{r}
set.seed(8848)
boost_model <- xgboost(
  data = as.matrix(Image_Train_KM[, 4:11]),
  label = as.matrix(Image_Train_KM$Cloud01),
  max.depth = 2,
  eta = 0.05,
  nthread = parallel::detectCores(),
  nrounds = 1000,
  objective = "binary:logistic",
  colsample_bytree=0.6,
  verbose=0
)

pred <- predict(boost_model, as.matrix(Image_Train_KM[, 4:11]))
t=pROC::coords(roc(Image_Train_KM$Cloud01, pred), "best", transpose = FALSE)

pred <- predict(boost_model, as.matrix(Image_Test_KM[, 4:11]))

# Cut ROC
# roc(Image_Test_Cut$Cloud01, pred, print.auc = T, plot = T, print.thres = 0.49478)
# pROC::coords(roc(Image_Test_Cut$Cloud01, pred), "best", best.method=c("closest.topleft"), transpose = FALSE)

prediction <- as.numeric(pred > t[[1]])
test_accuracy <- mean(prediction == Image_Test_KM$Cloud01)

summary(boost_model)


# var importance plot
xgb.importance(colnames(as.matrix(Image_Train_KM[, 4:11])), model = boost_model)
xgb.plot.importance(xgb.importance(colnames(as.matrix(Image_Train_KM[, 4:11])), model = boost_model))
```



```{r}
# only test accuracy
library(xgboost)
set.seed(8848)

# Tested multiple times, deptth = 4 is the best
# Cut
boost_model <- xgboost(
  data = as.matrix(Image_Train_Cut[, 4:11]),
  label = as.matrix(Image_Train_Cut$Cloud01),
  max.depth = 4,
  eta = 0.01,
  nthread = parallel::detectCores(),
  nrounds = 3,
  objective = "binary:logistic"
)


pred <- predict(boost_model, as.matrix(Image_Test_Cut[, 4:11]))
prediction <- as.numeric(pred > 0.5)
mean(prediction == Image_Test_Cut$Cloud01)


# KM
boost_model <- xgboost(
  data = as.matrix(Image_Train_KM[, 4:11]),
  label = as.matrix(Image_Train_KM$Cloud01),
  max.depth = 4,
  eta = 0.01,
  nthread = parallel::detectCores(),
  nrounds = 3,
  objective = "binary:logistic"
)


pred <- predict(boost_model, as.matrix(Image_Test_KM[, 4:11]))
prediction <- as.numeric(pred > 0.5)
mean(prediction == Image_Test_KM$Cloud01)
```

```{r}
# boosting prediction plot
Image_all_KM$Cloud01 <- factor(ifelse(Image_all_KM$Cloud == "-1", "0", Image_all_KM$Cloud))
levels(Image_all_KM$Cloud01) <- c("0", "1")

boost_model <- xgboost(
  data = as.matrix(Image_all_KM[, 4:11]),
  label = as.matrix(Image_all_KM$Cloud01),
  max.depth = 4,
  eta = 0.01,
  nthread = parallel::detectCores(),
  nrounds = 3,
  objective = "binary:logistic"
)

all_pred <- predict(boost_model, as.matrix(Image_all_KM[, 4:11]))
Image_all_KM$all_prediction <- as.factor(as.numeric(all_pred > 0.5))

# roc_boosting <- roc(Image_all_KM$Cloud01, all_pred)
# plot(roc_boosting)
# pROC::coords(roc_boosting, best.method=c("youden"))


pROC::coords(roc_boosting, "best", transpose = FALSE)
```

```{r, fig.width = 12, fig.height = 4, out.width="85%"}
# draw plots
Palette_2 <- c("gray", "white")

map1 <- Image_all_KM %>%
  filter(image == "1") %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = all_prediction), size = 0.5) +
  theme_dark() + labs(title = "Image 1", x = "X Coordinate", y = "Y Coordinate", color = "Prediction") +
  scale_colour_manual(values = Palette_2)

map2 <- Image_all_KM %>%
  filter(image == "2") %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = all_prediction), size = 0.5) +
  theme_dark() + labs(title = "Image 2", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_2)

map3 <- Image_all_KM %>%
  filter(image == "3") %>%
  ggplot() + geom_point(aes(x = X_coord, y = Y_coord, color = all_prediction), size = 0.5) +
  theme_dark() + labs(title = "Image 3", x = "X Coordinate", y = "Y Coordinate") +
  scale_colour_manual(values = Palette_2)

map_legend <- lemon::g_legend(map1 +  guides(colour = guide_legend(nrow = 1)))

grid.arrange(arrangeGrob(map1 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         map2 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         map3 + ggeasy::easy_center_title() + theme(legend.position="none"),
                         nrow = 1),
             map_legend, nrow = 2, heights = c(10, 1))
```

# Old boost code

```{r}
gbm_train_plot = function(d = 2){
# note that we need to make sure gbm labels are between 0 and 1, and integers
  gbm_cls <- gbm(
    as.integer(Cloud01) - 1 ~ NDAI + logSD + CORR + DF + CF + BF + AF + AN,
    data = Image_Train_Cut,
    interaction.depth = d,
    distribution = "bernoulli",
    n.trees = 500
  )
summary(gbm_cls)
Bs <- seq(500) * 10
test_err <- rep(0, length(Bs))
for (i in seq(500)){
predicted_local <- predict(gbm_cls, Image_Test_Cut, "response",
n.trees=Bs[i])
predicted_y <- as.integer(predicted_local > 0.5)
true_y <- as.integer(Image_Test_Cut$Cloud01) -1
test_err[i] <- 1 - sum(predicted_y==true_y)/nrow(Image_Test_Cut)
}
gbm_err_df <- data.frame(Bs = Bs, test_err=test_err)
ggplot(data = gbm_err_df) +
geom_line(aes(x=Bs, y=test_err))+ ylim(0.1, 0.25)
}
gbm_train_plot(d=2)
```

## All result

```{r}
all_result <- rbind(result_logit_Cut, result_logit_KM,
      result_lda_Cut, result_lda_KM,
      result_qda_Cut, result_qda_KM,
      result_NB_Cut, result_NB_KM
      # ,result_BT_Cut, result_BT_KM
      )  %>% 
  as_tibble() %>% 
  mutate_at(c("fold 1","fold 2","fold 3","fold 4","fold 5","fold 6","fold 7","fold 8","fold 9","fold 10","CV Average", "Test Accuracy"), as.numeric) %>% 
  mutate(across(where(is.numeric), ~ round(., 4)))

# write.csv(all_result, file = "all_result.csv", row.names = F)
```

```{r}
read.csv("all_result.csv")

knitr::kable(all_result,
    caption = "CV Results and Test Accuracy based on Two ways of Data Splitting",
    booktabs = T ) %>%
  kableExtra::kable_styling(latex_options = "HOLD_position", position = "center", font_size = 8)
```



# Model thoughts

```{r logistic}

train_feature <- Image_Train_Cut[, 4:6]
train_label <- Image_Train_Cut[, 3]

fold_index = kmeans(train_feature, centers = 5, nstart = 25)
dat = cbind(train_feature, train_label)
dat$fold_index = fold_index$cluster

dat_CV_train = dat %>% filter(fold_index != 1)
dat_CV_test = dat %>% filter(fold_index == 1)

M1.glm <- glm(paste("train_label ~ ", paste(names(train_feature), collapse = "+" ), sep = "" ), data = dat_CV_train, family = 'binomial')

prob <- predict(M1.glm, dat_CV_test , type="response")

prediction <- ifelse(prob > 0.5, 1, 0)
mean(prediction != dat_CV_test$train_label)
```


```{r LDA}
LDA_model <-
  lda(Cloud ~ NDAI + SD + CORR,
    data = Image_1_train
  )
mean(predict(LDA_model, Image_1_test)$class != Image_1_test$Cloud)

# klaR::partimat(Cloud ~ NDAI + SD + CORR, data = Image_1_train, method = "lda")
```


```{r QDA}
QDA_model <-
  lda(Cloud ~ NDAI + SD + CORR,
    data = Image_1_train
  )
mean(predict(QDA_model, Image_1_test)$class != Image_1_test$Cloud)
```


```{r NB}
naive_bayes <- naiveBayes(Cloud ~ NDAI + SD + CORR, 
                data=Image_1_train)
mean(predict(naive_bayes, Image_1_test) != Image_1_test$Cloud)
```


## Testing

```{r}
final_boost_All <- xgboost(
  data = as.matrix(Image_Train_Cut[, 4:11]),
  label = as.matrix(Image_Train_Cut$Cloud01),
  max.depth = 4,
  eta = 0.05,
  nrounds = 3000,
  verbose = F, nthread = parallel::detectCores(),
  objective = "binary:logistic",
  colsample_bytree = 1,
)
xgb.plot.tree(model = final_boost_All, trees = 1000, plot_width = 1000, plot_height = 1000)

```

```{r}
library(MLmetrics)
source("CVmaster.R")
cv_result = CVmaster(
  classifier = "Logistic",
  train_coord = Image_Train_Cut[, 1:2],
  train_feature = Image_Train_Cut[, c(4:11)],
  train_label = Image_Train_Cut[, 3],
  K = 10,
  loss_fun = "F1 Score"
)
```