library("tidyverse")
library("class")
library("rpart")
library("rpart.plot")
library("viridis")

## 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=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
tmp_data_coverage=data%>%group_by(simkey,wireless,wakeupfor,seed)%>%mutate(coverage=sum(nDataRcv))%>%ungroup()%>%filter(isSender==1)%>%select(simkey,wireless,wakeupfor,seed,coverage)
data_seed_isSender=data%>%group_by(simkey,wireless,wakeupfor,seed,isSender)%>%summarize(energy_mean=mean(energy))%>%
    left_join(tmp_data_coverage,by=c("simkey","wireless","wakeupfor","seed"))%>%
    mutate(efficiency=energy_mean/coverage)%>%
    ungroup()
data_seed=data%>%group_by(simkey,wireless,wakeupfor,seed)%>%summarize(energy=sum(energy),coverage=sum(nDataRcv))%>%
    mutate(efficiency=energy/coverage)%>%
    ungroup()

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)
    data_seed=data_seed%>%filter(wakeupfor==wuf,wireless==wrl)
    data_ml=data_seed%>%select(-efficiency)%>%mutate(wireless=wireless_map[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%>%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=10)
    ## KNN analysis
    knn_cont_table=table(knn_predictions,test_set$simkey)
    knn_accuracy=round((sum(diag(knn_cont_table)/sum(rowSums(knn_cont_table))))*100)
    knn_prop_table=round(prop.table(knn_cont_table),digits=2)

    ## 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=round((sum(diag(tree_cont_table)/sum(rowSums(tree_cont_table))))*100)
    tree_prop_table=round(prop.table(tree_cont_table),digits=2)
     list(tibble(seed_max=seed_max,knn_accuracy=knn_accuracy,tree_accuracy=tree_accuracy))
    })
    ## 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=20, accuracy=10){
    ## Generate inputs
    result=tibble()
    for(i in seq(1,200,by=steps)){
        acc=generate_accuracy_for(ignore_hint=TRUE,seed=i,attempts_max=accuracy,wrl=wireless,wuf=wakeupfor)
        result=rbind(result,acc)
    }
    result%>%mutate(days=seed_max*3) # Since 3 policies (since ignore_hint=TRUE)
}

## Generate accuracy for each wireless and uptime
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,seed_max)%>%summarize(attempts_max=mean(attempts_max),days=mean(days),mean_knn_accuracy=mean(knn_accuracy),sd_knn_accuracy=sd(knn_accuracy),min_knn_accuracy=min(knn_accuracy),max_knn_accuracy=max(knn_accuracy),mean_tree_accuracy=mean(tree_accuracy),sd_tree_accuracy=sd(tree_accuracy),min_tree_accuracy=min(tree_accuracy),max_tree_accuracy=max(tree_accuracy))

## Result max
result_max=result_summary%>%group_by(wireless,wakeupfor)%>%summarize(max_knn_mean=max(mean_knn_accuracy),max_tree_mean=max(mean_tree_accuracy))

## Plot
ggplot(result_summary,aes(days,mean_knn_accuracy))+
    geom_errorbar(aes(ymin=min_knn_accuracy,ymax=max_knn_accuracy),width=5)+
    geom_boxplot(aes(ymin=min_knn_accuracy,ymax=max_knn_accuracy,
                                    middle=mean_knn_accuracy,
                                    upper=mean_knn_accuracy+sd_knn_accuracy,
                                    lower=mean_knn_accuracy-sd_knn_accuracy,group=days),stat="identity",fill="grey")+
    geom_line(size=1.1)+geom_point(size=3,pch=15)+xlab("Number of training days")+ylab("Mean KNN accuracy")+ggtitle("KNN Accuracy")+
    ylim(c(0,100))+
    facet_wrap(~wireless+wakeupfor)+
    geom_hline(data=result_max,aes(yintercept=max_knn_mean),color="red",size=1)+
    geom_text(data=result_max, geom="text",x=0,aes(y=max_knn_mean,label = max_knn_mean,vjust=-1),color="red")
ggsave("figures/days_knn.pdf")

ggplot(result_summary,aes(days,mean_tree_accuracy))+
    geom_errorbar(aes(ymin=min_tree_accuracy,ymax=max_tree_accuracy),width=5)+
    geom_boxplot(aes(ymin=min_tree_accuracy,ymax=max_tree_accuracy,
                                    middle=mean_tree_accuracy,
                                    upper=mean_tree_accuracy+sd_tree_accuracy,
                                    lower=mean_tree_accuracy-sd_tree_accuracy,group=days),stat="identity",fill="grey")+
    geom_line(size=1.1)+geom_point(size=3,pch=15)+xlab("Number of training days")+ylab("Mean tree accuracy")+ggtitle("TREE Accuracy")+
    ylim(c(0,100))+
    facet_wrap(~wireless+wakeupfor)+
    geom_hline(data=result_max,aes(yintercept=max_tree_mean),color="red",size=1)+
    geom_text(data=result_max, geom="text",x=0,aes(y=max_tree_mean,label = max_tree_mean,vjust=-1),color="red")
ggsave("figures/days_tree.pdf")