# convex regression spline code: hinge algorithm in R
# must have x=predictor,y=response,k=#interior knots then use the 
# following....  the x values must be sorted and distinct
# quadratic version only

cspl=function(x,y,k){
n=length(x)
sm=0.00000001
q1=(1:n*0+1)/sqrt(n)
q2=(x-mean(x))/sqrt(sum((x-mean(x))^2))
qm=cbind(q1,q2)
l=2
h=c(1,2)
r=diag(c(1,1),nrow=2)

knots=round((0:(k+1))*n/(k+1))
knots[1]=1
t=x[knots]
m=k+3
sigma=matrix(1:(m*n),nrow=m,ncol=n)
sigma[1,]=1:n*0+1
sigma[2,]=t(q2)

for(j in 1:(m-3)){
 	for(i in 1:knots[j]){
 		sigma[j+2,i] = 0
 	}
 	for(i in (knots[j]+1):knots[j+1]){
 		sigma[j+2,i] = (x[i]-t[j])^2 / 2  
 	}	
    for(i in (knots[j+1]+1):n){
    	sigma[j+2,i]=(t[j+1]-t[j])^2/2 +(x[i]-t[j+1])*(t[j+1]-t[j])
    }
}
for(i in 1:knots[k+1]){
 sigma[m,i] = 0
}
for(i in (knots[k+1]+1):n){
 sigma[m,i] = (x[i]-t[m-2])^2 / 2  
}	

id=diag(1:n*0+1,nrow=n)
proj=id-qm%*%t(qm)
sigma=proj%*%t(sigma)
sigma=t(sigma)
sigma[1,]=t(q1)
sigma[2,]=t(q2)

# add one edge to start

check=0
rhat=y-qm%*%t(qm)%*%y
b2=sigma%*%rhat

if(b2[rank(b2)==m]>sm){
obs=1:m
i=obs[rank(b2)==m]
l=l+1
qnew=sigma[i,]-qm%*%t(qm)%*%sigma[i,]
qnew=qnew/sqrt(sum(qnew^2))
qm=cbind(qm,qnew)
h=c(h,i)
r=t(qm)%*%t(sigma[h,])
}
if(b2[rank(b2)==m]<sm){
check=1
}
# LOOP starts here:

while(check==0){

# Fit data to current EDGES

a=t(qm)%*%y
 
# check if convex: 
#first find the b vector:
b=1:l*0
b[l]=a[l]/r[l,l]
for( j in (l-1):1){
b[j]=a[j]
  for(i in (j+1):l){
    b[j]=b[j]-r[j,i]*b[i]
  }
b[j]=b[j]/r[j,j]
}

#check to see if b positive

obs=2:l
i=obs[rank(b[2:l])==1]
if(b[i]<(-sm)){
# if not, remove hinge, make new q and r
c1=0
h=c(h[1:(i-1)],h[(i+1):l])
l=l-1
qm=qm[,1:(i-1)]
for( j in i:l){
qnew=sigma[h[j],]-qm%*%t(qm)%*%sigma[h[j],]
qnew=qnew/sqrt(sum(qnew^2))
qm=cbind(qm,qnew)
}
r=t(qm)%*%t(sigma[h,])
}
if(b[i]>(-sm)) {
c1=1
#
# now see if we need to add another hinge
#
theta=qm%*%t(qm)%*%y
rhat=y-theta
b2=sigma%*%rhat

# check to see if b2 negative

obs=1:m
i=obs[rank(b2)==m]
if(b2[i]>sm){
l=l+1
qnew=sigma[i,]-qm%*%t(qm)%*%sigma[i,]
qnew=qnew/sqrt(sum(qnew^2))
qm=cbind(qm,qnew)
h=c(h,i)
r=t(qm)%*%t(sigma[h,])
c2=0
}
if(b2[i]<sm){c2=1}
check=c1*c2
h
}
}
yhat=qm%*%a
yhat
}