library(rpart)
set.seed(100542)

target.lrn <- read.table("../lrn/num/472.dat",header=T,colClasses="numeric")
target.val <- read.table("../val/num/472.dat",header=T,colClasses="numeric")
target.tst <- read.table("../tst/num/472.dat",header=T,colClasses="numeric")

 
y.lrn <- target.lrn[,1]
y.val <- target.val[,1]
y.tst <- target.tst[,1]
y     <- c(y.lrn,y.val,y.tst)

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

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

rm(target.lrn,target.val,target.tst)

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

idx.lrn <- 1:n.lrn
idx.val <- (n.lrn+1):(n.lrn+n.val)
idx.tst <- (n.lrn+n.val+1):n

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="")
  fn.tst <- paste("../tst/",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.tst <- read.table(fn.tst,
               header=T,colClasses="character",blank.lines.skip=F)

    f <- c(f.lrn[,1],f.val[,1],f.tst[,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=="VA"|f=="VT"|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)
    f.name <- levels(f) 
    if(n.lev==2) f.name <- c(mod$feature[i],mod$feature[i])
    print(paste("  nlevels = ",n.lev))

    f <- model.matrix(y ~ f - 1)      # Note: Intercept removed.
    f <- f[,2:n.lev]                  # Note: First dummy deleted.
    f.name <- f.name[2:n.lev]         # Note: First name 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.tst<-read.table(fn.tst,
             header=T,colClasses="numeric",blank.lines.skip=F)

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

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

    f.name <- mod$feature[i]

    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)
      f.name <- c("DOB.0","DOB.1","DOB.2")
    }
  }

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

  prev.X <- X
  prev.X.names <- X.names

}

rm(prev.X,prev.X.names)
rm(f.lrn,f.val,f)

dimnames(X) <- list(NULL,X.names)


# Regression fit to continuous target

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

yhat <- y.fit$fitted.values

ehat <- y - mean(y)
mse.nul <- mean(ehat[idx.lrn]^2) 

ehat <- y.fit$residuals
mse.lrn <- mean(ehat[idx.lrn]^2)
mse.val <- mean(ehat[idx.val]^2)
mse.tst <- mean(ehat[idx.tst]^2)

print("Regression fit to continuous target")
print(paste("  mse.nul = ",mse.nul))
print(paste("  mse.lrn = ",mse.lrn))
print(paste("  mse.val = ",mse.val))
print(paste("  mse.tst = ",mse.tst))

source("psopts.r");

for (alpha in c(00,10,20,30,40,50,60,80)) {

  for (i in 1:n.lrn) {
    if (b[i]>0) {
      wts[i] = 1
    } else {
      wts[i] = 1.0/(1.0 + alpha/100.0)
    }
  }
  
  # Regression fit to binary target
  
  b.fit <- lm(b~X,weights=wts)       # Note: Validation sample not in fit.
  
  bhat <- mat.or.vec(n,1) 
  bhat[b.fit$fitted.values>0.5] <- 1
  
  ehat <- (b != bhat)
  
  cer.nul <- mean(b[idx.lrn]) 
  cer.lrn <- mean(ehat[idx.lrn])
  cer.val <- mean(ehat[idx.val])
  cer.tst <- mean(ehat[idx.tst])
  
  print(paste("Regression fit to binary target, alpha =",alpha))
  print(paste("  cer.nul = ",cer.nul))
  print(paste("  cer.lrn = ",cer.lrn))
  print(paste("  cer.val = ",cer.val))
  print(paste("  cer.tst = ",cer.tst))
  
  
  # Tree fit to binary target
  
  c <- as.factor(b)
  
  ctrl <- rpart.control(minsplit=20,maxsurrogate=0,maxdepth=25,cp=0.0001,xval=0)
  tree <- rpart(c~X,weights=wts,method="class",control=ctrl)
  
  chat <- tree$frame[tree$where,5] - 1
  
  c <- b
  ehat <- (c != chat)
  
  cer.nul <- mean(b[idx.lrn])
  cer.lrn <- mean(ehat[idx.lrn])
  cer.val <- mean(ehat[idx.val])
  cer.tst <- mean(ehat[idx.tst])  
  
  print(paste("Tree fit to binary target, alpha =",alpha))
  print(paste("  cer.nul = ",cer.nul))
  print(paste("  cer.lrn = ",cer.lrn))
  print(paste("  cer.val = ",cer.val))
  print(paste("  cer.tst = ",cer.tst))
  
  
  x0 <- c(0,n.lrn)/n.lrn
  y0 <- c(0,sum(y.lrn-0.68))
  
  idx <- seq(1,n.lrn,length=500)
  
  x1 <- (1:n.lrn)/n.lrn
  x1 <- x1[idx]
  
  idx.yhat <- sort.list(yhat[idx.lrn],decreasing=TRUE)
  y1 <- y.lrn[idx.yhat]-0.68
  y1 <- cumsum(y1)
  y1 <- y1[idx]
  
  idx.bhat <- sort.list(bhat[idx.lrn],decreasing=TRUE)  
  y2 <- y.lrn[idx.bhat]-0.68
  y2 <- cumsum(y2)
  y2 <- y2[idx]
  
  idx.chat <- sort.list(chat[idx.lrn],decreasing=TRUE)  
  y3 <- y.lrn[idx.chat]-0.68
  y3 <- cumsum(y3)
  y3 <- y3[idx]
  
  filename <- paste("cty_class_lif_lrn_",alpha,".eps",sep="")
  postscript(file=filename);
  
  plot(x=c(x0,x1,x1),y=c(y0,y1,y3),ylab="dollars",xlab="proportion",type="n")
  lines(x=x0,y=y0,col="green")
  lines(x=x1,y=y1,col="blue")
  lines(x=x1,y=y2,col="red")
  lines(x=x1,y=y3,col="orange")
  
  dev.off()
  
  x0 <- c(0,n.val)/n.val
  y0 <- c(0,sum(y.val-0.68))
  
  idx <- seq(1,n.val,length=500)
  
  x1 <- (1:n.val)/n.val
  x1 <- x1[idx]
  
  idx.yhat <- sort.list(yhat[idx.val],decreasing=TRUE)
  y1 <- y.val[idx.yhat]-0.68
  y1 <- cumsum(y1)
  y1 <- y1[idx]
  
  idx.bhat <- sort.list(bhat[idx.val],decreasing=TRUE)  
  y2 <- y.val[idx.bhat]-0.68
  y2 <- cumsum(y2)
  y2 <- y2[idx]
  
  idx.chat <- sort.list(chat[idx.val],decreasing=TRUE)  
  y3 <- y.val[idx.chat]-0.68
  y3 <- cumsum(y3)
  y3 <- y3[idx]
  
  filename <- paste("cty_class_lif_val_",alpha,".eps",sep="")
  postscript(file=filename);

  
  plot(x=c(x0,x1,x1),y=c(y0,y1,y3),ylab="dollars",xlab="proportion",type="n")
  lines(x=x0,y=y0,col="green")
  lines(x=x1,y=y1,col="blue")
  lines(x=x1,y=y2,col="red")
  lines(x=x1,y=y3,col="orange")
  
  dev.off()
  
  filename <- paste("cty_class_tree_",alpha,".eps",sep="")
  postscript(file=filename);

  plot.rpart(tree,uniform=T,margin=0.01)
# text.rpart(tree)

  dev.off()

}  

for (alpha in c(00,10,20,30,40,50,60,80)) {

  for (i in 1:n.lrn) {
    if (b[i]>0) {
      wts[i] = 1
    } else {
      wts[i] = 1.0/(1.0 + alpha/100.0)
    }
  }
  
  # Regression fit to binary target
  
  b.fit <- lm(b~X,weights=wts)       # Note: Validation sample not in fit.
  
  bhat <- b.fit$fitted.values
  

  # Tree fit to binary target
  
  c <- as.factor(b)
  
  ctrl <- rpart.control(minsplit=20,maxsurrogate=0,maxdepth=25,cp=0.0001,xval=0)
  tree <- rpart(c~X,weights=wts,method="anova",control=ctrl)
  
  chat <- tree$frame[tree$where,5]
  
  x0 <- c(0,n.lrn)/n.lrn
  y0 <- c(0,sum(y.lrn-0.68))
  
  idx <- seq(1,n.lrn,length=500)
  
  x1 <- (1:n.lrn)/n.lrn
  x1 <- x1[idx]
  
  idx.yhat <- sort.list(yhat[idx.lrn],decreasing=TRUE)
  y1 <- y.lrn[idx.yhat]-0.68
  y1 <- cumsum(y1)
  y1 <- y1[idx]
  
  idx.bhat <- sort.list(bhat[idx.lrn],decreasing=TRUE)  
  y2 <- y.lrn[idx.bhat]-0.68
  y2 <- cumsum(y2)
  y2 <- y2[idx]
  
  idx.chat <- sort.list(chat[idx.lrn],decreasing=TRUE)  
  y3 <- y.lrn[idx.chat]-0.68
  y3 <- cumsum(y3)
  y3 <- y3[idx]
  
  filename <- paste("cty_prob_lif_lrn_",alpha,".eps",sep="")
  postscript(file=filename);
  
  plot(x=c(x0,x1,x1),y=c(y0,y1,y3),ylab="dollars",xlab="proportion",type="n")
  lines(x=x0,y=y0,col="green")
  lines(x=x1,y=y1,col="blue")
  lines(x=x1,y=y2,col="red")
  lines(x=x1,y=y3,col="orange")
  
  dev.off()
  
  x0 <- c(0,n.val)/n.val
  y0 <- c(0,sum(y.val-0.68))
  
  idx <- seq(1,n.val,length=500)
  
  x1 <- (1:n.val)/n.val
  x1 <- x1[idx]
  
  idx.yhat <- sort.list(yhat[idx.val],decreasing=TRUE)
  y1 <- y.val[idx.yhat]-0.68
  y1 <- cumsum(y1)
  y1 <- y1[idx]
  
  idx.bhat <- sort.list(bhat[idx.val],decreasing=TRUE)  
  y2 <- y.val[idx.bhat]-0.68
  y2 <- cumsum(y2)
  y2 <- y2[idx]
  
  idx.chat <- sort.list(chat[idx.val],decreasing=TRUE)  
  y3 <- y.val[idx.chat]-0.68
  y3 <- cumsum(y3)
  y3 <- y3[idx]
  
  filename <- paste("cty_prob_lif_val_",alpha,".eps",sep="")
  postscript(file=filename);
  
  plot(x=c(x0,x1,x1),y=c(y0,y1,y3),ylab="dollars",xlab="proportion",type="n")
  lines(x=x0,y=y0,col="green")
  lines(x=x1,y=y1,col="blue")
  lines(x=x1,y=y2,col="red")
  lines(x=x1,y=y3,col="orange")
  
  dev.off()
  
  filename <- paste("cty_prob_tree_",alpha,".eps",sep="")
  postscript(file=filename);

  plot.rpart(tree,uniform=T,margin=0.01)
# text.rpart(tree)

  dev.off()

}