######
## f(x) Regression Example
######
set.seed(100)
x=runif(150,0,2*pi)
f<-function(x){sin(2*x)+log(exp(x)+2)}
y=f(x)+rnorm(150,,0.1)

plot(x,y,pch=16)
psi1=2.1
psi2=4.15
abline(v=c(psi1,psi2))

xgrid=seq(0,2*pi,length=500)
n=length(xgrid)
ygrid=vector(length=n)


####Compute (get h)

##Linear (intercpt)
linearInterceptGenerator<-function(x)
{
   as.matrix(cbind(as.numeric(x<psi1),as.numeric(x>=psi1 & x < psi2),as.numeric(x>=psi2)))
}

linearInterceptPredictor <-function(x0)
{
  cbind(as.numeric(x0<psi1), as.numeric(x0>=psi1 & x0<psi2), as.numeric(x0>=psi2)) 
}

## Linear (discontinuous)
linearDiscontinuousGenerator<-function(x)
{
  x=as.matrix(x)
  myColumn = x[,1]
  cbind(linearInterceptGenerator(x), as.numeric(x<psi1)*myColumn, as.numeric(x>=psi1 & x < psi2)*myColumn, as.numeric(x>=psi2)*myColumn)
}

linearDiscontinuousPredictor <-function(x0)
{
  cbind(linearInterceptPredictor(x0) ,as.numeric(x0<=psi1)*x0, as.numeric(x0>=psi1 & x0< psi2)*x0, as.numeric(x0>=psi2)*x0)
}
## Linear (continuous)
linearContinuousGenerator<-function(x)
{
  x=as.matrix(x)
  as.matrix(cbind(rep(1, dim(x)[1]), x[,1], (x[,1]-psi1) * as.numeric(x<=psi1), (x[,1]-psi2)*as.numeric(x<=psi2)))
}

linearContinuousPredictor <-function(x0)
{
  cbind(1, x0, (x0-psi1) * as.numeric(x0<=psi1), (x0-psi2) * as.numeric(x0<=psi2))
}

##Quadratic (continuous)
quadraticContinuousGenerator<-function(x)
{
  x=as.matrix(x)
  as.matrix(cbind(rep(1, dim(x)[1]), x[,1], x * x,  (x[,1] - psi1) * (x[,1] - psi1) * as.numeric(x<=psi1), (x[,1] - psi2) * (x[,1] - psi2) * as.numeric(x<=psi2)))
}

quadraticContinuousPredictor <-function(x0)
{
  cbind(1, x0, x0 * x0, (x0 - psi1) * (x0 - psi1) * as.numeric(x0 <= psi1), (x0 - psi2) * (x0 - psi2) * as.numeric(x0 <= psi2))
}

##Cubic (continuous)
cubicContinuousGenerator<-function(x)
{
  x=as.matrix(x)
  as.matrix(cbind(rep(1, dim(x)[1]), x[,1], x * x * x,  (x[,1] - psi1) * (x[,1] - psi1) * (x[,1] - psi1) * as.numeric(x<=psi1), (x[,1] - psi2) * (x[,1] - psi2) * (x[,1] - psi2) * as.numeric(x<=psi2)))
}

cubicContinuousPredictor <-function(x0)
{
  cbind(1, x0, x0 * x0 * x0, (x0 - psi1) * (x0 - psi1) * (x0 - psi1) * as.numeric(x0 <= psi1), (x0 - psi2) * (x0 - psi2) * (x0 - psi2) * as.numeric(x0 <= psi2))
}

##bs

##stop
h = cubicContinuousGenerator(x)
bls=solve(t(h)%*%h,t(h)%*%y)

for(i in 1:n){
  ## Set-up prediction (h0)
  x0=xgrid[i]

  ##Linear (intercept)
  ##h0 = linearInterceptPredictor(x0)

  ##Linear (discontinuous)
  ##h0 = linearDiscontinuousPredictor(x0)

  ##Linear (continuous)
  ##h0 = linearContinuousPredictor(x0)

  ##Quadratic (continuous)
  ##h0 = quadraticContinuousPredictor(x0)

  ##Cubic (continuous)
  h0 = cubicContinuousPredictor(x0)
  
  ##stop

  ##predict
  ygrid[i]=as.vector(h0%*%bls)### prediction
}

lines(xgrid,f(xgrid),col=c("blue"))
ind<-xgrid< psi1
lines(xgrid[ind],ygrid[ind],col=c("red"))

ind<-xgrid> psi1 & xgrid<psi2
lines(xgrid[ind],ygrid[ind],col=c("red"))

ind<-xgrid> psi2
lines(xgrid[ind],ygrid[ind],col=c("red"))