#-*- R -*-

# generation of populations predicted from data

library(MASS)
library(class)
library(nnet)

#decrease len if you have little memory.

decplot <- function(xp, yp, Z.a, Z.b, Z.c, t)
{
    plot(Pop[, 1], Pop[, 2], xlim=c(-6.0,10), ylim=c(-6,10), type = "n", xlab = "x1", ylab = "x2")
	title(t)
    for(il in 1:3) {
        set <- Pop$y==levels(Pop$y)[il]
        text(Pop[set, 1], Pop[set, 2],
             labels = as.character(Pop$y[set]), col = 2 + il) }
    zp <- Z.a - pmax(Z.b, Z.c)
    contour(xp, yp, matrix(zp, np), add = T, levels = 0, labex = 0)
    zp <- Z.c - pmax(Z.a, Z.b)
    contour(xp, yp, matrix(zp, np), add = T, levels = 0, labex = 0)
    invisible()
}

# total number of datapoints

N <- 500

# total number of iterations

j <- 1

# read in data - with class information

c1<-read.table("ex20a.dat")
names(c1)<-c("x1","x2")

c2<-read.table("ex20b.dat")
names(c2)<-c("x1","x2")

c3<-read.table("ex20c.dat")
names(c3)<-c("x1","x2")

# compute covariance matrices (without class info)

c1.cov<-cov(cbind(c1$x1,c1$x2))
c2.cov<-cov(cbind(c2$x1,c2$x2))
c3.cov<-cov(cbind(c3$x1,c3$x2))

# compute the means

c1.m <- diag(c(mean(c1$x1),mean(c1$x2)))
c2.m <- diag(c(mean(c2$x1),mean(c2$x2)))
c3.m <- diag(c(mean(c3$x1),mean(c3$x2)))

# set the covariance and mean
ones <- matrix( c(rep(1,N), rep(1,N)), ncol=2)

# rotate the data
c1.r <- chol(c1.cov)
c2.r <- chol(c2.cov)
c3.r <- chol(c3.cov)

# set up ploting parameters 

if( j > 6 )
	j<- 1

if( j == 1)
	par(mfcol=c(1,1)) 

if( j == 2 )
	par(mfcol=c(2,1))

if( j == 3 )  
	par(mfcol=c(3,1))

if( j == 4 )
	par(mfcol=c(2,2))

if( j == 5 )
	par(mfcol=c(2,3))

if( j == 6 )
	par(mfcol=c(2,3))

# loop through process i times

for(i in 1:j)
{
	# generate the data: normal distributions

	pop1 <- data.frame(rnorm(N),rnorm(N))
	names(pop1)<-c("x1","x2")
	pop2 <- data.frame(rnorm(N),rnorm(N))
	names(pop2)<-c("x1","x2")
	pop3 <- data.frame(rnorm(N),rnorm(N))
	names(pop3)<-c("x1","x2")

	pop1.r <- as.matrix(pop1) %*% c1.r
	pop2.r <- as.matrix(pop2) %*% c2.r
	pop3.r <- as.matrix(pop3) %*% c3.r

	# shift the data
	pop1.r.s <- pop1.r + (ones %*% c1.m)
	pop2.r.s <- pop2.r + (ones %*% c2.m)
	pop3.r.s <- pop3.r + (ones %*% c3.m)

	# add in the class information

	pop1.r.s.c = data.frame(pop1.r.s, rep("a",length(pop1.r.s[,1])))
	names(pop1.r.s.c)<-c("x1","x2","y")
	pop2.r.s.c = data.frame(pop2.r.s, rep("b",length(pop2.r.s[,1])))
	names(pop2.r.s.c)<-c("x1","x2","y")
	pop3.r.s.c = data.frame(pop3.r.s, rep("c",length(pop3.r.s[,1])))
	names(pop3.r.s.c)<-c("x1","x2","y")

	Pop <- data.frame(rbind(pop1.r.s.c, pop2.r.s.c, pop3.r.s.c))

	# Pop[, 3] would return just the third column
	# Pop[,-3] returns everything but the third column

	p <- data.frame(Pop$x1,Pop$x2,class.ind(Pop$y))
	names(p)<-c("x1","x2","a","b","c")

	xp <- seq(-6.0, 10.0, length = 100); np <- length(xp)
	yp <- seq(-6.0, 10.0, length = 100)

	pt <- expand.grid(x1 = xp, x2 = yp)

	g.a <- lm( a ~ x1 + x2, data = p)
	g.b <- lm( b ~ x1 + x2, data = p)
	g.c <- lm( c ~ x1 + x2, data = p)

	Z.a <- predict(g.a, pt)
	Z.b <- predict(g.b, pt)
	Z.c <- predict(g.c, pt)

	decplot(xp, yp, Z.a, Z.b, Z.c, "Regression of Response Matrix")

}