target.lrn <- read.table("../lrn/num/472.dat",header=T,colClasses="numeric")
target.val <- read.table("../val/num/472.dat",header=T,colClasses="numeric")
 
y.lrn <- target.lrn[,1]
y.val <- target.val[,1]
y     <- c(y.lrn,y.val)

n.lrn <- length(y.lrn)
n.val <- length(y.val)
n     <- length(y)

rm(target.lrn,target.val)

wts <- mat.or.vec(n,1) ; for (i in 1:n.lrn) wts[i]=1

mod <- read.table("cty_mod.txt",
         header=F,colClasses="character",col.names=c("file","feature","type"))

n.mod <- length(mod$file)

first.time <- TRUE

for (i in 1:n.mod) {

  fn.lrn <- paste("../lrn/",mod$type[i],"/",mod$file[i],".dat",sep="")
  fn.val <- paste("../val/",mod$type[i],"/",mod$file[i],".dat",sep="")
  print(mod$feature[i])

  if (mod$type[i]=="chr") {

    f.lrn <- read.table(fn.lrn,
               header=T,colClasses="character",blank.lines.skip=F)
    f.val <- read.table(fn.val,
               header=T,colClasses="character",blank.lines.skip=F)

    f <- c(f.lrn[,1],f.val[,1])

    if (mod$feature[i]=="STATE") {
      f[f=="AS"|f=="DC"|f=="DE"|f=="MA"|f=="ME"|f=="NH"] <- "S1"
      f[f=="OH"|f=="RI"|f=="VI"|f=="WV"]                 <- "S1"
      f[f=="AA"|f=="AE"|f=="AP"|f=="CT"|f=="GU"|f=="MD"] <- "S2"
      f[f=="NJ"|f=="NY"|f=="PA"|f=="PA"|f=="VA"|f=="VT"] <- "S2"
      f[f=="WY"]                                         <- "S2"  
      f[f=="AK"|f=="UT"|f=="MS"]                         <- "S3"
      f[f=="NE"|f=="ND"]                                 <- "S4"
      f[f=="SD"|f=="SC"]                                 <- "S5"
    } 
    
    f <- as.factor(f)

    n.lev <- nlevels(f)
    print(paste("  nlevels = ",n.lev))

    f <- model.matrix(y ~ f - 1)   # Note: Intercept removed.
    f <- f[,2:ncol(f)]             # Note: First dummy deleted.

  } else {
    
    f.lrn<-read.table(fn.lrn,
             header=T,colClasses="numeric",blank.lines.skip=F)
    f.val<-read.table(fn.val,
             header=T,colClasses="numeric",blank.lines.skip=F)

    f <- c(f.lrn[,1],f.val[,1])

    f[is.na(f)] <- 0

    if (mod$feature[i]=="DOB") { 
      d <- f ; d[d>0] <- 1         # Note: Dummy for missing DOB
      f <- cbind(d,f,f^2)          # Note: Quadratic term added to DOB.
      rm(d)
    }
  }

  if (first.time) {
    X <- f
    first.time <- FALSE
  } else {
    X <- cbind(prev.X,f)
  }

  prev.X <- X

  rm(f.lrn,f.val,f)

  fit <- lm(y~X,weight=wts)        # Note: Validation sample not in fit.

  ehat.lrn <- fit$residuals[1:n.lrn]
  ehat.val <- fit$residuals[(n.lrn+1):n]
  ehat     <- fit$residuals

  sse.lrn <- sum(ehat.lrn^2) 
  sse.val <- sum(ehat.val^2) 
  sse     <- sum(ehat^2)     

  mse.lrn <- sse.lrn/(n.lrn - fit$rank) 
  mse.val <- sse.val/n.val  
  mse     <- mse.lrn*(fit$rank/n.lrn) + mse.val 

  rm(ehat.lrn,ehat.val,ehat)
  aov <- anova.lm(fit)

  print(paste("  rank    = ",fit$rank))
  print(paste("  P-value = ",aov[1,5]))
  print(paste("  mse.lrn = ",mse.lrn))
  print(paste("  mse.val = ",mse.val))
  print(paste("  mse     = ",mse))  
}

print(summary.lm(fit))

x0 <- c(0,n.lrn)
y0 <- c(0,sum(y.lrn-0.68))

x1 <- (1:n.lrn)
y1 <- fit$fitted.values[1:n.lrn]
y1 <- y1-0.68
y1 <- sort(y1)
y1 <- y1[n.lrn:1]
y1 <- cumsum(y1)

idx <- 1:n.lrn

print(paste("maximum profit in learning sample is ", max(y1)))
print(paste("maximum occurs at ", idx[y1==max(y1)]))

idx <- seq(1,n.lrn,length=200)

x1 <- x1[idx]
y1 <- y1[idx]

source("psopts.r");
postscript(file="cty_mod.eps");

plot(x=c(x0,x1),y=c(y0,y1),ylab="dollars",xlab="size of mailing",type="n")
lines(x=x0,y=y0,col="green")
lines(x=x1,y=y1,col="red")

dev.off()