#  Code to test H_1:Ab>=0 versus H_2: unconstrained, using simulated
# null distribution of (Ab)_{min}
##  Model is y=Xb+e, X is full column rank n x k design matrix 
#  The matrix amat is the constraint matrix A; this must be
#    m by k and irreducible (Meyer 1999)
# nsim is the number of simulations used to estimate distribution
# returns constrained fit, unconstrained fit, p-value

conetest=function(y,xmat,amat,nsim){
	n=length(y)
	sm=1e-6
	m=length(xmat)/n
	xpx=t(xmat)%*%xmat
	umat=chol(xpx)
	uinv=solve(umat)
	atil=amat%*%uinv
	z=t(uinv)%*%t(xmat)%*%y
	ans=coneproj(z,atil)
	bhat=uinv%*%ans$phihat
	fhatc=xmat%*%bhat
	bmat=solve(xpx)%*%t(xmat)
	btil=bmat%*%y
	fhatuc=xmat%*%btil
	sighat=sqrt(sum((y-fhatc)^2)/(n-ans$dim))
	amin=min(amat%*%btil)
	if(amin<(-sm)){
		smin=1:nsim
		for(isim in 1:nsim){
			ysim=fhatc+rnorm(n)*sighat
			bsim=bmat%*%ysim
			smin[isim]=min(amat%*%bsim)
		}
		pval=sum(amin>smin)/nsim
	}
	if(amin>=0){pval=1}
	ans=new.env()
	ans$cfit=fhatc
	ans$ucfit=fhatuc
	ans$pval=pval
	ans
}	

# code for cone projection

#################################################
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
			}
		}
	}
	ans=new.env()
	ans$dim=n-sum(h)
	ans$phihat=y-theta
	ans
}