library(MASS)

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

ntw.names <- names(ntw.lrn)

y <- ntw.lrn$connection.type

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

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"]))

real <- c(0.1,0.8,0.1)

fit <- lda(cbind(f23,f29),y,prior=real) 
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,2)
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];  
  }
}

pred <- predict(fit,grid,prior=real)

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_real.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()