library(ggplot2)
source("defines.R")
source("gaussfilt.R")

# Cyclus Measurements
cyclus <- read.csv("cyclusintervals.csv")
# Mechanisches Gangverhältnis
gear_ratio <- cyclus$ErgoFrontGear[1]/cyclus$ErgoBackGear[1]

# SRM Measurements
# 0: Cadence
# 1: Torque
srm <- read.csv("srmintervals.csv")
srm_cadence <- plyr::rename(subset(srm, Type==0, select=c("Time", "Value")), c("Value"="Cadence"))
srm_torque <- plyr::rename(subset(srm, Type==1, select=c("Time", "Value")), c("Value"="Torque"))

# Arduino Measurements
# 0: Photo Sensor (back, 25 teeth, Value: time difference)
# 1: Reed Switch (front, Value: time difference)
# 2: Reed Switch (back, Value: time difference)
arduino <- read.csv("arduinointervals.csv")
arduino_photoBack <- plyr::rename(subset(arduino, Type==0, select=c("Time", "Value")), c("Value"="TimeDiff"))
arduino_photoBack$RPS <- 1/25 / arduino_photoBack$TimeDiff
arduino_photoBack$Cadence <- arduino_photoBack$RPS / gear_ratio * 60

# Cyclus2 and SRM power
# with SRM Cadence
srm_torque$PowerSRMCad <- srm_torque$Torque * 2 * pi * approx(x=srm_cadence$Time, y=srm_cadence$Cadence, xout=srm_torque$Time, yleft = 0, yright = 0)[[2]] / 60
# with Arduino photo cadence
srm_torque$PowerArduinoCad <- srm_torque$Torque * 2 * pi * approx(x=arduino_photoBack$Time, y=arduino_photoBack$Cadence, xout=srm_torque$Time, yleft = 0, yright = 0)[[2]] / 60
# filtered
srm_torque$PowerSRMCadFiltered <- gaussfilt(srm_torque$Time, srm_torque$PowerSRMCad, 0.25)

xmin <- 195
xmax <- 205

ggplot() +
  geom_line(data=srm_torque, aes(x=Time, y=PowerSRMCad, color="SRM"), size = 1) + 
  geom_line(data=srm_torque, aes(x=Time, y=PowerArduinoCad, color="SRM with Arduino"), size = 1, linetype = "dashed") + 
  geom_line(data=cyclus, aes(x=Time, y=Power, color="Cyclus2"), size = 2) + 
  geom_line(data=srm_torque, aes(x=Time, y=PowerSRMCadFiltered, color="SRM (filtered)"), size = 2) + 
  scale_colour_manual(name = "Device", 
                      values = c("Cyclus2" = "blue", "SRM with Arduino"="yellow", "SRM"="red", "SRM (filtered)"="black"),
                      labels = c("Cyclus2"=paste0("Cyclus2 (mean=",round(mean(cyclus$Power),2),")"),
                                 "SRM"=paste0("SRM (mean=",round(mean(srm_torque$PowerSRMCad),1),")"),
                                 "SRM with Arduino"=paste0("SRM with Arduino (mean=",round(mean(srm_torque$PowerArduinoCad),1),")"),
                                 "SRM (filtered)"=paste0("SRM filtered (mean=",round(mean(srm_torque$PowerSRMCadFiltered),1),")"))) +
  plot_style +# coord_cartesian(xlim = c(xmin, xmax), ylim = c(0, 500)) +
  labs(title = "Power")

# Compute error between SRM and Cyclus2 power

# Error functions
# Root-Mean-Square Error
rmse <- function(x1, x2)
{
  return(sqrt(mean((x1-x2)^2)))
}
# Mean Absolute Percentage Error
mape <- function(x1, x2)
{
  return(100*mean(abs((x1-x2)/x1)))
}

# eliminate zeros for mape
cyclus0 <- subset(cyclus, Power!=0)
srm0 <- subset(srm_torque, PowerSRMCad!=0)

srmpower <- approx(x=srm0$Time, y=srm0$PowerSRMCad, xout=cyclus0$Time, yleft = 0, yright = 0)[[2]]
srmpowerfilt <-approx(x=srm0$Time, y=srm0$PowerSRMCadFiltered, xout=cyclus0$Time, yleft = 0, yright = 0)[[2]] 
mape_cs <- mape(cyclus0$Power, srmpower)
mape_csf <- mape(cyclus0$Power, srmpowerfilt)
rmse_cs <- rmse(cyclus0$Power, srmpower)
rmse_csf <- rmse(cyclus0$Power, srmpowerfilt)

print(paste0("MAPE to unfiltered: ", round(mape_cs,2), " %, MAPE to filtered: ", round(mape_csf, 2)," %."))
print(paste0("RMSE to unfiltered: ", round(rmse_cs,2), " W, RMSE to filtered: ", round(rmse_csf, 2), " W."))