# $Id: gradasc.R,v 1.1 2007/11/22 04:17:18 kaip Exp $
#
# Copyright (c) 2007 Kai Puolamaki <Kai.Puolamaki@iki.fi>
#
# Permission to use, copy, modify, and distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# copyright notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
#
#
# T-61.3050 Machine Learning: Basic Principles
# For lecture 10/2007
# 


gradasc <- function(x0,func,dfunc,eta=.1,maxt=100) {
# In real life you should not use this function, but some more
# advanced one, see help(optim).

  path <- matrix(NA,nrow=(maxt+1),ncol=length(x0))
  path[1,] <- x0
  x <- x0
  t <- 1
  repeat {
    dx <- dfunc(x)
    x <- x-eta*dx
    cat(sprintf("gradasc: step %d, x=(%.3f,%.3f) f(x)=%.3f\n",t,x[1],x[2],func(x)))
    path[t+1,] <- x
    t <- t+1
    if(t>maxt || sum((eta*dx)^2)<1e-7) break
  }
  list(x=x,path=path)
}

myfunc <- function(x) (x[1]+x[2])^2+(x[2]-1)^2
dmyfunc <- function(x) c(2*(x[1]+x[2]),2*(x[1]+x[2])+2*(x[2]-1)) #gradient

ga1 <- gradasc(c(0,0),myfunc,dmyfunc,eta=.1)
ga2 <- gradasc(c(0,0),myfunc,dmyfunc,eta=.4)



pdf("gradasc.pdf")
x <- seq(-2,0,by=.01)
y <- seq(0,2,by=.01)
z <- matrix(NA,nrow=length(x),ncol=length(y))
for(i in 1:length(x)) {
  for(j in 1:length(y)) {
    z[i,j] <- myfunc(c(x[i],y[j]))
  }
}
plot(c(-2,0),c(0,2),type="n",
     main="Convergence of GRADASC",
     xlab=expression(theta[1]),
     ylab=expression(theta[2]))
contour(x,y,z,add=TRUE)
lines(ga1$path)
points(ga1$path)
lines(ga2$path,col="red")
points(ga2$path,col="red")
points(-1,1,col="green",pch=19)
legend(x="topright",
       legend=c(expression(eta == 0.1),expression(eta == 0.4)),
       lty="solid",
       pch=21,
       col=c("black","red"),
       bg="white")
dev.off()