hold-peak.R (4825B)
1 source("R/common.R") 2 3 select_intervals <- function (table="hold_peak_samples",diff=1) { 4 db <- dbConnect(RSQLite::SQLite(), "samples/hold-peak.sqlite3") 5 tmp <- dbGetQuery(db, 6 sprintf("WITH series(t,t0,t1) AS ( 7 SELECT timestamp AS t, 8 LAG(timestamp,1,NULL) OVER win AS t0, 9 LEAD(timestamp,1,NULL) OVER win AS t1 10 FROM %s 11 WINDOW win AS (ORDER BY timestamp) 12 ORDER BY timestamp) 13 SELECT t0 AS end, t AS start FROM series WHERE t0+:diff<t OR t0 IS NULL 14 UNION 15 SELECT t,t1 FROM series WHERE t+:diff<t1 OR t1 IS NULL 16 ORDER BY t0", table), params=list(diff=diff)) 17 dbDisconnect(db) 18 n <- length(tmp$start) 19 data.frame(start=tmp$start[c(1:(n-1))], end=tmp$end[c(2:n)]) 20 } 21 22 select_hold_peak_samples <- function (timestamp, time_delta) { 23 db <- dbConnect(RSQLite::SQLite(), "samples/hold-peak.sqlite3") 24 tmp <- dbGetQuery(db, 25 "SELECT timestamp,speed 26 FROM hold_peak_samples 27 WHERE timestamp BETWEEN :t0 AND (:t0+:dt) 28 ORDER BY timestamp", 29 params = list(t0 = timestamp, dt=time_delta)) 30 dbDisconnect(db) 31 tmp 32 } 33 34 # normalized rmse 35 normalized_rmse <- function(estimated, actual) { 36 s_max <- max(actual) 37 s_min <- min(actual) 38 sqrt(mean((estimated - actual)**2)) / (s_max - s_min) 39 } 40 41 rmse <- function(estimated, actual) { 42 sqrt(mean((estimated - actual)**2)) 43 } 44 45 hp_intervals <- select_intervals() 46 lc_intervals <- select_intervals(table="load_cell_samples",diff=2) 47 #print(hp_intervals) 48 #print(lc_intervals) 49 # select all data to remove the trend 50 51 result <- data.frame(rmse=numeric(0),start=numeric(0),end=numeric(0),urmse=numeric(0)) 52 for (i in 1:nrow(hp_intervals)) { 53 t0 <- hp_intervals[i,"start"] 54 t1 <- hp_intervals[i,"end"] 55 tmp <- lc_intervals[lc_intervals$start <= t0 & t0 <= lc_intervals$end,] 56 #print(tmp) 57 lc_t0 <- tmp$start 58 lc_t1 <- tmp$end 59 all_lc_velocity <- select_samples(lc_t0, lc_t1-lc_t0, 60 url="samples/hold-peak.sqlite3", 61 table="load_cell_samples") 62 lc_velocity <- all_lc_velocity[all_lc_velocity$timestamp>=t0 & all_lc_velocity$timestamp<=t1,] 63 hp_velocity <- select_hold_peak_samples(t0, t1-t0) 64 velocity <- merge(lc_velocity, hp_velocity, by=c("timestamp")) 65 velocity$timestamp = velocity$timestamp - min(velocity$timestamp) 66 # scale samples 67 t <- c(1:length(velocity$speed.x)) 68 model <- nls(speed.y ~ a*speed.x, data=velocity) 69 velocity$speed.x <- velocity$speed.x * summary(model)$coefficients[[1]] 70 lo <- loess.smooth(velocity$timestamp, velocity$speed.x, span=0.5) 71 new_row <- nrow(result)+1 72 result[new_row,"rmse"] <- normalized_rmse(velocity$speed.x, velocity$speed.y) 73 result[new_row,"urmse"] <- rmse(velocity$speed.x, velocity$speed.y) 74 result[new_row,"start"] <- t0 75 result[new_row,"end"] <- t1 76 scatter.smooth(velocity$timestamp, velocity$speed.x) 77 #plot(velocity$timestamp, velocity$speed.x) 78 lines(velocity$timestamp, velocity$speed.y, col='red') 79 #lines(lo) 80 make_directory(file.path("build","hold-peak")) 81 filename <- file.path("build","hold-peak",sprintf("%d",t0)) 82 write.table(data.frame(timestamp=velocity$timestamp, 83 speed_lc=velocity$speed.x, 84 speed_hp=velocity$speed.y), 85 filename, row.names=FALSE, quote=FALSE) 86 write("\n",file=filename,append=TRUE) 87 write.table(data.frame(timestamp=lo$x,speed=lo$y), filename, 88 row.names=FALSE, quote=FALSE, append=TRUE) 89 } 90 91 print(result[result$rmse==max(result$rmse),]) 92 print(result[result$rmse==min(result$rmse),]) 93 sorted_result <- result[order(result$rmse),] 94 print(sorted_result) 95 96 gnuplot <- function (file, str, ...) { 97 cat(sprintf(str, ...), file=file, append=TRUE) 98 } 99 100 gnuplot_all <- function (filename, prefix, rmse_column, row_min, row_max, row_longest) { 101 f <- filename 102 cat("",file=f) 103 gnuplot(f, paste(prefix,"_rmse_min = %f\n",sep=""), row_min[rmse_column]) 104 gnuplot(f, paste(prefix,"_rmse_min_timestamp = %d\n",sep=""), as.integer(row_min["start"])) 105 gnuplot(f, paste(prefix,"_rmse_max = %f\n",sep=""), row_max[rmse_column]) 106 gnuplot(f, paste(prefix,"_rmse_max_timestamp = %d\n",sep=""), as.integer(row_max["start"])) 107 gnuplot(f, paste(prefix,"_rmse_longest = %f\n",sep=""), row_longest[rmse_column]) 108 gnuplot(f, paste(prefix,"_rmse_longest_timestamp = %d\n",sep=""), as.integer(row_longest["start"])) 109 } 110 111 gnuplot_all("build/gnuplot/hold_peak_rmse.gnuplot", "speed", "rmse", 112 result[result$rmse == min(result$rmse),], 113 result[result$rmse == max(result$rmse),], 114 result[2,])