#monotone regression spline code: hinge algorithm in R
# two parallel curves -- monotone increasing
#must have x= continuous predictor, dum=dummy predictor
#y=response,k=#interior knots then use the #following....  
# the x values must be sorted
# quadratic version only
#  uses golden ratio search

# Make the sigma vectors and first q-vectors

mspla=function(x,dum,y,k){
n=length(x)
sm=0.00000001
gold=(3-sqrt(5))/2
q1=(1:n*0+1)/sqrt(n)
q=q1

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

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


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

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

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

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

# start the search for the dummy variable parameter

a=min(y)-max(y)
c=max(y)-min(y)
b=a+(c-a)*gold

ysend=y-a*dum
yhata=mspl(x,ysend,sigma)
ssea=sum((ysend-yhata)^2)
ysend=y-b*dum
yhatb=mspl(x,ysend,sigma)
sseb=sum((ysend-yhatb)^2)
ysend=y-c*dum
yhatc=mspl(x,ysend,sigma)
ssec=sum((ysend-yhatc)^2)

while(abs(c-a)>.01){

if(c-b>b-a){d=c*gold+b*(1-gold)}
if(c-b<b-a){d=a*gold+b*(1-gold)}
ysend=y-d*dum
yhatd=mspl(x,ysend,sigma)
ssed=sum((ysend-yhatd)^2)

if(d<b){
	if(ssed<=sseb){
		c=b
		ssec=sseb
		b=d
		sseb=ssed
	}
	if(ssed>sseb){
		a=d
		ssea=ssed
	}
}
if(d>b){
	if(ssed<=sseb){
		a=b
		ssea=sseb
		b=d
		sseb=ssed
	}
	if(ssed>sseb){
		c=d
		ssec=ssed
	}
}
}
ysend=y-b*dum
yh=mspl(x,ysend,sigma)
ans=1:(n+1)*0
ans[1:n]=yh
ans[n+1]=b
ans
}


# Here is the hinge code:

mspl=function(x,y,sigma){
n=length(x)
m=length(sigma)/n
sm=0.00000001
l=1
h=c(1)
r=1
q=sigma[1,]/sqrt(sum(sigma[1,]^2))

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

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

while(check==0){

# Fit data to current EDGES

a=t(q)%*%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
q=q[,1:(i-1)]
for( j in i:l){
qnew=sigma[h[j],]-q%*%t(q)%*%sigma[h[j],]
qnew=qnew/sqrt(sum(qnew^2))
q=cbind(q,qnew)
}
r=t(q)%*%t(sigma[h,])
}
if(b[i]>(-sm)) {
c1=1
#
# now see if we need to add another hinge
#
theta=q%*%t(q)%*%y
rhat=y-theta
b2=sigma%*%rhat

# check to see if b2 negative

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