library(RandomFields)

######################################################################
## here, some examples for the suggestions in practikum.R are given ##
######################################################################


x <- y <- seq(-10, 10, 0.2)

## delay
model <- RMdelay(RMstable(alpha = 1.9, scale = 2), s = c(4, 4))
plot(model, dim = 2, xlim = c(-8, 8))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu, zlim = 'joint')


## schur: lecture by Sudipto Banerjee, p.6
T <- matrix(c(3,1,1,2), ncol=2)
rho <- RMexp(var=5, scale=2)
C <- RMschur(M=T, RMmatrix(M=c(1,1), rho))
plot(C, dim = 2, xlim = c(-3, 3))
simu <- RFsimulate(C, x, y, grid = TRUE)
plot(simu)


## bigneiting
model <- RMbigneiting(kappa=1, mu=1, s=c(6,4), sred12=1,
                      gamma=c(1,2,3), cdiag=c(1,2), rhored=1)
plot(model, dim = 2, xlim = c(-3, 3))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu)


## divergence free
model <- RMdivfree(RMmatern(nu = 5), scale = 4)
plot(model, dim = 2, xlim = c(-3, 3))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu, select.variables = list(1, 2 : 3, 4))


## general vector fields
model <- RMvector(RMgauss(), scale=0.3)
plot(model, dim = 2, xlim = c(-0.5, 0.5))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu, select.variables = list(1:2))


## Kolmogorov's stochastic model of turbulence
x <- y <- seq(-2, 2, len = 20)
model <- RMkolmogorov()
simu <- RFsimulate(model, x, y, z = 0, grid = TRUE)
plot(simu, select.variables = list(1 : 2), col = c("red"))
plot(model, dim = 3, xlim = c(-3, 3), MARGIN = 1 : 2, 
            fixed.MARGIN = 1.0, main = "")


## Multivariate model with anisotropy
M1 <- c(0.9, 0.43)
M2 <- c(0.43, 0.9)
A1 <- RMangle(angle = pi / 4, diag = c(0.1, 0.5))
A2 <- RMangle(angle = 0, diag = c(0.1, 0.5))
x <- y <- seq(-4, 4, 0.2)
model <- RMmatrix(M = M1, RMgengneiting(kappa = 0, mu = 2, Aniso = A1)) +
          RMmatrix(M = M2, RMgengneiting(kappa = 3, mu = 2, Aniso = A2))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu)

## Another model of Anisotropy: space transformation (Guttorp & Sampson, 1992)
Aniso <- RMuser(fctn=x^2, vdim=2)
model <- RMbiwm(nudiag=c(1, 2), nured=1, rhored=0.95, cdiag=c(1, 5),
                 s=c(1, 1, 2), Aniso=Aniso)
x <- y <- seq(-10, 10, 0.2)
simu <- RFsimulate(model, x, y, grid=TRUE)
plot(simu)



## Hierarchical modelling
M1 <- c(0.9, 0.43)
M2 <- c(0.43, 0.9)
model <- RMmatrix(M = M1, RMwhittle(nu = exp())) +
         RMmatrix(M = M2, RMwhittle(nu = 2, scale = 2))
repeat {
  plot(model, dim = 2, xlim = c(-3, 3))
  simu <- RFsimulate(model, x, y, grid = TRUE)
  readline("Press return")
  plot(simu)
  readline("Press return or Ctrl-C")
}



## Extensions of the LMC model
RFoptions(seed = 0, always_close_screen=FALSE)
## compare
M1 <- c(0.9, 0.43)
M2 <- c(0.43, 0.9)
submodel <- RMbiwm(nudiag = c(0.3, 2), nured = 1, rhored = 1, cdiag = c(1, 1),
                     s = c(1, 2, 2))
model <- RMmatrix(M = cbind(M1, M2), submodel)
plot(model, dim = 2, xlim = c(-3, 3))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu)

## with the genuine LMC model
model <- RMmatrix(M = M1, RMwhittle(nu = 0.3)) +
         RMmatrix(M = M2, RMwhittle(nu = 2, scale = 2))
plot(model, dim = 2, xlim = c(-3, 3))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu)


## LMC with 3 components
M1 <- c(0.8, 0.5, 0.3316)
M2 <- c(0.3316, 0.8, 0.5)
M3 <- c(0.5, 0.3316, 0.8)
model <- RMmatrix(M = M1, RMwhittle(nu = 0.3)) +
         RMmatrix(M = M2, RMwhittle(nu = 2, scale = 2)) +
         RMmatrix(M = M3, RMgauss(scale=3))
plot(model, dim = 2, xlim = c(-3, 3))
simu <- RFsimulate(model, x, y, grid = TRUE)
plot(simu)