##############################################################################################
#' @title Calculate standard deviation and stand error at different aggregation stage: from minutely to hourly, to mean diurnal cycle and monthly average

#' @author
#' Ke Xu \email{xuke2012abroad@gmail.com} \cr

#' @description 
#' Function defintion. Calculate standard deviation and stand error at different aggregation stage: from minutely to hourly, to mean diurnal cycle and monthly average

#' @param Currently none

#' @return Currently none

#' @references Currently none

#' @keywords standard deviation, standard error, aggregation, aggregate

#' @examples Currently none

#' @seealso Currently none

#' @export

# changelog and author contributions / copyrights
#   Ke Xu (TBD)
#     original creation
##############################################################################################



#this needs to go into the function call:
# train.data <- read.csv(paste("../kxu_NEON_140912/tower_ERF/out/Phase2-footprint-2011_JA-122-150225-DEM-w/H_en/H_en_test.data_160224.csv", sep=""))
# str(train.data) 


sdse <- function(
  DOY = train.data$DOY - 6/24,   #input DOY time in CST time
  hour = train.data$hour,   #input hour in CST time
  meanMinu = train.data$response, #input minutely mean
  sdMinu = train.data$sd   #input minutely standard deviation
  
  ){
  
  
  data <- list()
  data$DOY <- DOY
  data$hour <- hour
  data$meanMinu <- meanMinu
  
  data$varMinu <- sdMinu
  data$seMinu_2 <- data$varMinu/(freq * window_size * 60)
  data$seMinu <- sqrt(data$seMinu_2)
  

  ################  hourly averaged  #####################
 
  dates <- sort(unique(floor(data$DOY)))
  
  data$date <- floor(data$DOY)
  
  
  data_hourly <- list()
  flag <- 0
  for (dd in unique(data$date)){
    # dd <- unique(data$date)[1]
    for(hh in 0:23){
      #hh <- 18
      if(length(which(floor(data$hour)%%24 == hh & data$date == dd)) == 0){
        tmp <- NaN
        tmp_var <- NaN
        tmp_seHour_2 <- NaN
      } else {
        tmp <-  mean(data$meanMinu[which(floor(data$hour) == hh & data$date == dd)], na.rm = TRUE)
        tmp_var <- stats::var(data$meanMinu[which(floor(data$hour) == hh & data$date == dd)], na.rm = TRUE) + mean(data$varMinu[which(floor(data$hour) == hh & data$date == dd)], na.rm=T)
        tmp_seHour_2 <- (stats::var(data$meanMinu[which(floor(data$hour) == hh & data$date == dd)], na.rm = TRUE) + mean(data$seMinu_2[which(floor(data$hour) == hh & data$date == dd)])) / 54
      }
      
      flag <- flag + 1
      tmp_1 <- data.frame(dd, hh, tmp, tmp_var, tmp_seHour_2)
      
      
      #dimnames(eper_out)[[1]][1] <- whr_indi_minu
      dimnames(tmp_1)[[2]] <- c("day", "hour", "meanHour" , "varHour", "seHour_2")
      
      if(flag == 1){
        data_hourly <- tmp_1  
      } else {
        data_hourly <- rbind(data_hourly, tmp_1)
      }
      
    }
  }
  
  
  data_hourly$seHour <- sqrt(data_hourly$seHour_2)
  
  
  #############  mean diurnal cycle  #####################
  data_diu <- list()
  #mean
  data_diu$mn <- sapply(0:23, function(tt)  mean(data_hourly$meanHour[which(floor(data_hourly$hour)%%24 == tt)], na.rm = TRUE))
  #data_diu$MED <- sapply(0:23, function(tt)  MEDmad(data_hourly$meanHour[which(floor(data_hourly$hour)%%24 == tt)])[1])
  
  
  #sd
  ##To calculate data_diu$varMdc = data_diu$var_meanHour + data_diu$mean_varHour
  data_diu$varMdc <- sapply(0:23, function(tt)  (stats::var(data_hourly$meanHour[which(floor(data_hourly$hour)%%24 == tt)], na.rm = TRUE) +
                                                  mean(data_hourly$varHour[which(floor(data_hourly$hour)%%24 == tt)], na.rm = TRUE))
  )
  
  data_diu$num <- sapply(0:23, function(tt) length(which(floor(data_hourly$hour)%%24 == tt & !is.na(data_hourly$meanHour))))
  
  ##ensemble error:
  data_diu$seMdc_2 <- sapply(0:23, function(tt) {(stats::var(data_hourly$meanHour[which(floor(data_hourly$hour)%%24 == tt)], na.rm = TRUE) + mean(data_hourly$seHour_2[which(floor(data_hourly$hour)%%24 == tt)], na.rm = TRUE))/data_diu$num[tt+1]}
  )
  
  data_diu$seMdc <- sqrt(data_diu$seMdc_2)
  
  
  
  ###############  monthly averaged ######################
  Mon <- list()

  Mon$meanMon <- sum(data_diu$mn, na.rm=T)/24
  Mon$varMon <- mean(data_diu$varMdc, na.rm=T) + stats::var(data_diu$mn, na.rm=T)
  Mon$sdMon <- sqrt(Mon$varMon)
  Mon$seMon_2 <- (mean(data_diu$seMdc_2, na.rm=T) + stats::var(data_diu$mn, na.rm=T))/length(which(!is.na(data_diu$mn)))
  Mon$seMon <- sqrt(Mon$seMon_2)
  
  return(Mon)
}

#this needs to go into the function call:
#sdse()