# Additive shape-restricted regression
#  y=f_1(x1)+f-2(x2)+...+f_L(x_L)+Za+e
#  xmat is n x L full column rank -- continuous, constrained predictors
#  zmat must be full column rank and column space contains constant vectors
#  [xmat | zmat] must be full column rank
#  if there are no covariates, make zmat=matrix(1,ncol=1,nrow=n)
#  the matrix xinterp is an M x L matrix of points at which f(x1,...x_L) is predicted
#   where f(x1,...x_L)=f_1(x1)+f-2(x2)+...+f_L(x_L).

#  the shapes are specified in a vector of L integers, where
#  
#    1=increasing
#    2=decreasing
#    3=convex
#    4=concave
#    5=increasing convex
#    6=decreasing convex
#    7=increasing concave
#    8=decreasing concave

#  for example if L=3 and E(y) is increasing in all three predictors, shapes=(1,1,1)

#  skip=TRUE means that the code does not check the irreducibility of the edges.
#  faster but might not terminate correctly

#  the function returns:
#         theta, the estimated f(x1,...x_L) at the x design points, and
#         alpha, the estimated coefficients for the Z columns
#         thint, the estimated f(x1,...,x_L) at the rows of xinterp
#
addshape=function(y,xmat,zmat,shapes,xinterp,skip){
	n=length(y)
	p=length(zmat)/n
	capl=length(xmat)/n
	sm=1e-8
#  create delta matrix
	delta=matrix(nrow=(n-1)*capl,ncol=n)
	belong=1:((n-1)*capl)*0
	nd=0
	for(l in 1:capl){
		del=makedelta(xmat[,l],shapes[l])
		nl=length(del)/n
		delta[(nd+1):(nd+nl),]=del
		belong[(nd+1):(nd+nl)]=l
		nd=nd+nl
	}
	delta=delta[1:nd,]
	belong=belong[1:nd]
	conv=shapes>2&shapes<5
	vmat=cbind(zmat,xmat[,conv])
	pv=length(vmat)/n
#  get irreducible set of edges
	if(!skip){
	keep=1:nd>0
	for(i in 1:nd){
		cols=1:nd!=i&keep
		others=delta[cols,]
		ans=hingep(delta[i,],others,vmat)
		if(max(abs(delta[i,]-ans$yhat))<sm){keep[i]=FALSE}
	}
	print("finished irred")
	delta=delta[keep,]
	belong=belong[keep]
	nd=sum(keep)
	}
	proj=hingep(y,delta,vmat)
	alpha=proj$coefs[1:pv]
	bvec=proj$coefs[(pv+1):(pv+nd)]
	vhat=vmat%*%alpha
	theta=t(delta)%*%bvec+vhat
#  interpolate!
	thetahat=matrix(nrow=n,ncol=capl)
	for(l in 1:capl){
		thetahat[,l]=t(delta[belong==l,])%*%bvec[belong==l]
		if(conv[l]){thetahat[,l]=thetahat[,l]+xmat[,l]*alpha[p+sum(conv[1:l])]}
	}
	thetahat[,1]=thetahat[,1]+alpha[1]
	capm=length(xinterp)/capl
	thint=1:capm*0
	obs=1:n
	for(m in 1:capm){
		xo=xinterp[m,]
		for(l in 1:capl){
			if(xo[l]<max(xmat[,l])){
				bigger=min(xmat[xmat[,l]>xo[l],l])
			}else{bigger=max(xmat[,l])}
			if(xo[l]>min(xmat[,l])){
				smaller=max(xmat[xmat[,l]<=xo[l],l])
			}else{smaller=min(xmat[,l])}
			i1=min(obs[xmat[,l]==bigger])
			i2=min(obs[xmat[,l]==smaller])
			al=(xo[l]-xmat[i1,l])/(xmat[i2,l]-xmat[i1,l])
			thint[m]=thint[m]+(1-al)*thetahat[i1,l]+al*thetahat[i2,l]
		}
	}
	ans=new.env()
	ans$theta=theta
	ans$alpha=alpha
	ans$thint=thint
	ans
}

######################################################
#  make cone edges
######################################################
makedelta=function(x,sh){
	n=length(x)
	xs=sort(x)
	xu=1:n*0
# find unique x values
	xs=round(xs*1e12)/1e12
	xu=unique(xs)
	n1=length(xu)
	sm=1e-8
# make bmat --equality constraints
	obs=1:n
	if(n1<n){
		bmat=matrix(0,nrow=n-n1,ncol=n)
		row=0
		for(i in 1:n1){
			cobs=obs[x==xu[i]]
			nr=length(cobs)
			if(nr>1){
				for(j in 2:nr){
					row=row+1
					bmat[row,cobs[1]]=-1;bmat[row,cobs[j]]=1
				}
			}
		}	
	}
#  increasing or decreasing
	if(sh<3){
		amat=matrix(0,nrow=n1-1,ncol=n)
		for(i in 1:(n1-1)){
			c1=min(obs[abs(x-xu[i])<sm]);c2=min(obs[abs(x-xu[i+1])<sm])
			amat[i,c1]=-1;amat[i,c2]=1
		}
		if(sh==2){amat=-amat}
	}else if(sh==3|sh==4){
#  convex or concave
		amat=matrix(0,nrow=n1-2,ncol=n)
		for(i in 1:(n1-2)){
			c1=min(obs[x==xu[i]]);c2=min(obs[x==xu[i+1]]);c3=min(obs[x==xu[i+2]])
			amat[i,c1]=xu[i+2]-xu[i+1];amat[i,c2]=xu[i]-xu[i+2];amat[i,c3]=xu[i+1]-xu[i]
		}
		if(sh==4){amat=-amat}
	}else if(sh>4){
		amat=matrix(0,nrow=n1-1,ncol=n)
		for(i in 1:(n1-2)){
			c1=min(obs[x==xu[i]]);c2=min(obs[x==xu[i+1]]);c3=min(obs[x==xu[i+2]])
			amat[i,c1]=xu[i+2]-xu[i+1];amat[i,c2]=xu[i]-xu[i+2];amat[i,c3]=xu[i+1]-xu[i]
		}
		if(sh==5){
			c1=min(obs[x==xu[1]]);c2=min(obs[x==xu[2]])
			amat[n1-1,c1]=-1;amat[n1-1,c2]=1
		}
		if(sh==6){
			c1=min(obs[x==xu[n1]]);c2=min(obs[x==xu[n1-1]])
			amat[n1-1,c1]=-1;amat[n1-1,c2]=1
		}
		if(sh==7){
			amat=-amat
			c1=min(obs[x==xu[n1]]);c2=min(obs[x==xu[n1-1]])
			amat[n1-1,c1]=1;amat[n1-1,c2]=-1
		}
		if(sh==8){
			amat=-amat
			c1=min(obs[x==xu[1]]);c2=min(obs[x==xu[2]])
			amat[n1-1,c1]=1;amat[n1-1,c2]=-1
		}
	}
	if(n1<n){
		wmat=matrix(0,nrow=n,ncol=n1)
		for(i in 1:n1){wmat[abs(x-xu[i])<sm,i]=1}
		atil=amat%*%wmat
		delta=t(wmat%*%t(atil)%*%solve(atil%*%t(atil)))
	}else{delta=solve(amat%*%t(amat))%*%amat}
	dr=length(delta)/n
	if(sh>2&sh<5){
		pr1=cbind(1:n*0+1,x)
		prmat=pr1%*%solve(t(pr1)%*%pr1)%*%t(pr1)
		for(i in 1:dr){delta[i,]=delta[i,]-t(prmat%*%delta[i,])}
	}else{
		for(i in 1:dr){delta[i,]=delta[i,]-mean(delta[i,])}
	}
	for(i in 1:dr){delta[i,]=delta[i,]/sqrt(sum(delta[i,]^2))}
	delta
}

#######################################################
# cone projection code: hinge algorithm in R
#  rows of delta are generators of the cone
# the linear space contained in the constraint set
# is spanned by columns of vmat (full col rank)
######################################################
hingep=function(y,delta,vmat){
	n=length(y)
	m=length(delta)/n
	p=length(vmat)/n
	sigma=matrix(nrow=m+p,ncol=n)
	sigma[1:p,]=t(vmat)
	sigma[(p+1):(m+p),]=delta
	sm=1e-8
	h=1:(m+p)<0
	h[1:p]=TRUE
	obs=1:(m+p)
	check=0
	theta=vmat%*%solve(t(vmat)%*%vmat)%*%t(vmat)%*%y
	b2=sigma%*%(y-theta)
	if(max(b2)>sm){
		i=min(obs[b2==max(b2)])
		h[i]=TRUE
	}
	if(max(b2)<sm){
		check=1;theta=1:n*0
	}
	nrep=0
	while(check==0){
		nrep=nrep+1
		xmat=sigma[h,]
		if(length(xmat)==n){xmat=matrix(sigma[h,],ncol=n)}
		a=solve(xmat%*%t(xmat))%*%xmat%*%y
		if(min(a[(p+1):length(a)])<(-sm)){
			avec=1:(m+p)*0
	 		avec[h]=a
			i=max(obs[avec==min(avec[(p+1):(p+m)])])
			h[i]=FALSE
			check=0
		}
		if(min(a[(p+1):length(a)])>(-sm)) {
			check=1
			theta=t(xmat)%*%a
			b2=sigma%*%(y-theta)
#			obs=1:m
			if(max(b2)>sm){
				i=min(obs[b2==max(b2)])
				check=0
				h[i]=TRUE
			}
		}
	}
ans=new.env()
ans$yhat=theta
avec=1:(m+p)*0
avec[h]=a
ans$coefs=avec
ans$nrep=nrep
ans
}