## library(igraph)
library(rpart)

######
## GAM example
######
set.seed(3)
n=50

x1<-runif(n,7*pi,13*pi)  ##7,13
x2<-runif(n,0,2*pi)
f<-function(x,y){
  f1<-sin(2*y)+pi/2
  f2<-exp(x/pi)
  f3<-1/(1+cos(y))
  eta=f1/f3+log(1+f2) 
  eta-mean(eta)
}
eta=f(x1,x2)+rnorm(n,0,1)
y=round( exp(eta)/(1+exp(eta)))


####
## Fit your GAM
####

## compute the starting values
rowLength = length(y)
alphaHat = log(mean(y)/(1-mean(y)))
myf1 = matrix(0, rowLength)
myf2 =matrix(0, rowLength)
trErr=te.er=matrix(0, rowLength)


## start iterating
for(i in 1:250){
	myEta1 = alphaHat + (myf1)
	myEta2 = alphaHat + (myf2)
	myP1 = 1/(1+exp(-myEta1))
	myP2 = 1/(1+exp(-myEta2))
	myZ1 = myEta1 + (y - myP1)/(myP1*(1-myP1))	
	myZ2 = myEta2 + (y - myP2)/(myP2*(1-myP2))	
	# now construct weights
	myW1 = myP1*(1-myP1)
	myW2 = myP2*(1-myP2)
	# fit an additive model to z's using w's
	h.x1=rpart(data=data.frame(y=myZ1,x1),control="rpart.control",minsplit=10,maxdepth=10,cp=-1, weights=myW1) # fit tree	
	h.x2=rpart(data=data.frame(y=myZ2,x2),control="rpart.control",minsplit=10,maxdepth=10,cp=-1, weights=myW2) # fit tree
# fit the tree and calculate the model
myf1=predict(h.x1,data.frame(x1))
myf2=predict(h.x2,data.frame(x2))
}

####
## Plot the true border
####
plot(x1,x2,col=y+2, pch=16, cex=2, xlab="x1", ylab="x2", ylim=c(0,8))##,xlim=c(min(x1),45))
x11<-seq(7*pi,13*pi,length=100)
x22<-seq(0,8,length=100)
z=outer(x11,x22,FUN=f)
contour(x11,x22,z,levels=0,add=TRUE,lty=ltys[1],method="edge",drawlabels=FALSE,col="purple",lwd=2)

####
## Use your gam to predict the x11, x22 grid
####
lines(myf1+myf2)


###
## Plot your prediction Border
####

title("Binary Simulated Data")