##----------------------------------------------------------
if(1) {
## load library
library(tree)
}
#--------------------------------------------------
if(0) {#boston data
library(MASS) ## a library of example datasets
attach(Boston)
n = dim(Boston)[1]
ddf = data.frame(x=lstat,y=medv)
oo=order(ddf$x)
ddf = ddf[oo,]
}
#--------------------------------------------------
if(1) { #simulated data
set.seed(99)
sigma=1
n=500
x=sort(-2 + 4*runif(n))
y = x^3+sigma*rnorm(n)
ddf=data.frame(x,y)
}
#--------------------------------------------------
if(1) {
par(mfrow=c(1,1))
plot(x,y)
readline("done?")
}
#--------------------------------------------------
if(1) {#do bootstrap loop
B=400
n=nrow(ddf)
nn = rep(0,B)
fmat=matrix(0,n,B)
set.seed(99)

tsleep1 = 10
tsleep2 = .5


par(mfrow=c(1,2))
for(i in 1:B) {
   ii = sample(1:n,n,replace=TRUE)
   nn[i] = length(unique(ii))
   bigtree = tree(y~x,ddf[ii,],mindev=.0002)
   #print(length(unique(bigtree$where)))
   temptree = prune.tree(bigtree,best=30)
   #print(length(unique(temptree$where)))
   fmat[,i]=predict(temptree,ddf)

   plot(ddf$x,ddf$y)
   iiu = unique(ii)
   nu = length(iiu)
   nvec=rep(0,nu)
   for(j in 1:nu) nvec[j] = sum(ii == iiu[j])
   points(ddf$x[iiu],ddf$y[iiu],col="red",pch=5,cex=nvec)
   #lines(ddf$x,fmat[,i],col=i,lwd=2)
   lines(ddf$x,fmat[,i],col='black',lwd=3)
   plot(temptree,type="uniform")
   if((i%%100)==0) {
      cat('i: ',i,'\n')
      print(table(nvec))
      cat("percent data in bootstrap sample:",length(iiu)/n)
      readline("go?")
   }
   if(i<5) {
      Sys.sleep(tsleep1)
   } else {
      Sys.sleep(tsleep2)
   }
   
}
readline("done?")
}

#--------------------------------------------------
### plot nn
par(mfrow=c(1,1))
plot(nn/n,xlab='bootstrap sample',ylab='percent obs in sample')
abline(h=mean(nn/n),col='blue',lty=3)
title(main=paste0('mean is ',round(mean(nn/n),2)))

#--------------------------------------------------
if(1) {#plot it
par(mfrow=c(1,1))
plot(ddf$x,ddf$y)
efit = apply(fmat,1,mean)
lines(ddf$x,efit,col='blue',lwd=4)
if(1) {
   lines(x,x^3,col="red")
}
for(i in 1:B) {
   lines(ddf$x,fmat[,i],col=i,lty=1)
   lines(ddf$x,efit,col='blue',lwd=4)
   Sys.sleep(.1)
}
}