arma-thesis

Simulation modelling of irregular waves for marine object dynamics programmes
git clone https://git.igankevich.com/arma-thesis.git
Log | Files | Refs | LICENSE

common.R (9672B)


      1 source(file.path("R", "waves.R"))
      2 source(file.path("R", "transform.R"))
      3 
      4 arma.qqplot_grid <- function (dir, params, titles, ...) {
      5   wave_params <- arma.load_wave_parameters(dir, params)
      6   i <- 1
      7   for (name in names(wave_params)) {
      8     arma.qqplot(wave_params[[name]], 100, titles[[i]], ...)
      9     i <- i + 1
     10   }
     11 }
     12 
     13 arma.qqplot_grid_adj <- function (dir, params, titles, adj, ...) {
     14   wave_params <- arma.load_wave_parameters(dir, params)
     15   i <- 1
     16   for (name in names(wave_params)) {
     17     ttl = list(title=titles[[i]], adjust=adj)
     18     arma.qqplot(wave_params[[name]], 100, ttl, ...)
     19     i <- i + 1
     20   }
     21 }
     22 
     23 arma.wavy_plot_matrix <- function (x, y, z, t, ...) {
     24 	nrz <- nrow(z)
     25 	ncz <- ncol(z)
     26 	# Create a function interpolating colors in the range of specified colors
     27 	jet.colors <- colorRampPalette( c("blue", "green") )
     28 	# Generate the desired number of colors from this palette
     29 	nbcol <- 100
     30 	color <- jet.colors(nbcol)
     31 	# Compute the z-value at the facet centres
     32 	zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
     33 	# Recode facet z-values into color indices
     34 	facetcol <- cut(zfacet, nbcol)
     35 	persp(x, y, z, phi=30, theta=30, col=color[facetcol], ...)
     36 }
     37 
     38 arma.to_matrix <- function (data, t) {
     39 	slice <- data[data$t == t,]
     40 	x <- unique(slice$x)
     41 	y <- unique(slice$y)
     42 	z <- with(slice, {
     43 		n <- sqrt(length(z))
     44 		out <- matrix(nrow=n, ncol=n)
     45 		out[cbind(x, y)] <- z
     46 		out
     47 	})
     48 	list(x=x,y=y,z=z)
     49 }
     50 
     51 arma.wavy_plot <- function (data, t, ...) {
     52 	with(arma.to_matrix(data, t), arma.wavy_plot_matrix(x, y, z, ...))
     53 }
     54 
     55 arma.acf_plot <- function (data, t, ...) {
     56 	ifft <- function (arr) { fft(arr, inverse=TRUE) }
     57 	with(
     58 		arma.to_matrix(data, t),
     59 		function () {
     60 			# compute ACF using Wiener---Khinchin theorem
     61 			n <- nrow(z)*ncol(z)
     62 			z <- Re(ifft(abs(fft(z))^2))/n/n
     63 			# slice 1/4th of the matrix
     64 			n1 <- floor(nrow(z)/2)
     65 			n2 <- floor(ncol(z)/2)
     66 			x <- x[1:n1]
     67 			y <- y[1:n2]
     68 			z <- z[1:n1,1:n2]
     69 			arma.wavy_plot_matrix(x, y, z, ...)
     70 		}
     71 	)()
     72 }
     73 
     74 arma.skew_normal_1_plot <- function(x, params) {
     75   data <- mapply(
     76     function (s, k) arma.skew_normal_1(x, s, k),
     77     params$skewness,
     78     params$kurtosis
     79   )
     80   plot.new()
     81   plot.window(xlim=range(x),ylim=range(data))
     82   axis(1); axis(2); box()
     83   for (i in seq_len(ncol(data))) {
     84     d <- data[,i]
     85     lines(x, d, lty=paste(params$linetypes[[i]]))
     86   }
     87   title(xlab="x", ylab="y")
     88 }
     89 
     90 
     91 arma.skew_normal_2_plot <- function(x, params) {
     92   data <- mapply(
     93     function (a) arma.skew_normal_2(x, a),
     94     params$alpha
     95   )
     96   plot.new()
     97   plot.window(xlim=range(x),ylim=range(data))
     98   axis(1); axis(2); box()
     99   for (i in seq_len(ncol(data))) {
    100     d <- data[,i]
    101     lines(x, d, lty=paste(params$linetypes[[i]]))
    102   }
    103   title(xlab="x", ylab="y")
    104 }
    105 
    106 arma.fmt <- function(x, ndigits) {
    107   format(round(x, ndigits), nsmall=ndigits)
    108 }
    109 
    110 arma.plot_partitions <- function() {
    111 	library("rgl")
    112 	part_alpha <- 0.2
    113 	part_col <- "grey"
    114 	part_size <- 2
    115 	sc <- 0.8
    116 	off <- part_size * (1-sc)/2
    117 	bcol <- "grey"
    118 	balpha <- 1.0
    119 	sc2 <- 1.0 - sc
    120 	off2 <- part_size * sc/2
    121 	ccol <- "red"
    122 	calpha <- 1.0
    123 	# whole parts
    124 	a1 <- cube3d(color=part_col, alpha=part_alpha)
    125 	a2 <- cube3d(color=part_col, alpha=part_alpha)
    126 	a3 <- cube3d(color=part_col, alpha=part_alpha)
    127 	shade3d(translate3d(a1, 0*part_size, 0, 0))
    128 	shade3d(translate3d(a2, 1*part_size, 0, 0))
    129 	shade3d(translate3d(a3, 2*part_size, 0, 0))
    130 	# stripped parts
    131 	b1 <- scale3d(cube3d(color=bcol, alpha=balpha), sc, sc, sc)
    132 	b2 <- scale3d(cube3d(color=bcol, alpha=balpha), sc, sc, sc)
    133 	b3 <- scale3d(cube3d(color=bcol, alpha=balpha), sc, sc, sc)
    134 	shade3d(translate3d(b1, 0 + off, off, off))
    135 	shade3d(translate3d(b2, 2 + off, off, off))
    136 	shade3d(translate3d(b3, 4 + off, off, off))
    137 	# overlap intervals
    138 	c1 <- scale3d(cube3d(color=ccol, alpha=calpha), sc2, 1, 1)
    139 	c2 <- scale3d(cube3d(color=ccol, alpha=calpha), sc2, 1, 1)
    140 	c3 <- scale3d(cube3d(color=ccol, alpha=calpha), sc2, 1, 1)
    141 	shade3d(translate3d(c1, 0 + off2, 0, 0))
    142 	shade3d(translate3d(c2, 2 + off2, 0, 0))
    143 	shade3d(translate3d(c3, 4 + off2, 0, 0))
    144 }
    145 
    146 arma.plot_ramp_up_interval <- function(label="Ramp-up interval") {
    147 	zeta <- read.csv(file.path("build", "arma-benchmarks", "verification-orig", "standing_wave", "zeta.csv"))
    148 	t <- round(mean(zeta$t))
    149 	res <- arma.wavy_plot(zeta, t, scale=FALSE)
    150 	library("grDevices")
    151 	ax <- 7
    152 	ay <- 7
    153 	my <- max(zeta$y)
    154 	lines(trans3d(
    155 		c(0,  0, ax, ax, 0),
    156 		c(0, my, my,  0, 0),
    157 		c(0,  0,  0,  0, 0),
    158 		pmat=res
    159   ), col="red", lwd=3)
    160   text(trans3d(0, my/2, max(zeta$z)*1.5, pmat=res), label, col="red", font=2)
    161 	from <- trans3d(0, my/2, max(zeta$z)*1.4, pmat=res)
    162 	to <- trans3d(0, my/2, max(zeta$z)*0.05, pmat=res)
    163 	arrows(from$x, from$y, to$x, to$y, lwd=2, angle=10, length=0.1, col="red")
    164 }
    165 
    166 arma.cube <- function (x, y, z) {
    167   c <- cube3d(alpha=0.2)
    168   c$material$lwd <- 4
    169 #  c$material$front <- 'line'
    170   c$material$back <- 'line'
    171   c$material$col <- '#ffffff'
    172   shade3d(translate3d(c, x, y, z))
    173 }
    174 
    175 arma.arrow <- function (x1,x2,ps) {
    176   arrow3d(
    177     c(x1[1], x1[2], x1[3])*ps,
    178     c(x2[1], x2[2], x2[3])*ps,
    179     type="rotation", col="grey", s=1/7)
    180 }
    181 
    182 arma.plot_ar_cubes <- function(nx, ny, nz) {
    183   library("rgl")
    184   part_size <- 2
    185   # generate cubes
    186 #  for (i in c(1:nx)) {
    187 #    for (j in c(1:ny)) {
    188 #      for (k in c(1:nz)) {
    189 #        arma.cube(i*part_size, j*part_size, k*part_size)
    190 #      }
    191 #    }
    192 #  }
    193   # generate arrows
    194   for (i in c(1:nx)) {
    195     for (j in c(1:ny)) {
    196       for (k in c(1:nz)) {
    197         m1 <- min(i+1, nx)
    198         m2 <- min(j+1, ny)
    199         m3 <- min(k+1, nz)
    200         for (l in c(i:m1)) {
    201           for (m in c(j:m2)) {
    202             for (n in c(k:m3)) {
    203               v1 <- c(i,j,k)
    204               v2 <- c(l,m,n)
    205               if (!(i==l && j==m && k==n)) {
    206                 arma.arrow(v1, v2, part_size)
    207               }
    208             }
    209           }
    210         }
    211       }
    212     }
    213   }
    214   # rotate the camera
    215   view3d(45,30)
    216 }
    217 
    218 arma.plot_ar_cubes_2d_v2 <- function (nx, ny, xlabel, ylabel, args) {
    219   if (!('arrow_args' %in% names(args))) {
    220     args$arrow_args <- list(lwd=3,angle=7)
    221   }
    222   if (!('adj_x' %in% names(args))) {
    223     args$adj_x <- 0.3
    224   }
    225   if (!('adj_y' %in% names(args))) {
    226     args$adj_y <- 0.4
    227   }
    228   part_size <- 2
    229   plot.new()
    230   plot.window(xlim=c(0,nx*part_size),ylim=rev(c(0,ny*part_size)))
    231   par(pty="s")
    232   char_idx <- 1
    233   adj_x <- args$adj_x
    234   adj_y <- args$adj_y
    235   for (i in c(1:nx)-1) {
    236     for (j in c(1:ny)-1) {
    237       x0 <- i*part_size
    238       y0 <- j*part_size
    239       x1 <- x0 + part_size
    240       y1 <- y0 + part_size
    241       rect(x0, y0, x1, y1)
    242       text(x0+part_size/2, y0+part_size/2, LETTERS[char_idx], cex=1.5)
    243       m1 <- max(i-1, 0)
    244       m2 <- max(j-1, 0)
    245       # draw arrows denoting AR dependencies
    246       for (l in c(m1:i)) {
    247         for (m in c(m2:j)) {
    248           if (!(l==i && m==j)) {
    249             xl <- l*part_size
    250             ym <- m*part_size
    251             x00 <- x0
    252             y00 <- y0
    253             if (xl != x0) {
    254               xl <- xl + adj_x
    255               x00 <- x00 - adj_x
    256             }
    257             if (ym != y0) {
    258               ym <- ym + adj_y
    259               y00 <- y00 - adj_y
    260             }
    261             do.call(arrows, c(list(
    262               x0=x00 + part_size/2,
    263               y0=y00 + part_size/2,
    264               x1=xl + part_size/2,
    265               y1=ym + part_size/2),
    266               args$arrow_args))
    267           }
    268         }
    269       }
    270       char_idx <- char_idx + 1
    271     }
    272   }
    273   if (!('no_axes' %in% names(args)) | !args$no_axes) {
    274     axis(3, at=c(0:nx)*part_size, labels=c(0:nx))
    275     axis(2, at=c(0:ny)*part_size, labels=c(0:ny))
    276     mtext(xlabel, side=3, line=3)
    277     mtext(ylabel, side=2, line=3)
    278   }
    279 }
    280 
    281 arma.plot_ar_cubes_2d <- function (nx, ny, xlabel, ylabel) {
    282   arma.plot_ar_cubes_2d_v2(nx, ny, xlabel, ylabel, list())
    283 }
    284 
    285 arma.plot_factory_vs_openmp <- function(...) {
    286   args <- list(...)
    287   perf <- read.csv(file.path("data", "performance", "factory-vs-openmp.csv"))
    288   scale <- 10 ** args$power
    289   x <- perf$nt * perf$nx * perf$ny / scale
    290   plot.new()
    291   plot.window(xlim=range(x),ylim=range(perf[c("openmp", "factory")]))
    292   pts <- pretty(x)
    293   axis(1, at=pts, labels=sapply(pts, function(x) {as.expression(bquote(.(x) %.% 10 ** .(args$power)))}))
    294   axis(2)
    295   box()
    296   lines(x, perf$openmp, lty="solid")
    297   lines(x, perf$factory, lty="dashed")
    298   title(xlab=args$xlab, ylab=args$ylab)
    299 }
    300 
    301 arma.plot_factory_vs_openmp_overlap <- function(...) {
    302   args <- list(...)
    303   openmp <- read.csv(file.path("data", "performance", "overlap-openmp.csv"), na.strings="")
    304   factory <- read.csv(file.path("data", "performance", "overlap-factory.csv"), na.strings="")
    305   openmp$t <- (openmp$t - min(openmp$t)) / args$scale
    306   factory$t <- (factory$t - min(factory$t)) / args$scale
    307   plot.new()
    308   plot.window(xlim=range(c(factory$t, openmp$t)),ylim=range(0, 5))
    309   axis(1)
    310   axis(2, at=c(1, 3), labels=args$labels, las=1, hadj=1)
    311   # OpenMP
    312   lines(openmp$t, rep.int(3, length(openmp$t)))
    313   openmp_pts <- openmp[!is.na(openmp$mark),]
    314   openmp_y <- rep.int(3, length(openmp_pts$t))
    315   points(openmp_pts$t, openmp_y)
    316   text(openmp_pts$t, openmp_y, labels=openmp_pts$mark, pos=c(3, 3, 1, 1))
    317   # Factory
    318   lines(factory$t, rep.int(1, length(factory$t)))
    319   factory_pts <- factory[!is.na(factory$mark),]
    320   factory_y <- rep.int(1, length(factory_pts$t))
    321   points(factory_pts$t, factory_y)
    322   text(factory_pts$t, factory_y, labels=factory_pts$mark, pos=c(3, 1, 3, 1))
    323   title(xlab=args$xlab)
    324 }
    325 
    326 # a workaround for a bug in ascii package
    327 # which does not honour "decimal.mark" argument
    328 arma.print_ascii <- function (obj) {
    329 	decimal.mark <- options("arma.mark")
    330 	if (is.null(decimal.mark)) {
    331 		decimal.mark <- "."
    332 	}
    333 	replacement <- paste('\\1', decimal.mark, '\\2', sep='')
    334 	cat(gsub('([0-9]+)\\.([0-9]+)', replacement, capture.output(obj$show.org())), sep="\n")
    335 }