# 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)
#
#  nloop is the number of simulations to calculate the mixing distribution 

#  the function returns:
#         theta, the estimated f(x1,...x_L) at the x design points, and
#         alpha, the estimated coefficients for the Z columns
#  pvalue for: test of H_0: f_1=...=f_L=0
### versus H_1: full model, 
###  assuming normal errors

mtest=function(y,xmat,zmat,shapes,nloop){
	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,]
	conv=shapes>2&shapes<5
	vmat=cbind(zmat,xmat[,conv])
	pv=length(vmat)/n
	pmat=zmat%*%solve(t(zmat)%*%zmat)%*%t(zmat)
	proj=hingep(y,delta,vmat)
	alpha=proj$coefs[1:pv]
	bvec=proj$coefs[(pv+1):(pv+nd)]
	vhat=vmat%*%alpha
	th1=t(delta)%*%bvec+vhat
	th0=pmat%*%y
	sse0=sum((y-th0)^2)
	sse1=sum((y-th1)^2)
	bstat=(sse0-sse1)/sse0
#  mixing coefficients
	mdist=1:(m+1)*0
	for(i in 1:nloop){
		ysend=rnorm(n)
		prmix=hingep(ysend,delta,vmat)
		l=prmix$df-2
		mdist[l]=mdist[l]+1
	}
	mdist=mdist/nloop
	obs=1:(m+1)
	mz=max(obs[mdist>0])

	pval = mdist[1]
	for(d in 1:mz){
		alp=(d+pv-p)/2; bet=(n-d-pv)/2
		addl=pbeta(bstat,alp,bet)*mdist[d+1]
		pval = pval+addl
	}
	pval=1-pval
	ans=new.env()
	ans$pval=pval
	ans$yhat=yhat
	ans
}


######################################################
#  make cone edges
######################################################
makedelta=function(x,sh){
	n=length(x)
	xs=sort(x)
	xu=1:n*0
# find unique x values
	nu=1;xu[1]=xs[1];sm=1e-12
	for(i in 2:n){
		if(xs[i]>xs[i-1]+sm){
			nu=nu+1
			xu[nu]=xs[i]
		}
	}
	n1=nu
# 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[x==xu[i]]);c2=min(obs[x==xu[i+1]])
			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[nr]]);c2=min(obs[x==xu[nr-1]])
			amat[n1-1,c1]=-1;amat[n1-1,c2]=1
		}
		if(sh==7){
			amat=-amat
			c1=min(obs[x==xu[nr]]);c2=min(obs[x==xu[nr-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[x==xu[i],i]=1}
		atil=amat%*%wmat
		delta=t(wmat%*%t(atil)%*%solve(atil%*%t(atil)))
		dr=length(delta)/n
	}else{delta=t(amat)%*%solve(amat%*%t(amat))}
	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
	if(check==1){xmat=sigma[h,];a=solve(xmat%*%t(xmat))%*%xmat%*%y}
	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$df=sum(h)
ans$coefs=avec
ans$nrep=nrep
ans
}