########## INFORMATIONS ##########
# This file is made to study online classification
# So, each pair (wireless,wakeupfor) has its classification models (knn and decision tree)
##################################

library("tidyverse")
options(dplyr.summarise.inform = FALSE)
library("class")
library("rpart")
library("rpart.plot")
library("viridis")
library("MLmetrics")

## Simulation Parameters:
## simkey         {baseline,extended,hint,hintandextended}
## wireless       {lora,nbiot}
## wakeupfor      {60s,180s}
## seed           [1,200]
## node           on[0,12]
## isSender       {0,1}
## dataSize       {1MB}

## Metrics:
## energy         [0,+inf)
## nDataRcv       [0,+inf)

nseed=200
nwakeupfor=2
nwireless=2
nsimkey=4
nsimulations=nseed*nwakeupfor*nwireless*nsimkey # Must be 3200

## Load data
data=suppressMessages(read_csv("../CCGRID2022.csv"))%>%distinct() # Note that in the data experiment wireless=="lora",seed==1,wakeupfor==60,simkey=="baseline" is present 2 times in the CSV file
data_seed=data%>%group_by(simkey,wireless,wakeupfor,seed)%>%summarize(energy=sum(energy),coverage=sum(nDataRcv))%>%
    mutate(efficiency=energy/coverage)%>%
    ungroup()

F1_Score2=function(truth, pred){
    result=sapply(c("baseline","extended","hint","hintandextended"),function(c){
        cur_truth=truth[truth==c]
        cur_pred=pred[truth==c]
        col=paste0("f1_",c)
        score=F1_Score(cur_truth,cur_pred)
        if(is.nan(score)){score=0}
        list(tibble(!!col:=score))
    })
    do.call("cbind",result)
}

generate_accuracy_for=function(ignore_hint=FALSE,seed_max=200,attempts_max=2,wrl="lora",wuf=180) {
    attempts=seq(1,attempts_max)
    results=sapply(attempts,function(attempt){
        ## Prepare data for traning
        set.seed(1+attempt) # Reproducibility
        wireless_map=c("lora"=1,"nbiot"=2)
        cur_data_seed=data_seed%>%filter(wakeupfor==wuf,wireless==wrl)
        data_ml=cur_data_seed%>%select(-efficiency)%>%mutate(wireless=wireless_map[cur_data_seed$wireless])
        if(ignore_hint){
            data_ml=data_ml%>%filter(simkey!="hint")
        }
        train_set=data_ml%>%filter(seed<=seed_max)%>%select(-seed)     # train data on seed_max*3 days
        test_set=data_ml%>%suppressMessages(anti_join(train_set))%>%select(-seed)        # build test_sed excluding training set

        ## KNN training
        knn_predictions=knn(train=train_set%>%select(-simkey),test=test_set%>%select(-simkey),cl=train_set$simkey,k=min(10,NROW(train_set)))
        ## KNN analysis
        knn_cont_table=table(knn_predictions,test_set$simkey)
        knn_accuracy=(sum(diag(knn_cont_table)/sum(rowSums(knn_cont_table))))
        knn_prop_table=round(prop.table(knn_cont_table),digits=2)
        knn_f1_score=F1_Score2(test_set$simkey,knn_predictions)
        
        ## Decision tree
        tree=rpart(
            simkey ~ wireless + wakeupfor + energy + coverage,
            data=train_set,
            method="class",
            minsplit=60,
            minbucket=1)
        tree_predictions=predict(tree,newdata=test_set%>%select(-simkey),type="class")
        tree_cont_table=table(tree_predictions,test_set$simkey)
        tree_accuracy=(sum(diag(tree_cont_table)/sum(rowSums(tree_cont_table))))
        tree_prop_table=round(prop.table(tree_cont_table),digits=2)
        tree_f1_score=F1_Score2(test_set$simkey,tree_predictions)

        ## Format data
        result_data=tibble(seed_max=seed_max,model=c("knn","tree"),accuracy=c(knn_accuracy,tree_accuracy))
        result_data=cbind(result_data,rbind(knn_f1_score,tree_f1_score))
        list(result_data)
    })
    ## Prints
    results=do.call("rbind",results)
    results%>%mutate(seed_max=seed_max,attempts_max=attempts_max,wireless=wrl,wakeupfor=wuf)

}

generate_accuracy = function(wireless,wakeupfor,steps=10, accuracy=10,ignore_hint=TRUE){
    npolicies=4
    if(ignore_hint){npolicies=npolicies-1}
    ## Generate inputs
    result=tibble()
    for(i in seq(1,160,by=steps)){ # We stop at 80% of the data (this way test set is at least 20%)
        acc=generate_accuracy_for(ignore_hint=ignore_hint,seed=i,attempts_max=accuracy,wrl=wireless,wuf=wakeupfor)
        result=rbind(result,acc)
    }
    result%>%mutate(days=seed_max*npolicies,months=days/30) # Since 3 policies (since ignore_hint=TRUE)
}

# Generate accuracy for each wireless and uptime
if(F){ # Toggle to train
accuracy=rbind(generate_accuracy("lora",60),
               generate_accuracy("lora",180),
               generate_accuracy("nbiot",60),
               generate_accuracy("nbiot",180))
}
## Summarize
result_summary=accuracy%>%group_by(wireless,wakeupfor,months,model)%>%
    summarize(
        mean_accuracy=mean(accuracy),sd_accuracy=sd(accuracy),min_accuracy=min(accuracy),max_accuracy=max(accuracy),
        mean_f1_baseline=mean(f1_baseline),sd_f1_baseline=sd(f1_baseline),min_f1_baseline=min(f1_baseline),max_f1_baseline=max(f1_baseline),
        mean_f1_hint=mean(f1_hint),sd_f1_hint=sd(f1_hint),min_f1_hint=min(f1_hint),max_f1_hint=max(f1_hint),
        mean_f1_extended=mean(f1_extended),sd_f1_extended=sd(f1_extended),min_f1_extended=min(f1_extended),max_f1_extended=max(f1_extended),
        mean_f1_hintandextended=mean(f1_hintandextended),sd_f1_hintandextended=sd(f1_hintandextended),min_f1_hintandextended=min(f1_hintandextended),max_f1_hintandextended=max(f1_hintandextended))

## Result max
metrics_peak=result_summary%>%group_by(wireless,wakeupfor,model)%>%
    summarize(max_accuracy=max(mean_accuracy))


ggplot(data=result_summary%>%mutate(model=ifelse(model=="knn","KNN","DT")),aes(linetype=model))+
    geom_line(aes(months,mean_f1_baseline,color="Baseline"),size=1.2)+
    geom_line(aes(months,mean_f1_extended,color="Extended"),size=1.2)+
    geom_line(aes(months,mean_f1_hintandextended,color="Hintandextended"),size=1.2)+labs(color="Classes",linetype="Model")+
    facet_wrap(~wireless+wakeupfor)+
    scale_x_continuous(breaks = seq(0, 15, by = 1))+
    scale_y_continuous(breaks = seq(0, 1, by = 0.1))+
    theme_bw()+
    theme(panel.grid.minor = element_blank(),
          legend.position = c(0.9,0.74),
          legend.margin = margin(2,4,2,4),
          legend.spacing=unit(-0.2,"cm"),
          legend.box.margin=margin(1,1,1,1),
          legend.box.background = element_rect(fill = "white", color = "black",size=0.8))+
    xlab("Training months")+ylab("Classes F1-Score")+
    scale_color_viridis(discrete=TRUE,end=0.7)
ggsave("figures/months_f1-score.pdf",width=8.5,height=6)
stopifnot(1)

## Plot Merge Accuracy
ggplot(data=result_summary%>%mutate(model=ifelse(model=="knn","KNN","DT")),aes(months,mean_accuracy,color=model,shape=model))+
    geom_point(size=3)+geom_line(size=1.2)+xlab("Training months")+ylab("Model accuracy")+labs(color="Models",shape="Models")+
    scale_x_continuous(breaks = seq(0, 15, by = 1))+
    scale_y_continuous(breaks = seq(0, 1, by = 0.1))+
    facet_wrap(~wireless+wakeupfor)+
    theme_bw()+
    theme(panel.grid.minor = element_blank(),
          legend.position = c(0.38,0.68),
          legend.background = element_rect(fill = "white", color = "black",size=0.8))+
    scale_color_viridis(discrete=TRUE,end=0.7)
ggsave("figures/months_accuracy.pdf",width=8.5,height=6)


## Plot accuracy + F1-Score
sapply(c("knn","tree"),function(grp){
    data=result_summary%>%filter(model==grp)
    plot=ggplot(data,aes(months,mean_accuracy))+
        geom_ribbon(aes(ymin=mean_accuracy-sd_accuracy,ymax=mean_accuracy+sd_accuracy),alpha=0.2,color=NA)+
        geom_line(size=1.1)+geom_point(size=3)+
        geom_point(data=data%>%drop_na(mean_f1_baseline),aes(months,mean_f1_baseline,color="baseline"))+geom_line(data=data%>%drop_na(mean_f1_baseline),aes(months,mean_f1_baseline,color="baseline"))+
        geom_point(data=data%>%drop_na(mean_f1_extended),aes(months,mean_f1_extended,color="extended"))+geom_line(data=data%>%drop_na(mean_f1_extended),aes(months,mean_f1_extended,color="extended"))+
        geom_point(data=data%>%drop_na(mean_f1_hintandextended),aes(months,mean_f1_hintandextended,color="hintandextended"))+geom_line(data=data%>%drop_na(mean_f1_hintandextended),aes(months,mean_f1_hintandextended,color="hintandextended"))

    if(any(!is.na(data$mean_f1_hint))){
        plot=plot+geom_point(data=data%>%drop_na(mean_f1_hint),aes(months,mean_f1_hint,color="hint"))+geom_line(data=data%>%drop_na(mean_f1_hint),aes(months,mean_f1_hint,color="hint"))
    }

    plot=plot+xlab("Number of months")+ylab(paste("Mean",grp,"accuracy"))+ggtitle(paste(grp,"accuracy"))+
#        geom_hline(data=metrics_peak%>%filter(model==grp),aes(yintercept=max_accuracy),color="red",size=1)+
#        geom_text(data=metrics_peak%>%filter(model==grp),x=0,aes(y=max_accuracy,label = round(max_accuracy,digits=1),vjust=-1),color="red")+
        facet_wrap(~wireless+wakeupfor)+scale_x_continuous(breaks = seq(0, max(result_summary$months), by = 1))+ylim(c(NA,1))+labs(color="F1 Score")+
    theme(legend.position="top")
    ggsave(paste0("figures/months_",grp,".pdf"),width=9,height=8)
    print(plot)
})