commit 8b541d32182547c2d64cf781720c78d1d72b9e5b
parent 6c53bbcaaf1cd4f409d0ab1c96bd67e352c47140
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Tue, 6 Apr 2021 13:07:04 +0300
Bivariate distribution.
Diffstat:
7 files changed, 156 insertions(+), 91 deletions(-)
diff --git a/Makefile b/Makefile
@@ -56,13 +56,13 @@ build/daily-stats.svg: gnuplot/daily-stats.gnuplot build/daily-stats
build/daily-rmse.svg: gnuplot/daily-rmse.gnuplot build/daily-rmse
gnuplot -d $<
-build/gnuplot/direction-dist-min.svg: gnuplot/direction-dist.gnuplot
+build/gnuplot/direction-dist-min.svg: gnuplot/direction-dist.gnuplot build/direction/1616497200-hist-xy
gnuplot -d -c gnuplot/direction-dist.gnuplot \
- build/direction/1616497200-hist-xy build/gnuplot/direction-dist-min.svg
+ build/direction/1616497200-hist-xy build/gnuplot/direction-dist-min.svg 0.024544
-build/gnuplot/direction-dist-max.svg: gnuplot/direction-dist.gnuplot
+build/gnuplot/direction-dist-max.svg: gnuplot/direction-dist.gnuplot build/direction/1616792400-hist-xy
gnuplot -d -c gnuplot/direction-dist.gnuplot \
- build/direction/1616792400-hist-xy build/gnuplot/direction-dist-max.svg
+ build/direction/1616792400-hist-xy build/gnuplot/direction-dist-max.svg 0.154791
#build/main.zip: build/gnuplot/*.eps main.tex
build/main.zip: main.tex
diff --git a/R/anal.R b/R/anal.R
@@ -1,4 +1,3 @@
-library(DBI)
library(fitdistrplus)
library(plotrix)
library(sn)
@@ -21,12 +20,6 @@ unit <- function(x) {
ifelse(length == 0, x, x/length)
}
-time_string_to_timestamp <- function(str) {
- as.integer(as.POSIXct(strptime(str, datetime_format)))
-}
-
-speed <- function(x,y,z=0) { sqrt(x**2 + y**2 + z**2) }
-direction <- function(x,y) { atan2(y,x) }
# retain the original sign
power <- function(x,y) {
ifelse(abs(x) < 1, sign(x)*abs(x)**y, x**y)
@@ -388,7 +381,6 @@ timestamps <- seq(time_start, time_end, time_delta)
#i <- 0
#for (timestamp in timestamps) {
cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='remove') %dopar% {
- db <- dbConnect(RSQLite::SQLite(), "samples/load-cell.sqlite3")
cumulative_data <- data.frame(speed=numeric(0),direction=numeric(0),
x_mode=numeric(0),abs_x_mode=numeric(0),x_sd=numeric(0),
y_mean=numeric(0),
@@ -402,29 +394,13 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
acf_y_a=numeric(0),acf_y_b=numeric(0),acf_y_c=numeric(0),
timestamp=numeric(0),
x_rmse=numeric(0), y_rmse=numeric(0), z_rmse=numeric(0))
- velocity = dbGetQuery(db,
- "SELECT x,y,z FROM samples WHERE timestamp BETWEEN :t0 AND (:t0 + :dt) ORDER BY timestamp",
- 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",
- params = list(t0 = timestamp, dt=time_delta,
- xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z)))
+ velocity <- select_samples(timestamp, time_delta)
if (nrow(velocity) == 0) { return(NULL) }
filename <- sprintf("build/rplot-%d.pdf", timestamp)
print(filename)
pdf(filename)
- old_mean_x <- max(abs(velocity$x)) - mean(velocity$x)
- coefficients <- calibration_coefficients()
- velocity$x = sample_to_speed(velocity$x, coefficients$x1, coefficients$x2)
- velocity$y = sample_to_speed(velocity$y, coefficients$y1, coefficients$y2)
- velocity$z = sample_to_speed(velocity$z, coefficients$z1, coefficients$z2)
x_mode <- getmode(velocity$x)
y_mode <- getmode(velocity$y)
- velocity$speed = speed(velocity$x,velocity$y,velocity$z)
- velocity$direction = direction(velocity$x,velocity$y)
new_row = nrow(cumulative_data)+1
cumulative_data[new_row, "timestamp"] = timestamp
cumulative_data[new_row, "x_mode"] = x_mode
@@ -483,25 +459,23 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
path <- file.path('build', 'direction')
dir.create(path, recursive=TRUE)
filename <- file.path(path, sprintf("%d-hist-xy", timestamp))
- #
- hist_xy <- hist2d(velocity$x, velocity$y, same.scale=TRUE, show=TRUE)
- grid_xy <- expand.grid(x=hist_xy$x,y=hist_xy$y)
- #
# add corner points with mean speed
- #x_max = max(abs(velocity$x))
- #y_max = max(abs(velocity$y))
- x_max = max(abs(grid_xy$x))
- y_max = max(abs(grid_xy$y))
+ x_max = max(abs(velocity$x))
+ y_max = max(abs(velocity$y))
s = mean(velocity$speed)
r = sqrt(x_max**2 + y_max**2)
- #write.table(data.frame(x=c(velocity$x/r,r,-r,-r,r),
- # y=c(velocity$y/r,r,r,-r,-r),
- # r=c(velocity$speed,s,s,s,s)),
- # filename, row.names=FALSE, quote=FALSE)
- write.table(data.frame(x=c(grid_xy$x/r,r,-r,-r,r),
- y=c(grid_xy$y/r,r,r,-r,-r),
- r=c(as.vector(hist_xy$counts),0,0,0,0)),
+ write.table(data.frame(x=c(velocity$x/r,r,-r,-r,r),
+ y=c(velocity$y/r,r,r,-r,-r),
+ r=c(velocity$speed,s,s,s,s)),
filename, row.names=FALSE, quote=FALSE)
+ #hist_xy <- hist2d(velocity$x, velocity$y, same.scale=TRUE, show=TRUE)
+ #grid_xy <- expand.grid(x=hist_xy$x,y=hist_xy$y)
+ #x_max = max(abs(grid_xy$x))
+ #y_max = max(abs(grid_xy$y))
+ #write.table(data.frame(x=c(grid_xy$x/r,r,-r,-r,r),
+ # y=c(grid_xy$y/r,r,r,-r,-r),
+ # r=c(as.vector(hist_xy$counts),0,0,0,0)),
+ # filename, row.names=FALSE, quote=FALSE)
polar.plot(direction_hist$direction,rad2deg(direction_hist$angle),
start=45,clockwise=TRUE,rp.type="polygon")
par(mar=old_mar)
@@ -511,7 +485,6 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
cumulative_data[new_row, "y_mean"] = mean(velocity$y)
#i <- i + 1
#setTxtProgressBar(pb, i)
- dbDisconnect(db)
dev.off()
cumulative_data
}
diff --git a/R/common.R b/R/common.R
@@ -1,3 +1,5 @@
+library(DBI)
+
# https://stackoverflow.com/questions/13023274/how-to-do-printf-in-r
printf <- function(...) cat(sprintf(...))
@@ -54,3 +56,34 @@ sample_to_speed <- function(x, c1, c2) {
x
}
+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)))
+ if (nrow(velocity) == 0) { return(data.frame()) }
+ # convert sensor values to wind speed in m/s
+ coefficients <- calibration_coefficients()
+ velocity$x = sample_to_speed(velocity$x, coefficients$x1, coefficients$x2)
+ velocity$y = sample_to_speed(velocity$y, coefficients$y1, coefficients$y2)
+ velocity$z = sample_to_speed(velocity$z, coefficients$z1, coefficients$z2)
+ # compute wind speed and direction
+ velocity$speed = speed(velocity$x,velocity$y,velocity$z)
+ velocity$direction = direction(velocity$x,velocity$y)
+ dbDisconnect(db)
+ velocity
+}
+
+time_string_to_timestamp <- function(str) {
+ as.integer(as.POSIXct(strptime(str, datetime_format)))
+}
diff --git a/R/direction-min-max.R b/R/direction-min-max.R
@@ -0,0 +1,41 @@
+library(gplots) # hist2d
+
+source("R/common.R")
+
+t_min_rmse <- 1616497200
+t_max_rmse <- 1616792400
+time_delta = 60*60*2
+for (timestamp in c(t_min_rmse, t_max_rmse)) {
+ velocity <- select_samples(timestamp, timestamp+time_delta)
+ path <- file.path('build', 'direction')
+ dir.create(path, recursive=TRUE)
+ filename <- file.path(path, sprintf("%d-hist-xy", timestamp))
+ print(filename)
+ # add corner points with mean speed
+ ## v1
+ #x_max = max(abs(velocity$x))
+ #y_max = max(abs(velocity$y))
+ #s = mean(velocity$speed)
+ #r = sqrt(x_max**2 + y_max**2)
+ #tmp <- data.frame(x=c(velocity$x/r,r,-r,-r,r),
+ # y=c(velocity$y/r,r,r,-r,-r),
+ # r=c(velocity$speed,s,s,s,s))
+ #tmp <- tmp[order(tmp$r),] # order by speed for gnuplot
+ #write.table(tmp, filename, row.names=FALSE, quote=FALSE)
+
+# write.table(data.frame(x=velocity$x,
+# y=velocity$y,
+# r=velocity$speed),
+# filename, row.names=FALSE, quote=FALSE)
+
+ # v2
+ hist_xy <- hist2d(velocity$x, velocity$y, same.scale=TRUE, show=TRUE, nbins=400)
+ grid_xy <- expand.grid(x=-hist_xy$x,y=-hist_xy$y)
+ x_max = max(abs(grid_xy$x))
+ y_max = max(abs(grid_xy$y))
+ r = sqrt(x_max**2 + y_max**2)
+ write.table(data.frame(x=c(grid_xy$x/r,r,-r,-r,r),
+ y=c(grid_xy$y/r,r,r,-r,-r),
+ r=c(as.vector(t(hist_xy$counts)),0,0,0,0)),
+ filename, row.names=FALSE, quote=FALSE)
+}
diff --git a/gnuplot/direction-dist.gnuplot b/gnuplot/direction-dist.gnuplot
@@ -1,6 +1,8 @@
load 'gnuplot/blues.pal'
set palette defined (0 '#4040c0', 1 '#c0c0c0', 2 '#c04040')
+#set palette defined (0 '#c0c0c0', 1 '#c04040')
set terminal svg size 450/2,450/2 font 'Free Serif, 10' enhanced round dynamic
+#set terminal png size 450/2,450/2 font 'Free Serif, 10' enhanced
set xtics nomirror out offset 0,0.25
set ytics nomirror out offset 0.5,0
#set border 2+8
@@ -13,13 +15,18 @@ unset key
#set xlabel 'υ_x'
#set ylabel 'υ_y'
set cbrange [0:*]
-set lmargin at screen 0.05
-set rmargin at screen 0.85
-set bmargin at screen 0.1
-set tmargin at screen 0.9
+set lmargin 2
+set tmargin 2.5
+set bmargin 1.5
+set rmargin 6
+set size square
+#set lmargin at screen 0.05
+#set rmargin at screen 0.85
+#set bmargin at screen 0.1
+#set tmargin at screen 0.9
-positive(x) = x<0 ? (1/0) : x
-negative(x) = x>0 ? (1/0) : x
+positive(x) = x<=25 ? NaN : x
+negative(x) = x>0 ? NaN : x
round(x) = floor(x+0.5)
length(x,y) = sqrt(x**2 + y**2)
angle(x,y) = atan2(y,x)
@@ -30,58 +37,69 @@ max(x,y) = x>y ? x : y
filename = ARG1
filename_tmp = sprintf('%s.tmp', filename)
filename_cnt = sprintf('%s.cnt', filename)
-filename_svg = ARG2
-stats filename using 1:3
-x_max = max(abs(STATS_min_x),abs(STATS_max_x))
-y_max = max(abs(STATS_min_y),abs(STATS_max_y))
-print(y_max)
+filename_out = ARG2
+rmse = ARG3
# Min. direction RMSE = 0.024544 (1616497200), mean direction = -2.399689
# Max. direction RMSE = 0.154791 (1616792400), mean direction = -2.440686
-# Mean direction RMSE = 0.056350ED6',\
+# Mean direction RMSE = 0.056350
-# http://www.gnuplotting.org/circular-heat-map/
-set table filename_tmp
-set dgrid3d 100,100
+## interpolation {{{
+## http://www.gnuplotting.org/circular-heat-map/
+#set table filename_tmp
+#set dgrid3d 100,100
set xrange [-1:1]
set yrange [-1:1]
-set colorbox user origin 0.925,0.1 size 0.03,0.8
-splot filename using 1:2:3
-unset table
-
-set table filename_cnt
-set contour base
-set cntrparam levels 5
-set cntrparam order 8 bspline
-#set cntrlabel onecolor format '%8.3g' font ',7' start 25 interval -1
-unset surface
-splot filename using 1:2:3 with lines
-unset table
-
-set output filename_svg
+#splot filename using 1:2:3
+#unset table
+## }}}
+#
+## contours {{{
+#set table filename_cnt
+#set contour base
+#set cntrparam levels 5
+#set cntrparam order 8 bspline
+##set cntrlabel onecolor format '%8.3g' font ',7' start 25 interval -1
+#unset surface
+#set dgrid3d 100,100
+#splot filename using 1:2:3 with lines
+#unset table
+## }}}
+set output filename_out
set multiplot
-plot filename_tmp using 1:2:(circle($1,$2,$3)) with image
-plot filename_cnt using (circle($1,$2,$1)):2 with lines lc '#404040'
-#plot 'build/1616792400-hist-xy.tmp' using 1:2 with points
+# polar grid {{{
# now plot the polar grid only
# https://stackoverflow.com/questions/18792461/gnuplot-2d-polar-plot-with-heatmap-from-3d-dataset-possible
r = 3
unset contour
-set border polar lw 5 lc rgb '#404040'
+set border polar lw 1 lc rgb '#404040' back
set surface
set angles degree
-set style line 11 lc rgb '#404040' lw 1 dashtype 2
-#set grid polar ls 11 nortics
+set style line 11 lc rgb '#404040' lw 0.5 dashtype 2
+set grid polar ls 11 nortics back
set polar
set rrange[0:r]
unset raxis
set rtics format '' scale 0
unset parametric
-set for [i=0:330:30] label at first (r+0.40)*cos(i), first (r+0.40)*sin(i)\
-center sprintf('%d', i)
+set for [i=0:330:30] label at first (r+0.40)*cos(i), first (r+0.40)*sin(i) center sprintf('%d', i)
plot NaN w l
-unset multiplot
+unset polar
+unset grid
+unset border
+# }}}
+#plot filename_tmp using 1:2:(circle($1,$2,$3)) with image
+set border lw 1 1+2+4+8
+unset border
+set colorbox user origin 0.900,0.1 size 0.03,0.8
+set cbtics offset -0.5,0
+s = 0.5
+set xrange [-1*s:1*s]
+set yrange [-1*s:1*s]
+set title offset 0,0.5 sprintf('RMSE = %.0f%%', round(rmse*100))
+#plot filename using 1:2:(positive($3)) with points pt 7 ps 0.25 linecolor palette z
+plot filename using 1:2:(positive($3)) with image
+#plot filename_cnt using (circle($1,$2,$1)):2 with lines lc '#404040'
-#splot 'build/direction/1616497200-hist-xy' using 1:2:3 with pm3d
diff --git a/main.tex b/main.tex
@@ -352,12 +352,12 @@ and high wind speeds. TODO
\begin{figure}
\centering
- \includegraphics{build/gnuplot/direction-dist-max.png}~
+ \includegraphics{build/gnuplot/direction-dist-max.png}~\phantom{x}~
\includegraphics{build/gnuplot/direction-dist-min.png}
- \caption{Per-axis wind velocity distributions fitted into Weibull
- distribution.
- Data intervals with the smallest error (the first row),
- data intervals with the largest error (the second row). Red line shows estimated
+ \caption{Bivariate \(x\) and \(y\) velocity distribution
+ in polar coordinates. Data interval with the largest error (top left), data
+ interval with the smallest error (top right).
+ Red line shows estimated
probability density of positive wind speed projections, blue line shows estimated
probability density of negative wind speed projections and circles denote observed
probability density of wind speed projections.
diff --git a/samples/load-cell.sqlite3 b/samples/load-cell.sqlite3
@@ -1 +1 @@
-../.git/annex/objects/JP/V6/SHA256E-s128065536--42760ec619c2410c0b6302f3149aad6229950b8ad1f604823b7ad36969597987/SHA256E-s128065536--42760ec619c2410c0b6302f3149aad6229950b8ad1f604823b7ad36969597987-
\ No newline at end of file
+../.git/annex/objects/PP/jm/SHA256E-s169267200--fa33d5863bdcf6313e381e3f04835b45e2d97a55365e70375d1a0650423a309a/SHA256E-s169267200--fa33d5863bdcf6313e381e3f04835b45e2d97a55365e70375d1a0650423a309a+
\ No newline at end of file
