commit b01df3952464949a914cb211effff1f21bf59758
parent 71bc63d5d5d2d701c5710ea1647d5cd420d56243
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Mon, 26 Apr 2021 14:05:18 +0300
Compare VTB and anemometer.
Diffstat:
| R/anal.R | | | 307 | +------------------------------------------------------------------------------ |
| R/common.R | | | 19 | +++++++++++++++++++ |
| R/distributions.R | | | 273 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| R/vtestbed.R | | | 73 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
4 files changed, 366 insertions(+), 306 deletions(-)
diff --git a/R/anal.R b/R/anal.R
@@ -7,244 +7,7 @@ library(gplots) # hist2d
registerDoParallel(cores=detectCores())
source("R/common.R")
-
-datetime_format <- "%FT%H:%M"
-
-# https://stackoverflow.com/questions/32370485/convert-radians-to-degree-degree-to-radians
-rad2deg <- function(rad) {(rad * 180) / (pi)}
-deg2rad <- function(deg) {(deg * pi) / (180)}
-
-# https://stackoverflow.com/questions/27455948/in-r-whats-the-simplest-way-to-scale-a-vector-to-a-unit-vector
-unit <- function(x) {
- length = sqrt(sum(x**2))
- ifelse(length == 0, x, x/length)
-}
-
-# retain the original sign
-power <- function(x,y) {
- ifelse(abs(x) < 1, sign(x)*abs(x)**y, x**y)
- #sign(x)*abs(x)**y
- #abs(x)**y
-}
-
-# normalized rmse
-rmse <- function(estimated, actual) {
- s_max <- max(actual)
- s_min <- min(actual)
- sqrt(mean((estimated - actual)**2)) / (s_max - s_min)
-}
-
-dweibull2 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g=0) {
- #ifelse(x<g,
- # a1* abs(c1*b1)*(abs(b1*(x-g))**(c1-1))* exp(-(abs(b1*(x-g))**c1)),
- # a1* abs(c2*b2)*(abs(b2*(x-g))**(c2-1))* exp(-(abs(b2*(x-g))**c2)))
- ifelse(x<g,
- a1*abs(b1*c1)*exp(-(abs(b1*(abs(x-g)-abs((c1-1)/c1)**(1/c1)))**c1)),
- a2*abs(b2*c2)*exp(-(abs(b2*(abs(x-g)-abs((c2-1)/c2)**(1/c2)))**c2)))
-}
-
-dweibull4 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g1=0,
- a3=1, b3=1, c3=1, a4=1, b4=1, c4=1, g2=0) {
- ifelse(x<0,
- dweibull2(x,a1,b1,c1,a2,b2,c2,g1),
- dweibull2(x,a3,b3,c3,a4,b4,c4,g2))
-}
-
-dvonmises <- function (x, x_mean=0, a=1, b=1, c=1) {
- a * exp(b*cos(c*(x-x_mean))) / (2*pi*besselI(b,0))
-}
-
-dgumbel <- function (x, location=0, scale=1) {
- z <- (x-location)/scale
- (1/scale) * exp(-(z + exp(-z)))
-}
-
-# https://www.tutorialspoint.com/r/r_mean_median_mode.htm
-# https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode
-getmode <- function(v) {
- #uniqv <- unique(v)
- #uniqv[which.max(tabulate(match(v, uniqv)))]
- as.numeric(names(sort(-table(v)))[1])
-}
-
-# https://stackoverflow.com/questions/27418461/calculate-the-modes-in-a-multimodal-distribution-in-r
-find_peaks_indices <- function (x) {
- modes <- c()
- for (i in 2:(length(x)-1)){
- if ((x[i] > x[i-1]) & (x[i] > x[i+1])) {
- modes <- c(modes,i)
- }
- }
- modes
-}
-
-find_real_peaks_indices <- function (x, max_modes=2) {
- n <- 1024
- #x_density <- density(x,n=n,adjust=0.125)
- #y <- x_density$y
- y <- x
- modes <- find_peaks_indices(y)
- # find the first N highest peaks
- modes_xy <- data.frame(x=modes,y=y[modes])
- modes_xy <- modes_xy[order(modes_xy$y, decreasing=TRUE),]
- modes_xy <- modes_xy[c(1:max_modes), ]
- # scale indices to [0,1]
- (modes_xy$x-1)/(n-1)
-}
-
-find_max_x <- function (x, y) {
- n <- length(y)-1
- indices <- c(0:n)
- index <- indices[y == max(y)]/n
- x[floor(index*(length(x)-1)+1)]
-}
-
-normalize_hist <- function(pdf_x) {
- #data.frame(density=pdf_x$density/max(abs(pdf_x$density)),x=pdf_x$mids)
- data.frame(density=pdf_x$density,x=pdf_x$mids)
-}
-
-fit_velocity_pdf <- function (xx, x, density, sign=-1, plot=TRUE, axis="x", dist="weibull2") {
- a0 <- 0
- b0 <- 0
- c0 <- 1
- c1 <- 3
- indices <- find_real_peaks_indices(x, max_modes=1)
- if (FALSE) {
- pdf <- normalize_hist(hist(x, plot=FALSE, breaks=1000))
- print(indices)
- indices = floor(indices * length(pdf$x))
- x_modes <- pdf$x[indices]
- print(x_modes)
- model <- nls(pdf$density ~ dweibull4(pdf$x,a1,b1,c1,a2,b2,c2,x_modes[[1]],a3,b3,c3,a4,b4,c4,x_modes[[2]]),
- data=pdf,
- start=list(a1=1,b1=1,c1=1,a2=1,b2=1,c2=1,
- a3=1,b3=1,c3=1,a4=1,b4=1,c4=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(a0,b0,c0,a0,b0,c0,a0,b0,c0,a0,b0,c0),
- upper=c(1000,1000,3,1000,2000,3,1000,1000,3,1000,2000,3))
- c = summary(model)$coefficients
- if (plot) {
- pdf_x_est <- dweibull4(pdf$x,c[[1]],c[[2]],c[[3]],c[[4]],c[[5]],c[[6]],x_modes[[1]],
- c[[7]],c[[8]],c[[9]],c[[10]],c[[11]],c[[12]],x_modes[[2]])
- plot(pdf$x, pdf$density, xlab=paste('Velocity', axis), ylab='PDF')
- lines(pdf$x, pdf_x_est, col='red')
- }
- c
- } else {
- #x_mode <- getmode(x)
- #pdf_x <- normalize_hist(hist(x, plot=FALSE, breaks=100))
- pdf_x <- data.frame(x=xx,density=x)
- c <- c()
- pdf_x_set <- c()
- if (dist == "weibull2") {
- x_mode <- find_max_x(pdf_x$x,pdf_x$density)
- #indices = floor(indices * length(pdf_x$x)) + 1
- #x_mode <- pdf_x$x[indices][[1]]
- c <- tryCatch(
- {
- model <- nls(pdf_x$density ~ dweibull2(pdf_x$x,a,b,c,d,e,f,g=x_mode),
- data=pdf_x,
- start=list(a=1,b=1,c=1,d=1,e=1,f=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(a0,b0,c0,a0,b0,c0), upper=c(1000,1000,c1,1000,2000,c1))
- summary(model)$coefficients
- },
- error=function(cond) {
- message(cond)
- c(1,1,1,1,1,1)
- })
- pdf_x_est <- dweibull2(pdf_x$x,c[[1]],c[[2]],c[[3]],c[[4]],c[[5]],c[[6]],x_mode)
- } else if (dist == "vonmises") {
- c <- tryCatch({
- model <- nls(density ~ dvonmises(2*pi*pdf_x$x/(max(pdf_x$x)),0,a,b,c),
- data=pdf_x,
- start=list(a=1,b=1,c=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(0,1,0), upper=c(1000,1000,pi))
- c <- summary(model)$coefficients
- print(c)
- },error=function(cond) {
- message(cond)
- c(1,1,1)
- })
- pdf_x_est <- dvonmises(2*pi*pdf_x$x/max(pdf_x$x),0,c[[1]],c[[2]],c[[3]])
- } else if (dist == "gumbel") {
- c <- tryCatch(
- {
- model <- nls(pdf_x$density ~ c*dgumbel(pdf_x$x,location=a,scale=b),
- data=pdf_x,
- start=list(a=1,b=1,c=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(-1000,0.1,c=1), upper=c(1000,1000,1000))
- c = summary(model)$coefficients
- },
- error=function(cond) {
- message(cond)
- c <- c(1,1,1,1)
- })
- pdf_x_est <- c[[3]]*dgumbel(pdf_x$x,location=c[[1]],scale=c[[2]])
- } else if (dist == "norm") {
- fit <- fitdist(pdf_x$x, "norm")
- c <- coef(fit)
- pdf_x_est <- dnorm(pdf_x$x,mean=c[[1]],sd=c[[2]])
- } else if (dist == "rayleigh") {
- c <- tryCatch(
- {
- model <- nls(pdf_x$density ~ b*dweibull(sign*pdf_x$x,scale=a,shape=2),
- data=pdf_x,
- start=list(a=1,b=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(0.1,0), upper=c(1000,1000))
- c = summary(model)$coefficients
- },
- error=function(cond) {
- message(cond)
- c <- c(1,1,1,1)
- })
- pdf_x_est <- c[[2]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=2)
- } else {
- x_mode <- find_max_x(pdf_x$x,pdf_x$density)
- c <- tryCatch(
- {
- model <- nls(pdf_x$density ~ c*dweibull(sign*pdf_x$x,scale=a,shape=b),
- data=pdf_x,
- start=list(a=1,b=1,c=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(0.1,0.1,0),
- upper=c(1000,1000,1000))
- c = summary(model)$coefficients
- },
- error=function(cond) {
- message(cond)
- c <- c(1,1,1,1)
- })
- pdf_x_est <- c[[3]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=c[[2]])
- }
- #if (plot) {
- # plot(sign*pdf_x$x, pdf_x$density, xlab=paste('Velocity', axis), ylab='PDF')
- # lines(density(sign*x), col='blue')
- # lines(sign*pdf_x$x, pdf_x_est, col='red')
- #}
- #if (sign < 0) {
- # print(pdf_x$x)
- # pdf_x$x = -rev(pdf_x$x)
- # pdf_x$density = rev(pdf_x$density)
- # pdf_x_est = rev(pdf_x_est)
- #}
- list(coefficients=c,
- actual_density=pdf_x$density,
- estimated_density=pdf_x_est,
- x=pdf_x$x,
- residual=max(abs(pdf_x$density-pdf_x_est)))
- }
-}
+source("R/distributions.R")
fit_velocity_pdf_bilateral <- function (velocity, axis="x", timestamp=0, dist="weibull") {
velocity_hist <- normalize_hist(hist(velocity, plot=FALSE, breaks=100))
@@ -277,45 +40,6 @@ fit_velocity_pdf_bilateral <- function (velocity, axis="x", timestamp=0, dist="w
rmse=rmse(estimated_density, actual_density))
}
-fit_velocity_pdf_v2 <- function (x, sign=-1, plot=TRUE, axis="x") {
- #fit <- fitdist(x, "weibull")
- x_mode <- getmode(x)
- fit <- fitdist(x-x_mode, "weibull2", start=list(a1=1, b1=1, c1=1, a2=1, b2=1, c2=1))
- if (plot) {
- plot(fit)
- }
- coef(fit)
-}
-
-fit_velocity_acf <- function (velocity, plot=TRUE, axis="x", timestamp=0, write=TRUE) {
- acf_x <- acf(velocity, type="covariance", plot=FALSE)
- acf_x <- data.frame(acf=acf_x$acf, lag=acf_x$lag)
- model <- nls(acf ~ a*exp(-abs(b*lag)**c),
- data=acf_x,
- start=list(a=1,b=1,c=1),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(0,0,1e-6), upper=c(1000,1000,10))
- c <- summary(model)$coefficients
- acf_est <- c[[1]]*exp(-abs(c[[2]]*acf_x$lag)**c[[3]])
- if (plot) {
- plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('ACF', axis))
- lines(acf_x$lag, acf_est, col='red')
- }
- if (write) {
- path <- file.path('build', 'acf', axis)
- make_directory(path)
- filename <- file.path(path, sprintf("%d", timestamp))
- write.table(data.frame(lag=acf_x$lag,actual=acf_x$acf,estimated=acf_est),
- filename, row.names=FALSE, quote=FALSE)
- }
- list(coefficients=c,
- rmse=rmse(acf_est,acf_x$acf),
- lag=acf_x$lag,
- actual=acf_x$acf,
- estimated=acf_est)
-}
-
fit_velocity_acf_bilateral <- function (velocity, axis="x", timestamp=0, plot=TRUE) {
v = velocity
model_1 <- fit_velocity_acf(v[v<=0], axis=axis, timestamp=timestamp, plot=FALSE, write=FALSE);
@@ -350,35 +74,6 @@ fit_velocity_acf_bilateral <- function (velocity, axis="x", timestamp=0, plot=TR
estimated=estimated)
}
-fit_direction <- function (theta, density, plot=TRUE) {
- pdf_x <- data.frame(x=theta,density=density);
- #print(pdf_x)
- scale <- max(abs(theta))
- c <- tryCatch({
- model <- nls(density ~ dvonmises(2*pi*x/scale-m,0,a,b,c),
- data=pdf_x,
- start=list(a=1,b=1,c=1,m=0),
- control=list(warnOnly=TRUE),
- algorithm="port",
- lower=c(0,1,0,-10),
- upper=c(1000,1000,1000,10))
- c <- summary(model)$coefficients
- },error=function(cond) {
- message(cond)
- c(1,1,1)
- })
- pdf_x_est <- dvonmises(2*pi*pdf_x$x/scale,c[[4]],c[[1]],c[[2]],c[[3]])
- if (plot) {
- plot(pdf_x$x, pdf_x$density, xlab='Direction', ylab='PDF')
- lines(pdf_x$x, pdf_x_est, col='red')
- }
- list(coefficients=c,
- actual_density=pdf_x$density,
- estimated_density=pdf_x_est,
- x=pdf_x$x,
- residual=max(abs(pdf_x$density-pdf_x_est)))
-}
-
fit_direction_bilateral <- function (direction, plot=TRUE, timestamp=0, dist="vonmises") {
hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000))
x = hist$x;
diff --git a/R/common.R b/R/common.R
@@ -96,6 +96,8 @@ select_samples <- function (timestamp, time_delta, url="samples/load-cell.sqlite
velocity
}
+datetime_format <- "%FT%H:%M"
+
time_string_to_timestamp <- function(str) {
as.integer(as.POSIXct(strptime(str, datetime_format)))
}
@@ -103,3 +105,20 @@ time_string_to_timestamp <- function(str) {
make_directory <- function (path) {
dir.create(path, recursive=TRUE, showWarnings=FALSE)
}
+
+# https://stackoverflow.com/questions/32370485/convert-radians-to-degree-degree-to-radians
+rad2deg <- function(rad) {(rad * 180) / (pi)}
+deg2rad <- function(deg) {(deg * pi) / (180)}
+
+# https://stackoverflow.com/questions/27455948/in-r-whats-the-simplest-way-to-scale-a-vector-to-a-unit-vector
+unit <- function(x) {
+ length = sqrt(sum(x**2))
+ ifelse(length == 0, x, x/length)
+}
+
+# normalized rmse
+rmse <- function(estimated, actual) {
+ s_max <- max(actual)
+ s_min <- min(actual)
+ sqrt(mean((estimated - actual)**2)) / (s_max - s_min)
+}
diff --git a/R/distributions.R b/R/distributions.R
@@ -0,0 +1,273 @@
+dweibull2 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g=0) {
+ #ifelse(x<g,
+ # a1* abs(c1*b1)*(abs(b1*(x-g))**(c1-1))* exp(-(abs(b1*(x-g))**c1)),
+ # a1* abs(c2*b2)*(abs(b2*(x-g))**(c2-1))* exp(-(abs(b2*(x-g))**c2)))
+ ifelse(x<g,
+ a1*abs(b1*c1)*exp(-(abs(b1*(abs(x-g)-abs((c1-1)/c1)**(1/c1)))**c1)),
+ a2*abs(b2*c2)*exp(-(abs(b2*(abs(x-g)-abs((c2-1)/c2)**(1/c2)))**c2)))
+}
+
+dweibull4 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g1=0,
+ a3=1, b3=1, c3=1, a4=1, b4=1, c4=1, g2=0) {
+ ifelse(x<0,
+ dweibull2(x,a1,b1,c1,a2,b2,c2,g1),
+ dweibull2(x,a3,b3,c3,a4,b4,c4,g2))
+}
+
+dvonmises <- function (x, x_mean=0, a=1, b=1, c=1) {
+ a * exp(b*cos(c*(x-x_mean))) / (2*pi*besselI(b,0))
+}
+
+dgumbel <- function (x, location=0, scale=1) {
+ z <- (x-location)/scale
+ (1/scale) * exp(-(z + exp(-z)))
+}
+
+# https://www.tutorialspoint.com/r/r_mean_median_mode.htm
+# https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode
+getmode <- function(v) {
+ #uniqv <- unique(v)
+ #uniqv[which.max(tabulate(match(v, uniqv)))]
+ as.numeric(names(sort(-table(v)))[1])
+}
+
+# https://stackoverflow.com/questions/27418461/calculate-the-modes-in-a-multimodal-distribution-in-r
+find_peaks_indices <- function (x) {
+ modes <- c()
+ for (i in 2:(length(x)-1)){
+ if ((x[i] > x[i-1]) & (x[i] > x[i+1])) {
+ modes <- c(modes,i)
+ }
+ }
+ modes
+}
+
+find_real_peaks_indices <- function (x, max_modes=2) {
+ n <- 1024
+ #x_density <- density(x,n=n,adjust=0.125)
+ #y <- x_density$y
+ y <- x
+ modes <- find_peaks_indices(y)
+ # find the first N highest peaks
+ modes_xy <- data.frame(x=modes,y=y[modes])
+ modes_xy <- modes_xy[order(modes_xy$y, decreasing=TRUE),]
+ modes_xy <- modes_xy[c(1:max_modes), ]
+ # scale indices to [0,1]
+ (modes_xy$x-1)/(n-1)
+}
+
+find_max_x <- function (x, y) {
+ n <- length(y)-1
+ indices <- c(0:n)
+ index <- indices[y == max(y)]/n
+ x[floor(index*(length(x)-1)+1)]
+}
+
+normalize_hist <- function(pdf_x) {
+ #data.frame(density=pdf_x$density/max(abs(pdf_x$density)),x=pdf_x$mids)
+ data.frame(density=pdf_x$density,x=pdf_x$mids)
+}
+
+fit_velocity_pdf <- function (xx, x, density, sign=-1, plot=TRUE, axis="x", dist="weibull2") {
+ a0 <- 0
+ b0 <- 0
+ c0 <- 1
+ c1 <- 3
+ indices <- find_real_peaks_indices(x, max_modes=1)
+ if (FALSE) {
+ pdf <- normalize_hist(hist(x, plot=FALSE, breaks=1000))
+ print(indices)
+ indices = floor(indices * length(pdf$x))
+ x_modes <- pdf$x[indices]
+ print(x_modes)
+ model <- nls(pdf$density ~ dweibull4(pdf$x,a1,b1,c1,a2,b2,c2,x_modes[[1]],a3,b3,c3,a4,b4,c4,x_modes[[2]]),
+ data=pdf,
+ start=list(a1=1,b1=1,c1=1,a2=1,b2=1,c2=1,
+ a3=1,b3=1,c3=1,a4=1,b4=1,c4=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(a0,b0,c0,a0,b0,c0,a0,b0,c0,a0,b0,c0),
+ upper=c(1000,1000,3,1000,2000,3,1000,1000,3,1000,2000,3))
+ c = summary(model)$coefficients
+ if (plot) {
+ pdf_x_est <- dweibull4(pdf$x,c[[1]],c[[2]],c[[3]],c[[4]],c[[5]],c[[6]],x_modes[[1]],
+ c[[7]],c[[8]],c[[9]],c[[10]],c[[11]],c[[12]],x_modes[[2]])
+ plot(pdf$x, pdf$density, xlab=paste('Velocity', axis), ylab='PDF')
+ lines(pdf$x, pdf_x_est, col='red')
+ }
+ c
+ } else {
+ #x_mode <- getmode(x)
+ #pdf_x <- normalize_hist(hist(x, plot=FALSE, breaks=100))
+ pdf_x <- data.frame(x=xx,density=x)
+ c <- c()
+ pdf_x_set <- c()
+ if (dist == "weibull2") {
+ x_mode <- find_max_x(pdf_x$x,pdf_x$density)
+ #indices = floor(indices * length(pdf_x$x)) + 1
+ #x_mode <- pdf_x$x[indices][[1]]
+ c <- tryCatch(
+ {
+ model <- nls(pdf_x$density ~ dweibull2(pdf_x$x,a,b,c,d,e,f,g=x_mode),
+ data=pdf_x,
+ start=list(a=1,b=1,c=1,d=1,e=1,f=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(a0,b0,c0,a0,b0,c0), upper=c(1000,1000,c1,1000,2000,c1))
+ summary(model)$coefficients
+ },
+ error=function(cond) {
+ message(cond)
+ c(1,1,1,1,1,1)
+ })
+ pdf_x_est <- dweibull2(pdf_x$x,c[[1]],c[[2]],c[[3]],c[[4]],c[[5]],c[[6]],x_mode)
+ } else if (dist == "vonmises") {
+ c <- tryCatch({
+ model <- nls(density ~ dvonmises(2*pi*pdf_x$x/(max(pdf_x$x)),0,a,b,c),
+ data=pdf_x,
+ start=list(a=1,b=1,c=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(0,1,0), upper=c(1000,1000,pi))
+ c <- summary(model)$coefficients
+ },error=function(cond) {
+ message(cond)
+ c(1,1,1)
+ })
+ pdf_x_est <- dvonmises(2*pi*pdf_x$x/max(pdf_x$x),0,c[[1]],c[[2]],c[[3]])
+ } else if (dist == "gumbel") {
+ c <- tryCatch(
+ {
+ model <- nls(pdf_x$density ~ c*dgumbel(pdf_x$x,location=a,scale=b),
+ data=pdf_x,
+ start=list(a=1,b=1,c=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(-1000,0.1,c=1), upper=c(1000,1000,1000))
+ c = summary(model)$coefficients
+ },
+ error=function(cond) {
+ message(cond)
+ c <- c(1,1,1,1)
+ })
+ pdf_x_est <- c[[3]]*dgumbel(pdf_x$x,location=c[[1]],scale=c[[2]])
+ } else if (dist == "norm") {
+ fit <- fitdist(pdf_x$x, "norm")
+ c <- coef(fit)
+ pdf_x_est <- dnorm(pdf_x$x,mean=c[[1]],sd=c[[2]])
+ } else if (dist == "rayleigh") {
+ c <- tryCatch(
+ {
+ model <- nls(pdf_x$density ~ b*dweibull(sign*pdf_x$x,scale=a,shape=2),
+ data=pdf_x,
+ start=list(a=1,b=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(0.1,0), upper=c(1000,1000))
+ c = summary(model)$coefficients
+ },
+ error=function(cond) {
+ message(cond)
+ c <- c(1,1,1,1)
+ })
+ pdf_x_est <- c[[2]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=2)
+ } else {
+ x_mode <- find_max_x(pdf_x$x,pdf_x$density)
+ c <- tryCatch(
+ {
+ model <- nls(pdf_x$density ~ c*dweibull(sign*pdf_x$x,scale=a,shape=b),
+ data=pdf_x,
+ start=list(a=1,b=1,c=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(0.1,0.1,0),
+ upper=c(1000,1000,1000))
+ c = summary(model)$coefficients
+ },
+ error=function(cond) {
+ message(cond)
+ c <- c(1,1,1,1)
+ })
+ pdf_x_est <- c[[3]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=c[[2]])
+ }
+ #if (plot) {
+ # plot(sign*pdf_x$x, pdf_x$density, xlab=paste('Velocity', axis), ylab='PDF')
+ # lines(density(sign*x), col='blue')
+ # lines(sign*pdf_x$x, pdf_x_est, col='red')
+ #}
+ #if (sign < 0) {
+ # print(pdf_x$x)
+ # pdf_x$x = -rev(pdf_x$x)
+ # pdf_x$density = rev(pdf_x$density)
+ # pdf_x_est = rev(pdf_x_est)
+ #}
+ list(coefficients=c,
+ actual_density=pdf_x$density,
+ estimated_density=pdf_x_est,
+ x=pdf_x$x,
+ residual=max(abs(pdf_x$density-pdf_x_est)))
+ }
+}
+
+fit_velocity_acf <- function (velocity, plot=TRUE, axis="x", timestamp=0, write=TRUE,
+ debug=FALSE) {
+ acf_x <- acf(velocity, type="covariance", plot=FALSE)
+ acf_x <- data.frame(acf=acf_x$acf, lag=acf_x$lag)
+ if (debug) {
+ plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('DEBUG ACF', axis))
+ }
+ model <- nls(acf ~ a*exp(-abs(b*lag)**c),
+ data=acf_x,
+ start=list(a=1,b=1,c=1),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(0,0,1e-6), upper=c(1000,1000,10))
+ c <- summary(model)$coefficients
+ acf_est <- c[[1]]*exp(-abs(c[[2]]*acf_x$lag)**c[[3]])
+ if (plot) {
+ plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('ACF', axis))
+ lines(acf_x$lag, acf_est, col='red')
+ }
+ if (write) {
+ path <- file.path('build', 'acf', axis)
+ make_directory(path)
+ filename <- file.path(path, sprintf("%d", timestamp))
+ write.table(data.frame(lag=acf_x$lag,actual=acf_x$acf,estimated=acf_est),
+ filename, row.names=FALSE, quote=FALSE)
+ }
+ list(coefficients=c,
+ rmse=rmse(acf_est,acf_x$acf),
+ lag=acf_x$lag,
+ actual=acf_x$acf,
+ estimated=acf_est)
+}
+
+fit_direction <- function (theta, density, plot=TRUE) {
+ pdf_x <- data.frame(x=theta,density=density);
+ #print(pdf_x)
+ scale <- max(abs(theta))
+ c <- tryCatch({
+ model <- nls(density ~ dvonmises(2*pi*x/scale-m,0,a,b,c),
+ data=pdf_x,
+ start=list(a=1,b=1,c=1,m=0),
+ control=list(warnOnly=TRUE),
+ algorithm="port",
+ lower=c(0,1,0,-10),
+ upper=c(1000,1000,1000,10))
+ c <- summary(model)$coefficients
+ },error=function(cond) {
+ message(cond)
+ c(1,1,1)
+ })
+ pdf_x_est <- dvonmises(2*pi*pdf_x$x/scale,c[[4]],c[[1]],c[[2]],c[[3]])
+ if (plot) {
+ plot(pdf_x$x, pdf_x$density, xlab='Direction', ylab='PDF')
+ lines(pdf_x$x, pdf_x_est, col='red')
+ }
+ list(coefficients=c,
+ actual_density=pdf_x$density,
+ estimated_density=pdf_x_est,
+ x=pdf_x$x,
+ residual=max(abs(pdf_x$density-pdf_x_est)))
+}
+
diff --git a/R/vtestbed.R b/R/vtestbed.R
@@ -0,0 +1,73 @@
+source("R/common.R")
+source("R/distributions.R")
+
+# the timestamp against which the simulation output is validated
+time_start = 1616936400
+time_end = time_start
+time_delta = 60*60*2
+
+data <- read.table("samples/vtestbed/statistics.log", sep=",", header=TRUE)
+#data <- data[c("time","wind_velocity_x","wind_velocity_y","wind_velocity_z")]
+#plot(data$time, data$wind_velocity_x)
+#plot(data$time, data$wind_velocity_y)
+#plot(data$time, data$wind_velocity_z)
+data$speed <- speed(data$wind_velocity_x,data$wind_velocity_y,data$wind_velocity_z)
+data$direction <- direction(data$wind_velocity_x,data$wind_velocity_y)
+#hist(data$speed)
+#hist(data$direction)
+
+
+do_fit_velocity <- function (velocity) {
+ velocity_hist <- normalize_hist(hist(velocity, plot=FALSE, breaks=100))
+ x = velocity_hist$x;
+ v = velocity_hist$density;
+ d = density(velocity,adjust=0.5)
+ dist <- fit_velocity_pdf(x, v, density=d$y, sign=1, axis="v", dist="weibull")
+ plot(x, dist$actual_density, xlab='Wind speed', ylab='PDF')
+ lines(d, col='green', lwd=2)
+ lines(dist$x, dist$estimated_density, col='red', lwd=2)
+ dist
+}
+
+do_fit_direction <- function (direction) {
+ hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000))
+ x = hist$x;
+ v = hist$density;
+ d = density(direction,adjust=0.5)
+ dist <- fit_direction(x, v, plot=FALSE)
+ plot(x, dist$actual_density, xlab='Direction', ylab='PDF')
+ lines(d, col='green', lwd=2)
+ lines(dist$x, dist$estimated_density, col='red', lwd=2)
+ dist
+}
+
+do_fit_acf <- function (velocity, axis="x") {
+ v = velocity
+ model <- fit_velocity_acf(v, axis=axis, plot=FALSE, write=FALSE);
+ plot(model$lag, model$actual, xlab='Lag', ylab=paste('ACF ', axis, sep=""))
+ lines(model$lag, model$estimated, col='red', lwd=2)
+ model
+}
+
+# https://stackoverflow.com/questions/19918985/r-plot-only-text
+banner <- function (str) {
+ plot(c(0, 1), c(0, 1), ann = F, bty = 'n', type = 'n', xaxt = 'n', yaxt = 'n')
+ text(x = 0.5, y = 0.5, str, cex = 1.6, col = "black")
+}
+
+banner("Virtual testbed")
+dist_speed <- do_fit_velocity(data$speed)
+dist_direction <- do_fit_direction(data$direction)
+#model_acf_x <- do_fit_acf(data$wind_velocity_x, axis="X")
+model_acf_y <- do_fit_acf(data$wind_velocity_y, axis="Y")
+#model_acf_z <- do_fit_acf(data$wind_velocity_z, axis="Z")
+
+banner("Anemometer")
+data_real <- select_samples(time_start, time_delta)
+dist_speed_real <- do_fit_velocity(data_real$speed)
+dist_direction_real <- do_fit_direction(data_real$direction)
+model_acf_x <- do_fit_acf(data_real$x[data_real$x<=0], axis="X")
+model_acf_y <- do_fit_acf(data_real$y[data_real$y<=0], axis="Y")
+model_acf_z <- do_fit_acf(data_real$z[data_real$z<=0], axis="Z")
+
+#print(data)
