#                                                                     
# R code to duplicate Figure~4.4 from the book ESLII
#
# Written by:                                                         
# --                                                                  
# John L. Weatherwax                2009-04-21                        
#                                                                     
# email: wax@alum.mit.edu                                             
#                                                                     
# Please send comments and especially bug reports to the              
# above email address.                                                
#                                                                     
#-----                                                                

source("load_vowel_data.R")
source("rda.R")

out    = load_vowel_data( TRUE, FALSE )
XTrain = out[[1]]
yTrain = out[[2]] 
XTest  = out[[3]]
yTest  = out[[4]] 

alphas = seq(0.0,1.0,length.out=100)
err_rate_train = c()
err_rate_test = c()
for( apha in alphas ){
  print( sprintf("Creating RDA classifier for alpha= %10.6f",apha) )  
  out = rda( XTrain, yTrain, XTest, yTest, apha )
  err_rate_train = c(err_rate_train, out[[2]])
  err_rate_test = c(err_rate_test, out[[4]])
}


#postscript("../../WriteUp/Graphics/Chapter4/dup_fig_4_7.eps", onefile=FALSE, horizontal=FALSE)

plot( alphas, err_rate_train, type="p", col="blue", ylim=range(c(err_rate_train,err_rate_test)), xlab="alpha", ylab="error rate", main="Duplication of Figure 4.7" )
lines( alphas, err_rate_test, type="p", col="red" )

#dev.off()

min_err_rate_spot = which.min( err_rate_test )
print( sprintf( "Min test error rate= %10.6f; alpha= %10.6f", err_rate_test[min_err_rate_spot], alphas[min_err_rate_spot] ) )


# run model selection with alpha and gamma models to combine to get Sigma_hat:
#
Nsamples = 100
alphas = seq(0.0,1.0,length.out=Nsamples)
gammas = seq(0.0,1.0,length.out=Nsamples)
err_rate_train = matrix( data=0, nrow=Nsamples, ncol=Nsamples )
err_rate_test = matrix( data=0, nrow=Nsamples, ncol=Nsamples )
for( ii in 1:Nsamples ){
  a = alphas[ii]
  for( jj in 1:Nsamples ){
    g = gammas[jj]
    print( sprintf("Creating RDA classifier for alpha= %10.6f, gamma= %10.6f",a,g) )  
    out = rda( XTrain, yTrain, XTest, yTest, a, gamma=g )
    err_rate_train[ii,jj] = out[[2]]
    err_rate_test[ii,jj] = out[[4]]
  }
}

inds = which( err_rate_test == min(err_rate_test), arr.ind=TRUE ); ii = inds[1]; jj = inds[2]
print( sprintf( "Min test error rate= %10.6f; alpha= %10.6f; gamma= %10.6f", err_rate_test[ii,jj], alphas[ii], gammas[jj] ) )