######################################################
#
# R implementation of a simple perceptron with
# visualization in 2D corresponding to exercise 1.2 of
# http://www.dbs.ifi.lmu.de/Lehre/MaschLernen/SS2009/index.html
#
# programming language R available at
#    http://cran.r-project.org/
#
# available in the CIP pool by
# $ R
#
#
######################################################

######################################################
#
# x   <- miss-classified pattern to be inserted
# y   <- class labels for pattern x
# w   <- vector of weights
# eta <- step size for each correction of a miss-
#         classified pattern
#
# performs one step of a pattern-based descent
#
######################################################
pStep<-function(x,y,w,eta=.1) {
  stopifnot(length(x) == length(w)-1);
  w[1] = w[1] + eta * y;
  for (j in 2:length(w)) {
    w[j] = w[j] + eta * y * x[j-1];
  }
  w;
}

######################################################
#
# w   <- vector of weights = separating hyperplane
# col <- color of straight line
# lwd <- thickness of straight line
#
# plots a seperating hyperplane w for two dimensions
#
######################################################
plotMe<-function(w,col=1,lwd=1) {
  if(length(w) == 3) {
    if (w[3]==0) {
      if (w[2]==0) {
        abline(v=0,col=col,lwd=lwd);
      } else {
        abline(v=-w[1]/w[2],col=col,lwd=lwd);
      }
    } else if (is.infinite(w[2])) {
      abline(v=0,col=col,lwd=lwd);
    } else if (is.infinite(w[3])) {
      abline(h=0,col=col,lwd=lwd);
    } else {
      abline(-w[1]/w[3],-w[2]/w[3],col=col,lwd=lwd);
    }
  }
}

######################################################
#
# m   <- matrix of feature vectors (stored in columns)
# eta <- step size for each correction of a miss-
#         classified pattern
#
# header for Plots resulting from runP() and gradD()
#
######################################################
plotHeader<-function(m,eta=0.1) {
  plot(m[1,],m[2,], xlab=expression(x[1]),
    ylab=expression(x[2]), xlim=range(-1,5),
    ylim=range(-1,5), col=ifelse(y==-1,2,3),
    main="Perceptron Iterations",
    sub=substitute(list(eta) == x, list(x=eta)))
  text(m[1,],m[2,], 1:dim(m)[2],
    col=ifelse(y==-1,2,3), pos=4)
  abline(h=0,lty=3)
  abline(v=0,lty=3)
}

######################################################
#
# m   <- matrix of feature vectors (stored in columns)
# y   <- vector of class labels for feature vectors;
#         these must be in {-1, 1}
# w   <- initialization vector of weights
# eta <- step size for each correction of a miss-
#         classified pattern
# ord <- order of repeated "random" pattern sampling
#         procedure; if == -1, random pattern selection
#
# pattern-based Perceptron implementation
#
######################################################
runP<-function(m,y,w=c(0,1,-1),eta=.1,ord=1:4) {
  plotMe(w,8,3)
  
  while (1) { # until convergence
    wtemp=w;
    for (i in ord) {
      if (i == -1) # random pattern
        i = sample(1:(dim(m)[2]),1)
      h<-sum(c(1,m[,i])*w);
      if (sign(h) == y[i]) {
        next; # correctly classified
      }
      w<-pStep(m[,i],y[i],w,eta)
      
      print (sprintf("%2d : % 3.3f,% 3.3f,% 3.3f",i,
              w[1],w[2],w[3]))
      k <- k+1;
      plotMe(w,cols[k %% length(cols)])
    }
    if (isTRUE(all.equal(wtemp,w))) {
      break;
    }
  }
  
  print(paste(k,"iterations"))
  plotMe(w,1,3)
  w
}

######################################################
#
# Gradient Descend variant of runP()
#
# result displays the number of iterations
#
######################################################
gradD<-function(m,y,w=c(0,1,-1),nu=.1) {
  plotMe(w,8,3)
  nullers <- rep(0,length(w));
  numIt<-0;
  while (1) {
    change<-rep(0,length(w));
    for (i in 1:length(y)) {
      h<-sum(c(1,m[,i])*w);
      if (sign(h) == y[i]) {
        next;
      }
      change[1] = change[1] + y[i];
      for (j in 2:length(w)) {
        change[j] = change[j] + y[i] * m[j-1,i];
      }
      k <- k+1;
    }
    if (isTRUE(all.equal(change,nullers))) {
      break;
    }
    numIt <- numIt+1;
    w <- w + nu * change;
    print (sprintf("Iter %2d - %3d : % 3.3f,% 3.3f,% 3.3f",
                    numIt,k,w[1],w[2],w[3]))
    plotMe(w,cols[numIt %% length(cols)])
  }
  print(paste(numIt,"iterations,",k,"used patterns"))
  plotMe(w,1,3)
  w
}

######################################################
#
# examplary function calls:
#
######################################################

eta<-.1
k<-0
cols<-rainbow(100) # colors for visualization
w<-c(0,1,-1)
m<-matrix(c(2,4, 1,.5, .5,1.5, 0,.5),nrow=2)
y<-c(1,1,-1,-1)


# 1-2 a)
plotHeader(m,eta)
runP(m,y,w,eta)

# output:
# [1] " 1 :  0.100, 1.200,-0.600"
# [1] 1 iterations
# [1]  0.1  1.2 -0.6
#
# 1-2 b)
#
# order {1,3,3,1,3,3}
eta<-.25
cols=rainbow(5)
cols=c(cols,1)
plotHeader(m,eta)
runP(m,y,w,eta)
legend("bottomright",legend=1:6,col=cols,lwd=c(rep(1,5),3))

# output:
# [1] " 1 :  0.250, 1.500, 0.000"
# [1] " 3 :  0.000, 1.375,-0.375"
# [1] " 3 : -0.250, 1.250,-0.750"
# [1] " 1 :  0.000, 1.750, 0.250"
# [1] " 3 : -0.250, 1.625,-0.125"
# [1] " 3 : -0.500, 1.500,-0.500"
# [1] 6 iterations
# [1] -0.5  1.5 -0.5

#
# order {1,4,3,2}
cols=rainbow(2)
cols=c(cols,1)
plotHeader(m,eta)
runP(m,y,w,.25, c(1,4,3,2))
legend("bottomright",legend=1:3,col=cols,lwd=c(rep(1,2),3))

# output:
# [1] " 1 :  0.250, 1.500, 0.000"
# [1] " 4 :  0.000, 1.500,-0.125"
# [1] " 3 : -0.250, 1.375,-0.500"
# [1] 3 iterations
# [1] -0.250  1.375 -0.500

# 1-2 c)
#
# eta .1
eta=.1
plotHeader(m,eta)
gradD(m,y,w,eta)

# output:
# [1] "Iter  1 -   1 :  0.100, 1.200,-0.600"
# [1] "1 iterations, 1 used patterns"
# [1]  0.1  1.2 -0.6

# eta .25
eta=.25
plotHeader(m,eta)
gradD(m,y,w,eta)

# output:
# [1] "Iter  1 -   1 :  0.250, 1.500, 0.000"
# [1] "Iter  2 -   3 : -0.250, 1.375,-0.500"
# [1] "2 iterations, 3 used patterns"
# [1] -0.250  1.375 -0.500