########## INFORMATIONS ##########
# This file is made to study online classification
# So, each pair (wireless,wakeupfor) has its classification models (knn and decision tree)
# Note that is the following error appears: object 'accuracy' not found
# it means you should toggle the boolean in the if condition in the code to generate the accuracy onbject
##################################

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 and prepare the 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_Score of multiclass vectors
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)
}

## Down scale data
reduce_days=function(data,every=6){data%>%filter(((days/3) %% every) == 0)}

## GGPlot Theme
th = function(option="D") {list(theme_bw(),
                                theme(legend.box.background = element_rect(fill = "white", color = "black",size=0.9),
                                      strip.background =element_rect(fill="#F5F5F5")),
                     scale_color_viridis(discrete=TRUE,option=option,end=0.95),
                     scale_fill_viridis(discrete=TRUE,option=option,end=0.95))}

## Train models
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_energy = function(wireless,wakeupfor,steps=1, accuracy=10,ignore_hint=TRUE){
    ## Setup variables
    npolicies=4
    data_seed_for_energy=data_seed%>%ungroup()
    if(ignore_hint){npolicies=npolicies-1;data_seed_for_energy=data_seed%>%filter(simkey!="hint")}
    data_seed_for_energy=data_seed_for_energy%>%filter(wireless==!!wireless,wakeupfor==!!wakeupfor)

    ## Generate inputs
    result=tibble()
    result_energy=tibble()
    for(i in seq(1,160,by=steps)){ # We stop at 80% of the data (this way test set is at least 20%)
        print(paste("Step",i))
        acc=generate_accuracy_for(ignore_hint=ignore_hint,seed=i,attempts_max=accuracy,wrl=wireless,wuf=wakeupfor)
        result_energy=rbind(result_energy,data_seed_for_energy%>%filter(seed<=i)%>%summarize(energy=sum(energy),seed_max=i,days=i*npolicies,setup=paste0(!!wireless," ",!!wakeupfor,"s"),wireless=!!wireless,wakeupfor=!!wakeupfor))
        result=rbind(result,acc)
    }
   list(accuracy=result%>%mutate(days=seed_max*npolicies), # Since 3 policies (since ignore_hint=TRUE)
        energy=result_energy)
}


########## Generate accuracy, energy, F1-Score and coverage data ##########
if(F){ # Toggle to train
    lora60=generate_accuracy_energy("lora",60)
    lora180=generate_accuracy_energy("lora",180)
    nbiot60=generate_accuracy_energy("nbiot",60)
    nbiot180=generate_accuracy_energy("nbiot",180)
    accuracy=rbind(lora60$accuracy,lora180$accuracy,nbiot60$accuracy,nbiot180$accuracy)
    coverage=data_seed%>%filter(simkey!="hint")%>%group_by(wireless,wakeupfor,seed)%>%summarize(coverage=sum(coverage))%>%mutate(days=seed*3)%>%filter(days %in% !!lora60$energy$days)%>%select(-seed)
    energy=rbind(lora60$energy,
                 lora180$energy,
                 nbiot60$energy,
                 nbiot180$energy)%>%left_join(coverage,by=c("wireless","wakeupfor","days"))
}


########## Build learning curve (Accuracy+F1_Scores) ##########
learning_curves=accuracy%>%group_by(wireless,wakeupfor,days,model)%>%
    summarize(
        ## Accuracy
        mean_accuracy=mean(accuracy),sd_accuracy=sd(accuracy),
        min_accuracy=min(accuracy),max_accuracy=max(accuracy),
        ## F1-Score Baseline
        mean_f1_baseline=mean(f1_baseline),sd_f1_baseline=sd(f1_baseline),
        min_f1_baseline=min(f1_baseline),max_f1_baseline=max(f1_baseline),
        ## F1-Score Hint
        mean_f1_hint=mean(f1_hint),sd_f1_hint=sd(f1_hint),min_f1_hint=min(f1_hint),
        max_f1_hint=max(f1_hint),
        ## F1-Score Extended
        mean_f1_extended=mean(f1_extended),sd_f1_extended=sd(f1_extended),
        min_f1_extended=min(f1_extended),max_f1_extended=max(f1_extended),
        ## F1-Score HintAndExtended
        mean_f1_hintandextended=mean(f1_hintandextended),sd_f1_hintandextended=sd(f1_hintandextended),
        min_f1_hintandextended=min(f1_hintandextended),max_f1_hintandextended=max(f1_hintandextended))


ggplot(data=learning_curves%>%mutate(model=ifelse(model=="knn","KNN","DT"))%>%reduce_days(6),aes(linetype=model))+
    geom_line(aes(days/30,mean_f1_baseline,color="Baseline"),size=1.2)+
    geom_line(aes(days/30,mean_f1_extended,color="Extended"),size=1.2)+
    geom_line(aes(days/30,mean_f1_hintandextended,color="Hintandextended"),size=1.2)+labs(color="Classes",linetype="Model")+
    facet_wrap(~wireless+wakeupfor)+
    scale_x_continuous(breaks = seq(0, max(learning_curves$days/30)))+
    scale_y_continuous(breaks = seq(0, 1, by = 0.1))+
    th()+
    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))+
    xlab("Training duration (months)")+ylab("Classes F1-Score")
ggsave("figures/days_f1-score.pdf",width=8.5,height=6)


## Plot Merge Accuracy
ggplot(data=learning_curves%>%mutate(model=ifelse(model=="knn","KNN","DT"))%>%reduce_days(3),
       aes(days/30,mean_accuracy))+
    geom_line(aes(linetype=model),size=1.2)+xlab("Training duration (months)")+ylab("Model accuracy")+labs(linetype="Models")+
    scale_x_continuous(breaks = seq(0, max(learning_curves$days/30)))+
    scale_y_continuous(breaks = seq(0, 1, by = 0.1))+
    facet_wrap(~wireless+wakeupfor)+
    th()+
    theme(panel.grid.minor = element_blank(),
          legend.position = c(0.38,0.68))
ggsave("figures/days_accuracy.pdf",width=8.5,height=6)


########## Energy and delta with raw policies ##########
## First we extended the number of seed to cover the entire duration of the training
data_seed_energy=rbind(data_seed,
                           data_seed%>%mutate(seed=seed+200),
                           data_seed%>%mutate(seed=seed+400),
                           data_seed%>%mutate(seed=seed+600)) # Almost same as if each experiment run 4 times more seed seed
## Compute the cumulative energy of each policies in each configuration accross the 4*200 days
data_seed_energy=data_seed_energy%>%group_by(wireless,wakeupfor,simkey)%>%mutate(energy=cumsum(energy),setup=paste0(wireless," ",wakeupfor,"s"),days=seed)
## Now filter the data
data_seed_energy=data_seed_energy%>%filter(days %in% !!energy$days)
## Compute the delta data
energy_coverage_delta=data_seed_energy%>%
    full_join(energy,by=c("days","wireless","wakeupfor"),suffix=c("","_training"))
energy_coverage_delta=energy_coverage_delta%>%group_by(wireless,wakeupfor)%>%summarize(delta_energy=energy-energy_training,simkey=simkey,days=days,delta_coverage=coverage-coverage_training)


write("wireless,wakeupfor,policy,slope,intercept,delta_coverage","figures/delta_energy_coverage.csv")
energy_coverage_delta%>%group_by(wireless,wakeupfor,simkey)%>%group_walk(function(data,grp){
    grp=as.list(grp)
    reg=lm(delta_energy ~ days,data)
    slope=round(as.numeric(reg$coefficients["days"]),digits=1)
    intercept=round(as.numeric(reg$coefficients[1]),digits=1)
    mean_delta_coverage=round(mean(data$delta_coverage),digits=1)
    print(paste0("Wireless=",grp$wireless," Wakeupfor=",grp$wakeupfor," Policy=",grp$simkey," Slope=",slope," Intercept=",intercept," Delta Coverage=",mean_delta_coverage))
    write(paste(grp$wireless,grp$wakeupfor,grp$simkey,slope,intercept,mean_delta_coverage,sep=","),"figures/delta_energy_coverage.csv",append=T)
})

ggplot(energy_coverage_delta,aes(days/30,delta_energy/1e3,color=simkey,shape=simkey))+
    geom_line(size=1.2)+ylab("Delta in energy (kJ)")+xlab("Training duration (months)")+
    facet_wrap(~wireless+wakeupfor,scale="free")+
    th()+theme(legend.position=c(0.61,0.93))+labs(color="Classes")
ggsave("figures/delta_energy_training.pdf",height=6,width=10)

ggplot(energy_coverage_delta,aes(days/30,delta_coverage,color=simkey))+
    geom_line(size=1.2)+ylab("Delta in coverage")+xlab("Training duration (months)")+
    facet_wrap(~wireless+wakeupfor,scale="free")+
    th()+theme(legend.position="top")+labs(color="Classes")
ggsave("figures/delta_coverage_training.pdf",width=9)

ggplot(data=energy%>%reduce_days(6),aes(days/30,energy/1e6,group=setup,fill=setup))+
    geom_bar(stat="identity",position="dodge")+
    labs(fill="Wireless and Uptime")+
    scale_x_continuous(breaks = seq(0, max(energy$days/30)))+
    xlab("Training duration (months)")+ylab("Energy consumption (MJ)")+
    th()+theme(legend.position=c(0.12,0.75))
ggsave("figures/days_energy.pdf",width=8.5,height=4)