quad_form_test.R
In emulator: Bayesian Emulation of Computer Programs

## ----set-options, echo = FALSE------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", dev = "png", fig.width = 7, fig.height = 3.5, message = FALSE, warning = FALSE)
options(width = 80, tibble.width = Inf)

## -----------------------------------------------------------------------------
library("emulator")


## -----------------------------------------------------------------------------
rcm <- function(row,col){matrix(rnorm(row*col)+1i*rnorm(row*col),row,col)}
M <- rcm(2,2)
x <- rcm(2,3)
y <- rcm(3,2)
x1 <- rcm(2,3)
y1 <- rcm(3,2)

## -----------------------------------------------------------------------------
tester <- function(a,b,TOL=1e-13){stopifnot(all(abs(a-b)< TOL))}

## -----------------------------------------------------------------------------
(jj1 <- Conj(t(x)))
(jj2 <- t(Conj(x)))
(jj3 <- ht(x))
tester(jj1,jj3)
tester(jj2,jj3)

## -----------------------------------------------------------------------------
(jj1 <- ht(x) %*% x1)
(jj2 <- cprod(x,x1))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- ht(x1) %*% x)
(jj2 <- cprod(x1,x))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- ht(x) %*% M %*% x)
(jj2 <- quad.form(M,x))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- ht(x) %*% solve(M) %*% x)
(jj2 <- quad.form(solve(M),x))
max(abs(jj1-jj2))

## -----------------------------------------------------------------------------
(jj1 <- ht(x) %*% M %*% x1)
(jj2 <- quad.3form(M,x,x1))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- y %*% M %*% ht(y1))
(jj2 <- quad.3tform(M,y,y1))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- y %*% M %*% ht(y))
(jj2 <- quad.tform(M,y))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- y %*% solve(M) %*% ht(y))
(jj2 <- quad.tform.inv(M,y))
tester(jj1,jj2)

## -----------------------------------------------------------------------------
(jj1 <- diag(ht(x) %*% M %*% x))
(jj2 <- diag(quad.form(M,x)))
(jj3 <- quad.diag(M,x))
tester(jj1,jj3)
tester(jj2,jj3)

## -----------------------------------------------------------------------------
(jj1 <- diag(y %*% M %*% ht(y)))
(jj2 <- diag(quad.tform(M,y)))
(jj3 <- quad.tdiag(M,y))
tester(jj1,jj3)
tester(jj2,jj3)

## -----------------------------------------------------------------------------
(jj1 <- diag(ht(x) %*% M %*% x1))
(jj2 <- diag(quad.3form(M,x,x1)))
(jj3 <- quad.3diag(M,x,x1))
tester(jj1,jj3)
tester(jj2,jj3)

## -----------------------------------------------------------------------------
(jj1 <- diag(y %*% M %*% ht(y1)))
(jj2 <- diag(quad.3tform(M,y,y1)))
(jj3 <- quad.3tdiag(M,y,y1))
tester(jj1,jj3)
tester(jj2,jj3)