#
# Contour plots of quadratic functions
#
# Code example for lecture: Numerical Optimization
# Author: Sangkyun Lee (sangkyun.lee@tu-dortmund.de)
#

library(rgl)

M = NULL
M[[1]] <- matrix(c(10,0,0,10), byrow=TRUE, nrow=2, ncol=2) # convex, axis-aligned
M[[2]] <- matrix(c(10,0,0,5), byrow=TRUE, nrow=2, ncol=2) # convex, axis-aligned
M[[3]] <- matrix(c(10,0,0,1), byrow=TRUE, nrow=2, ncol=2) # convex, axis-aligned, ill-conditioned
M[[4]] <- matrix(c(10,3,3,5), byrow=TRUE, nrow=2, ncol=2) # convex, not aligned to axes
M[[5]] <- matrix(c(-10,0,0,-5), byrow=TRUE, nrow=2, ncol=2) # concave
M[[6]] <- matrix(c(10,0,0,-10), byrow=TRUE, nrow=2, ncol=2) # neither convex nor concave

#
# Generate a 3D contour plot of the quadratic function f=(x,y)^T M (x,y), for a symmetric 2x2 matrix M.
# Also shows line segments representing eigenvectors, of lengths proportional to eigenvalues.
# The line segments are color-coded: blue (positive eigenvalue), cyan (negative)
#
plotQuad <- function(M, x, y) {
  
  f <- function(x,y) {
    z <- matrix(c(x,y), byrow=TRUE, nrow=2);
    return(diag(t(z) %*% M %*% z))
  }
  z <- outer(x, y, f)  
  
  #
  # This funcion has been adapted from http://statmath.wu.ac.at/research/talks/resources/tscu3dg.R
  #
  showsurface <- function(x, y, z) {
    persp3d(x,y,z, col="red", alpha=0.4, axes=TRUE)
    contours <- contourLines(x,y,z)
    for (i in 1:length(contours)) 
      with(contours[[i]], lines3d(x, y, level, col="darkred"))
  }
  open3d()
  showsurface(x,y,z)
  
  d <- eigen(M)
  w <- 2* abs(d$values) / sum(abs(d$values))
  lines3d(w[1]*c(min(x),max(x))*d$vectors[1,1], w[1]*c(min(y),max(y))*d$vectors[2,1], c(0,0), col=ifelse(is.complex(d$value[1]), "brown", ifelse(d$value[1]>0, "blue", "cyan")), lwd=2)
  lines3d(w[2]*c(min(x),max(x))*d$vectors[1,2], w[2]*c(min(y),max(y))*d$vectors[2,2], c(0,0), col=ifelse(is.complex(d$value[2]), "brown", ifelse(d$value[2]>0, "blue", "cyan")), lwd=2)
}

x <- seq(-1,1,len=40)
y <- seq(-1,1,len=40)

for(i in seq_len(length(M))) {
  plotQuad(M[[i]], x, y)  
  scan()
}