# convex regression code: hinge algorithm in R
#
creg=function(xs,ys){
	n=length(y)
	x=sort(xs)
	y=ys[order(xs)]
	m=n-2
	sm=1e-10;h=1:m<0;obs=1:m;check=0
	delta=matrix(0,nrow=n,ncol=m)
	x1=matrix(1,ncol=2,nrow=n)
	x1[,2]=x-mean(x)
	yhat1=x1%*%solve(t(x1)%*%x1)%*%t(x1)%*%y
	for(i in 3:n){
		for(j in i:n){delta[j,i-2]=x[j]-x[i-1]}
	}
	delta=t(delta-x1%*%solve(t(x1)%*%x1)%*%t(x1)%*%delta)
	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
			}
		}
	}
	theta+yhat1
}