commit f877df31ca13b4b7001976a1e7c783261983c5f2
parent 5d2c87e2a4c55fb981767c4e6715bfd1597b38fe
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Sun, 18 Apr 2021 00:54:33 +0300
Plot speed RMSE.
Diffstat:
| R/common.R | | | 36 | ++++++++++++++++++++++++------------ |
| R/hold-peak.R | | | 109 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
2 files changed, 133 insertions(+), 12 deletions(-)
diff --git a/R/common.R b/R/common.R
@@ -59,18 +59,30 @@ sample_to_speed <- function(x, c1, c2) {
speed <- function(x,y,z=0) { sqrt(x**2 + y**2 + z**2) }
direction <- function(x,y) { atan2(y,x) }
-select_samples <- function (timestamp, time_delta) {
- db <- dbConnect(RSQLite::SQLite(), "samples/load-cell.sqlite3")
- velocity = dbGetQuery(db,
- "SELECT x,y,z FROM samples WHERE timestamp BETWEEN :t0 AND (:t0 + :dt) ORDER BY timestamp, t",
- params = list(t0 = timestamp, dt=time_delta))
- velocity = dbGetQuery(db,
- "SELECT x,y,z FROM samples
- WHERE timestamp BETWEEN :t0 AND (:t0 + :dt)
- AND ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000
- ORDER BY timestamp, t",
- params = list(t0 = timestamp, dt=time_delta,
- xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z)))
+select_samples <- function (timestamp, time_delta, url="samples/load-cell.sqlite3",
+ table="samples") {
+ db <- dbConnect(RSQLite::SQLite(), url)
+ velocity <- NULL
+ if (is.null(timestamp)) {
+ velocity = dbGetQuery(db,
+ sprintf("SELECT x,y,z FROM %s ORDER BY timestamp, t", table))
+ velocity = dbGetQuery(db,
+ sprintf("SELECT timestamp,x,y,z FROM %s
+ WHERE ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000
+ ORDER BY timestamp, t", table),
+ params = list(xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z)))
+ } else {
+ velocity = dbGetQuery(db,
+ sprintf("SELECT x,y,z FROM %s WHERE timestamp BETWEEN :t0 AND (:t0 + :dt) ORDER BY timestamp, t", table),
+ params = list(t0 = timestamp, dt=time_delta))
+ velocity = dbGetQuery(db,
+ sprintf("SELECT timestamp,x,y,z FROM %s
+ WHERE timestamp BETWEEN :t0 AND (:t0 + :dt)
+ AND ABS(x-:xmean)<10000 AND ABS(y-:ymean)<10000 AND ABS(z-:zmean)<10000
+ ORDER BY timestamp, t", table),
+ params = list(t0 = timestamp, dt=time_delta,
+ xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z)))
+ }
if (nrow(velocity) == 0) { return(data.frame()) }
# convert sensor values to wind speed in m/s
coefficients <- calibration_coefficients()
diff --git a/R/hold-peak.R b/R/hold-peak.R
@@ -0,0 +1,109 @@
+source("R/common.R")
+
+select_intervals <- function (table="hold_peak_samples",diff=1) {
+ db <- dbConnect(RSQLite::SQLite(), "samples/hold-peak.sqlite3")
+ tmp <- dbGetQuery(db,
+ sprintf("WITH series(t,t0,t1) AS (
+ SELECT timestamp AS t,
+ LAG(timestamp,1,NULL) OVER win AS t0,
+ LEAD(timestamp,1,NULL) OVER win AS t1
+ FROM %s
+ WINDOW win AS (ORDER BY timestamp)
+ ORDER BY timestamp)
+ SELECT t0 AS end, t AS start FROM series WHERE t0+:diff<t OR t0 IS NULL
+ UNION
+ SELECT t,t1 FROM series WHERE t+:diff<t1 OR t1 IS NULL
+ ORDER BY t0", table), params=list(diff=diff))
+ dbDisconnect(db)
+ n <- length(tmp$start)
+ data.frame(start=tmp$start[c(1:(n-1))], end=tmp$end[c(2:n)])
+}
+
+select_hold_peak_samples <- function (timestamp, time_delta) {
+ db <- dbConnect(RSQLite::SQLite(), "samples/hold-peak.sqlite3")
+ tmp <- dbGetQuery(db,
+ "SELECT timestamp,speed
+ FROM hold_peak_samples
+ WHERE timestamp BETWEEN :t0 AND (:t0+:dt)
+ ORDER BY timestamp",
+ params = list(t0 = timestamp, dt=time_delta))
+ dbDisconnect(db)
+ tmp
+}
+
+# normalized rmse
+normalized_rmse <- function(estimated, actual) {
+ s_max <- max(actual)
+ s_min <- min(actual)
+ sqrt(mean((estimated - actual)**2)) / (s_max - s_min)
+}
+
+rmse <- function(estimated, actual) {
+ sqrt(mean((estimated - actual)**2))
+}
+
+hp_intervals <- select_intervals()
+lc_intervals <- select_intervals(table="load_cell_samples",diff=2)
+#print(hp_intervals)
+#print(lc_intervals)
+# select all data to remove the trend
+
+result <- data.frame(rmse=numeric(0),start=numeric(0),end=numeric(0),urmse=numeric(0))
+for (i in 1:nrow(hp_intervals)) {
+ t0 <- hp_intervals[i,"start"]
+ t1 <- hp_intervals[i,"end"]
+ tmp <- lc_intervals[lc_intervals$start <= t0 & t0 <= lc_intervals$end,]
+ #print(tmp)
+ lc_t0 <- tmp$start
+ lc_t1 <- tmp$end
+ all_lc_velocity <- select_samples(lc_t0, lc_t1-lc_t0,
+ url="samples/hold-peak.sqlite3",
+ table="load_cell_samples")
+ lc_velocity <- all_lc_velocity[all_lc_velocity$timestamp>=t0 & all_lc_velocity$timestamp<=t1,]
+ hp_velocity <- select_hold_peak_samples(t0, t1-t0)
+ velocity <- merge(lc_velocity, hp_velocity, by=c("timestamp"))
+ # scale samples
+ t <- c(1:length(velocity$speed.x))
+ reg <- lm(velocity$speed.x-velocity$speed.y~t)
+ #print(coef(reg))
+ new_row <- nrow(result)+1
+ result[new_row,"rmse"] <- normalized_rmse(velocity$speed.x, velocity$speed.y)
+ result[new_row,"urmse"] <- rmse(velocity$speed.x, velocity$speed.y)
+ result[new_row,"start"] <- t0
+ result[new_row,"end"] <- t1
+ velocity$speed.x <- velocity$speed.x - reg$fitted.values
+ #print(lc_velocity)
+ t_origin <- max(min(velocity$timestamp),min(velocity$timestamp))
+ plot(velocity$timestamp-t_origin, velocity$speed.x)
+ lines(velocity$timestamp-t_origin, velocity$speed.y)
+ dir.create(file.path("build","hold-peak"), recursive=TRUE, showWarnings=FALSE)
+ write.table(data.frame(timstamp=velocity$timestamp,
+ speed_lc=velocity$speed.x,
+ speed_hp=velocity$speed.y),
+ file.path("build","hold-peak",sprintf("%d",t0)),
+ row.names=FALSE, quote=FALSE)
+}
+
+print(result[result$rmse==max(result$rmse),])
+print(result[result$rmse==min(result$rmse),])
+sorted_result <- result[order(result$rmse),]
+
+gnuplot <- function (file, str, ...) {
+ cat(sprintf(str, ...), file=file, append=TRUE)
+}
+
+gnuplot_all <- function (filename, prefix, rmse_column, row_min, row_max, row_longest) {
+ f <- filename
+ cat("",file=f)
+ gnuplot(f, paste(prefix,"_rmse_min = %f\n",sep=""), row_min[rmse_column])
+ gnuplot(f, paste(prefix,"_rmse_min_timestamp = %d\n",sep=""), as.integer(row_min["start"]))
+ gnuplot(f, paste(prefix,"_rmse_max = %f\n",sep=""), row_max[rmse_column])
+ gnuplot(f, paste(prefix,"_rmse_max_timestamp = %d\n",sep=""), as.integer(row_max["start"]))
+ gnuplot(f, paste(prefix,"_rmse_longest = %f\n",sep=""), row_longest[rmse_column])
+ gnuplot(f, paste(prefix,"_rmse_longest_timestamp = %d\n",sep=""), as.integer(row_longest["start"]))
+}
+
+gnuplot_all("build/gnuplot/hold_peak_rmse.gnuplot", "speed", "rmse",
+ result[result$rmse == min(result$rmse),],
+ result[result$rmse == max(result$rmse),],
+ sorted_result[2,])
