## library(igraph)
library(rpart)
library(MASS)
library(cluster)
library(scatterplot3d)

######
## 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
####

## do the initialization
n = 50
tmp1 = rep(1,n)
tmpIdent = diag(n)
myX = cbind(tmp1,tmpIdent,tmpIdent)
myAlpha = 0
myF1 = rep(0,n)
myF2 = rep(0,n)

W1=exp(-as.matrix(daisy(as.matrix(x1)))^2/0.075)
W2=exp(-as.matrix(daisy(as.matrix(x2)))^2/0.075)
myDelta1=diag(apply(W1,1,sum))-W1
myDelta2=diag(apply(W2,1,sum))-W2
W=(W1+W2)/2

myBeta = cbind(myAlpha, t(myF1), t(myF2))



## start iterating
for(i in 1:50){
	myP = as.vector(exp(2*myX%*%t(myBeta))/(1+exp(2*myX%*%t(myBeta))))
	
	myAlpha = log(myP/(1-myP)) - myF1 - myF2
	myF1 = ginv(myDelta1)*(sum(myX*y) - sum(myX*myP)) 
	myF2 = ginv(myDelta2)*(sum(myX*y) - sum(myX*myP)) 
}

####
## Plot the true border
####
ltys = c(1,1)
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(myAlpha + myF1 + myF2, col="blue")


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

title("Binary Simulated Data")