library(MASS)

ntw.lrn <- read.table("../lrn/ntw_lrn.dat",header=T,sep=",")

ntw.names <- names(ntw.lrn)

y.lrn <- ntw.lrn$connection.type

f23 <- ntw.lrn$count
f29 <- ntw.lrn$same.srv.rate
f.lrn <- cbind(f23,f29,f23^2,f29^2,f23*f29)

rm(ntw.lrn)

m.neptune <- cbind(mean(f23[y=="neptune"]),mean(f29[y=="neptune"]))
m.normal <- cbind(mean(f23[y=="normal"]),mean(f29[y=="normal"]))
m.smurf <- cbind(mean(f23[y=="smurf"]),mean(f29[y=="smurf"]))

fit <- lda(f.lrn,y.lrn) 
print(fit$prior) 

size <- 100

f23.axis <- seq(from=min(f23),to=max(f23),length=100)
f29.axis <- seq(from=min(f29),to=max(f29),length=100)

grid <- mat.or.vec(size*size,5)
for (i in 1:size) {
  for (j in 1:size) {
    grid[size*(j-1)+i,1] <- f23.axis[i]
    grid[size*(j-1)+i,2] <- f29.axis[j]  
    grid[size*(j-1)+i,3] <- f23.axis[i]^2
    grid[size*(j-1)+i,4] <- f29.axis[j]^2  
    grid[size*(j-1)+i,5] <- f29.axis[j]*f23.axis[i]
  }
}

pred <- predict(fit,grid)

yhat <- pred$class

edge <- mat.or.vec(size*size,1)
edge <- as.logical(edge)
for (i in 1:size) { 
  for (j in 1:size) {
    cntr <- size*(j-1)+i
    left <- cntr - size
    down <- cntr - 1
    if ( (left > 0) && (down > 0) ) {
      edge[cntr] <- (yhat[left] != yhat[cntr]) || (yhat[down] != yhat[cntr])
    }
    if (i == 1) edge[cntr] <- FALSE
  }
}

source("psopts.r")
postscript(file="ntw_lda_driv.eps")
par.default <- par(no.readonly=TRUE) 
plot(grid,type='n',xlab=ntw.names[23],ylab=ntw.names[29])
  points(grid[yhat=="neptune",],pch='.',col="red")
  points(grid[yhat=="normal",],pch='.',col="green")
  points(grid[yhat=="smurf",],pch='.',col="blue")
  points(m.neptune,pch="o",col="red")
  points(m.normal,pch="o",col="green")
  points(m.smurf,pch="o",col="blue")
  points(grid[edge,],pch='+',col="black",cex=0.4)
dev.off()

ntw.val <- read.table("../val/ntw_val.dat",header=T,sep=",")

y.val <- ntw.val$connection.type

f23 <- ntw.val$count
f29 <- ntw.val$same.srv.rate
f.val <- cbind(f23,f29,f23^2,f29^2,f23*f29)

yhat <- predict(fit,f.val)
n <- length(y.val)
ehat <- (y.val != yhat$class)
cer <- sum(ehat)/n

print(cer)