# smooth decreasing density estimation using quadratic regression splines
# and least-squares criterion function
# user specified knots; support of density is range of knots so
# knots must span x-values
# mode is at smallest knot!

dspllsq=function(x,knots){
	n=length(x)
	t=knots-min(knots)
	x=x-min(knots)
	k=length(t)-2
	m=k+3
	delta=matrix(0,nrow=m,ncol=n)	
	delta[m,]=1:n*0+1
	xp=0:999/999*max(t)
	dp=matrix(0,nrow=m,ncol=1000)	

	for(j in 1:k){
	 	ind=x<=t[j]
	 	delta[j,ind] = 1
	 	ind=x>t[j]&x<=t[j+1]
	 	delta[j,ind] = 1 - (x[ind]-t[j])^2 / (t[j+2]-t[j]) / (t[j+1]-t[j])
	    ind=x>t[j+1]&x<=t[j+2]
		delta[j,ind] = (x[ind]-t[j+2])^2/(t[j+2]-t[j+1])/(t[j+2]-t[j])
	    ind=x>t[j+2]
	    delta[j,ind]=0
	}
	ind=x<=t[k+1]
	delta[k+1,ind]=1
	ind=x>t[k+1]
	delta[k+1,ind]=1-(x[ind]-t[k+1])^2/(t[k+2]-t[k+1])^2
	ind=x<t[2]
	delta[k+2,ind]=(t[2]-x[ind])^2/(t[2]-t[1])^2

# make grid for plotting smooth answers
	
	for(j in 1:k){
	 	ind=xp<=t[j]
	 	dp[j,ind] = 1
	 	ind=xp>t[j]&xp<=t[j+1]
	 	dp[j,ind] = 1 - (xp[ind]-t[j])^2 / (t[j+2]-t[j]) / (t[j+1]-t[j])
	    ind=xp>t[j+1]&xp<=t[j+2]
		dp[j,ind] = (xp[ind]-t[j+2])^2/(t[j+2]-t[j+1])/(t[j+2]-t[j])
	    ind=xp>t[j+2]
	    dp[j,ind]=0
	}
	ind=xp<=t[k+1]
	dp[k+1,ind]=1
	ind=xp>t[k+1]
	dp[k+1,ind]=1-(xp[ind]-t[k+1])^2/(t[k+2]-t[k+1])^2
	ind=xp<t[2]
	dp[k+2,ind]=(t[2]-xp[ind])^2/(t[2]-t[1])^2
	dp[m,]=1:1000*0+1
		
# create the H matrix

	hmat=matrix(nrow=m,ncol=m)
	for(j in 1:(m-3)){
		hmat[j,j] = t[j+1]-2*(t[j+1]-t[j])^2/(t[j+2]-t[j])/3 +(t[j+1]-t[j])**3/(t[j+2]-t[j])^2/5+(t[j+2]-t[j+1])^3/(t[j+2]-t[j])**2/5
	}
	hmat[m-2,m-2] = (8*t[m-1]+7*t[m-2])/15
	hmat[m-1,m-1]=(t[2]-t[1])/5
	hmat[m,m]=t[k+2]

	if(m>4){for(j in 1:(m-4)){
		hmat[j,j+1]=t[j+1] + ((t[j+2]-t[j+1])^2 - (t[j+1]-t[j])^2) / (t[j+2]-t[j])/3 -(t[j+2]-t[j+1])^3 / (t[j+3]-t[j+1])/(t[j+2]-t[j]) /30
	}}
	hmat[k,k+1]=t[k+1]+ ((t[k+2]-t[k+1])^2 - (t[k+1]-t[k])^2) / (t[k+2]-t[k])/3 -(t[k+2]-t[k+1])^2 /(t[k+2]-t[k]) /30
	
	
	if(m>4){for(j in 1:(m-4)){
		for(l in (j+2):(m-2)){
			hmat[j,l] = (t[j]+t[j+1]+t[j+2])/3
		}
	}}
	for(j in 1:(m-3)){
		hmat[j,m]=(t[j]+t[j+1]+t[j+2])/3
	}	
	for(j in 2:(m-2)){
		hmat[j,m-1]=(t[2]-t[1])/3
	}
	hmat[1,m-1]=t[2]/3-t[2]^2/30/t[3]
	hmat[m-2,m]=(2*t[m-1]+t[m-2])/3
	hmat[m-1,m]=(t[2]-t[1])/3

	for(i in 2:m){
		for(j in 1:(i-1)){
			hmat[i,j]=hmat[j,i]
		}
	}

# create the c-vector

	one=1:n*0+1/n
	cvec=delta%*%one
	
# cholesky decomposition and cone projection
	
	umat=chol(hmat)
	uinv=solve(umat)
	zvec=t(uinv)%*%cvec
	sol=coneproj(zvec,uinv)	
	bhat=uinv%*%sol$bhat
	theta=t(delta)%*%bhat	
	tspl=t(dp)%*%bhat
	
# compute lambda if area !=1

	if(!sol$h[m]){
		avec=hmat[m,]
		h=abs(bhat)>1e-10
		uj=umat[,h]
		pj=uj%*%solve(t(uj)%*%uj)%*%t(uj)
		lambda=(1-sum(avec*bhat))/(t(avec)%*%uinv%*%pj%*%t(uinv)%*%avec)
	
		zvec2=t(uinv)%*%(cvec+lambda*avec/2)
		sol2=coneproj(zvec2,uinv)	
		bhat2=uinv%*%sol2$bhat
		theta2=t(delta)%*%bhat2
		tspl2=t(dp)%*%bhat2
	}else{tspl2=tspl;theta2=theta}
	
	
	ans=new.env()
	ans$theta=theta2
	ans$tspl=tspl2
	ans$xp=xp+min(knots)
	ans$crit=t(bhat)%*%hmat%*%bhat-2*sum(cvec*bhat)
	ans
}


##############################################################
# cone projection code: hinge algorithm in R
#  find theta-hat to minimize || y-theta ||^2
#  subject to amat%*%theta >= 0
#
# amat must be irreducible
#
##############################################################
coneproj=function(y,amat){
	n=length(y);m=length(amat)/n
	sm=1e-8;h=1:m<0;obs=1:m;check=0
	delta=-amat;b2=delta%*%y
	if(max(b2)>sm){
		i=min(obs[b2==max(b2)])
		h[i]=TRUE
	}else{check=1;theta=1:n*0}
	while(check==0){
		xmat=matrix(delta[h,],ncol=n)
		a=solve(xmat%*%t(xmat))%*%xmat%*%y
		if(min(a)<(-sm)){
			avec=1:m*0;avec[h]=a
			i=min(obs[avec==min(avec)])
			h[i]=FALSE;check=0
		}else{
			check=1
			theta=t(xmat)%*%a
			b2=delta%*%(y-theta)
			if(max(b2)>sm){
				i=min(obs[b2==max(b2)])		
				h[i]=TRUE;check=0
			}
		}
	}
	bhat=y-theta
	ans=new.env()
	ans$bhat=bhat
	ans$h=!h
	ans$df=m-sum(h)
	ans
}