#'---
#'title: "Cube: 3D phase plane analysis in Grind"
#'author: "Rob de Boer"
#'output:
#'  html_document: default
#'---
#'
#' This script `cube.R` allows for 3D phase plane analysis, by defining the
#' functions `cube()` and `run3d()` that both use the `plot3D` library.
#' The call to `cube()` is similar to the call to `plane()`. 
#' The `cube.R` script can be downloaded from the Grind
#' [directory](https://tbb.bio.uu.nl/rdb/grindR/).
#' There are **three extra options**:  
#' 1. `shade="xyx"` allows one to shade nullclines in the `x`, `y`, and/or `z` direction.
#' The shading is off by default (`shade=""`).  
#' 2. `box=FALSE` suppresses the drawing of the upper edges of the cube.  
#' 3. `grid=5` defines that 5 lines are used in the shading.
#'
#' The viewpoint is defined by the defaults `theta=40` and `phi=40`,
#' where `theta` gives the azimuthal direction and `phi` the colatitude.
#' One can use `plotdev(theta=60,phi=40)` to change the viewpoint afterwards.
#' It is typically good top set `par(mar=c(0,0,0,0))` to reduce the margins.
#' 
#' The tutorial [chaos3d](https://tbb.bio.uu.nl/rdb/grindR/chaos3d.html)
#' provides and example of calling `cube()` and `plot3d()`.
#'

#' After also installing the `rgl` and the `plot3Drgl` R-packages, 
#' one can call the `plotrgl()` function to open a new window
#' in which the graphics can be rotated and zoomed with the mouse 
#' (try using the left and right buttons).
#' 
 
#' Finally, if one prefers to have the origin in the left hand back corner, 
#' with `x` pointing rightwards, `y` upwards, and `z` forwards, 
#' call `cube()` and `run3d()` with `(x=3, y=1, z=2, theta=130, ...)`.
#' 

cube_version <- "15-05-2024"
library(plot3D)

x_cube <- 1; xmin_cube <- 0.001; xmax_cube <- 1.05
y_cube <- 2; ymin_cube <- 0.001; ymax_cube <- 1.05
z_cube <- 3; zmin_cube <- 0.001; zmax_cube <- 1.05
log_cube <- ""

cube <- function(x=1, y=2, z=3, xmin=0.001, xmax=1.05, ymin=0.001, ymax=1.05, zmin=0.001, zmax=1.05, xlab="", ylab="", zlab="", log="", shade="", grid=5, npixels=150, state=s, parms=p, odes=model, time=0, eps=NULL, show=NULL, zero=TRUE, add=FALSE, box=FALSE, lwd=1, ...) {
  # Make a 3D-phase space with nullclines
  dots <- list(...)
  if (!is.null(dots)) {
    unknown <- names(dots[!names(dots) %in% c(args_run,args_plot,names(formals(contour3D)),"ticktype","nticks")])
    if (length(unknown)>0) warning(paste("Unknown argument(s):",unknown,sep=" "))
  }
  dots_run <- if (!is.null(dots)) dots[names(dots) %in% args_run] else NULL
  if (!is.null(eps)) print("Option eps is no longer supported")
  if (add) {
    x <- x_cube
    y <- y_cube
    z <- z_cube
    xmin <- xmin_cube; xmax <- xmax_cube
    ymin <- ymin_cube; ymax <- ymax_cube
    zmin <- zmin_cube; zmax <- zmax_cube
    log <- log_cube
  } else {
    if (!is.numeric(x)) x <- index(x,names(state))
    if (!is.numeric(y)) y <- index(y,names(state))
    if (!is.numeric(z)) z <- index(z,names(state))
    x_cube <<- x
    y_cube <<- y
    z_cube <<- z
    xmin_cube <<- xmin; xmax_cube <<- xmax
    ymin_cube <<- ymin; ymax_cube <<- ymax 
    zmin_cube <<- zmin; zmax_cube <<- zmax
    log_cube <<- log
  }
  ishows <- if (!is.null(show)) index(show, names(state)) else c(x, y, z)
  nvar <- length(state)
  if (zero) state[1:nvar] <- rep(0,nvar)
  logx <- ifelse(grepl('x',log), TRUE, FALSE)
  logy <- ifelse(grepl('y',log), TRUE, FALSE)
  logz <- ifelse(grepl('z',log), TRUE, FALSE)
  xc <- plane_coord(logx,xmin,xmax,npixels)
  yc <- plane_coord(logy,ymin,ymax,npixels)
  zc <- plane_coord(logz,zmin,zmax,npixels)
  xvar <- names(state)[x]; yvar <- names(state)[y]; zvar <- names(state)[z]
  if (xlab == "") xlab <- names(state)[x]
  if (ylab == "") ylab <- names(state)[y]
  if (zlab == "") zlab <- names(state)[z]
  npixels2 <- npixels^2
  vparms <- as.list(parms)
  vparms <- lapply(vparms,rep,vparms,npixels2)
  vstate <- as.list(state)
  #vstate<-lapply(vstate,rep,vstate,npixels2);vstate[[x]]<-0;vstate[[y]]<-0
  for (j in seq(1,nvar)) if (j!=x & j!=y & j!=z) vstate[[j]]<-rep(vstate[[j]],npixels2);
  
  #xzy
  dxzy <- list()
  nxzy <- ifelse(grepl('y',shade), grid+2, 2)
  pxzy <- plane_coord(logy,ymin,ymax,nxzy)
  vstate[[x]] <- rep.int(xc, npixels)
  vstate[[z]] <- rep.int(zc, rep.int(npixels, npixels))  #outer(xc,zc)
  for (i in seq(nxzy)) {
    vstate[[y]] <- rep.int(pxzy[i], npixels2)
    dxzy[[i]] <- odes(time,vstate,vparms)[[1]]
    dim(dxzy[[i]]) <- c(npixels,npixels,nvar)
  }
  
  #xyz
  dxyz <- list()
  nxyz <- ifelse(grepl('z',shade), grid+2, 2)
  pxyz <- plane_coord(logz,zmin,zmax,nxyz)
  vstate[[y]] <- rep.int(yc, rep.int(npixels, npixels))
  for (i in seq(nxyz)) {
    vstate[[z]] <- rep.int(pxyz[i], npixels2)
    dxyz[[i]] <- odes(time,vstate,vparms)[[1]]
    dim(dxyz[[i]]) <- c(npixels,npixels,nvar)
  }
  
  #yzx
  dyzx <- list()
  nyzx <- ifelse(grepl('x',shade), grid+2, 2)
  pyzx <- plane_coord(logx,xmin,xmax,nyzx)
  vstate[[y]] <- rep.int(yc, npixels)
  vstate[[z]] <- rep.int(zc, rep.int(npixels, npixels))
  for (i in seq(nyzx)) {
    vstate[[x]] <- rep.int(pyzx[i], npixels2)
    dyzx[[i]] <- odes(time,vstate,vparms)[[1]]
    dim(dyzx[[i]]) <- c(npixels,npixels,nvar)
  }
  
  addcon <- add
  for (i in ishows) {
    try(contour3D(x=xc,y=pxzy[1],z=zc,colvar=dxzy[[1]][,,i],levels=0,col=colors[i],add=addcon,addbox=box,lwd=lwd,xlim=c(xmin,xmax),ylim=c(ymin,ymax),zlim=c(zmin,zmax),xlab=names(state[x]),ylab=names(state[y]),zlab=names(state[z]),...),silent=TRUE)
    addcon <- TRUE
    for (j in seq(2,nxzy-1))
      try(contour3D(x=xc,y=pxzy[j],z=zc,colvar=dxzy[[j]][,,i],levels=0,col=colors[i],add=TRUE,addbox=FALSE,lwd=lwd,...),silent=TRUE)
    try(contour3D(x=xc,y=pxzy[nxzy],z=zc,colvar=dxzy[[nxzy]][,,i],levels=0,col=colors[i],add=TRUE,addbox=box,lwd=lwd,...),silent=TRUE)
    
    try(contour3D(x=xc,y=yc,z=pxyz[1],colvar=dxyz[[1]][,,i],levels=0,col=colors[i],add=TRUE,addbox=box,lwd=lwd,...),silent=TRUE)
    for (j in seq(2,nxyz-1))
      try(contour3D(x=xc,y=yc,z=pxyz[j],colvar=dxyz[[j]][,,i],levels=0,col=colors[i],add=TRUE,addbox=FALSE,lwd=lwd,...),silent=TRUE)
    try(contour3D(x=xc,y=yc,z=pxyz[nxyz],colvar=dxyz[[nxyz]][,,i],levels=0,col=colors[i],add=TRUE,addbox=box,lwd=lwd,...),silent=TRUE)
    
    try(contour3D(x=pyzx[1],y=yc,z=zc,colvar=dyzx[[1]][,,i],levels=0,col=colors[i],add=TRUE,addbox=box,lwd=lwd,...),silent=TRUE)
    for (j in seq(2,nyzx-1))
      try(contour3D(x=pyzx[j],y=yc,z=zc,colvar=dyzx[[j]][,,i],levels=0,col=colors[i],add=TRUE,addbox=FALSE,lwd=lwd,...),silent=TRUE)
    try(contour3D(x=pyzx[nyzx],y=yc,z=zc,colvar=dyzx[[nyzx]][,,i],levels=0,col=colors[i],add=TRUE,addbox=box,lwd=lwd,...),silent=TRUE)
  }
}

run3d <- function(x=1, y=2, z=3, xmin=0.001, xmax=1.05, ymin=0.001, ymax=1.05, zmin=0.001, zmax=1.05, log="", col=1, add=FALSE, state=s, parms=p, ...) {
  dots <- list(...)
  if (!is.null(dots)) {
    args_run3d <- c(names(formals(scatter3D)),"ticktype","nticks")
    unknown <- names(dots[!names(dots) %in% c(args_run,args_run3d)])
    if (length(unknown)>0) warning(paste("Unknown argument(s):",unknown,sep=" "))
    dots_run <- dots[names(dots) %in% args_run]
    dots_l3d <- dots[names(dots) %in% args_run3d]
  }else dots_run <- NULL
  if (add) {
    x <- x_cube
    y <- y_cube
    z <- z_cube
    xmin <- xmin_cube; xmax <- xmax_cube
    ymin <- ymin_cube; ymax <- ymax_cube
    zmin <- zmin_cube; zmax <- zmax_cube
    log <- log_cube
  } else {
    if (!is.numeric(x)) x <- index(x,names(state))
    if (!is.numeric(y)) y <- index(y,names(state))
    if (!is.numeric(z)) z <- index(z,names(state))
    x_cube <<- x
    y_cube <<- y
    z_cube <<- z
    xmin_cube <<- xmin; xmax_cube <<- xmax
    ymin_cube <<- ymin; ymax_cube <<- ymax 
    zmin_cube <<- zmin; zmax_cube <<- zmax
    log_cube <<- log
  }
  data <- do.call('run',c(list(state=state,parms=parms,timeplot=FALSE,table=TRUE),dots_run))
  if (!add) {
    do.call('lines3D',c(list(x=data[,1+x],y=data[,1+y],z=data[,1+z],xlim=c(xmin,xmax),ylim=c(ymin,ymax),zlim=c(zmin,zmax),xlab=names(state[x]),ylab=names(state[y]),zlab=names(state[z]),col=col),dots_l3d))
  } else {
    do.call('lines3D',c(list(x=data[,1+x],y=data[,1+y],z=data[,1+z],col=col,add=TRUE),dots_l3d))
  }
  f <- state
  f[1:length(f)] <- as.numeric(data[nrow(data),2:(length(s)+1)])
  return(f)
}

cat(paste("cube.R (",cube_version,") was sourced\n",sep=""))