##Ordinary Rejection sampling

al<-2.75
be<-1
xv<-seq(0,10,by=0.01)
yv<-dgamma(xv,al,be)
plot(xv,yv,type="l",lwd=2,xlab="x")

x0<-2*(al-1)

M<-(2/gamma(al))*x0^(al-1)*exp(-x0/2)

y0<-M*dexp(xv,1/2)

plot(xv,yv,type="l",lwd=2,xlab="x",ylim=range(0,M/2),ylab="p(x)")
lines(xv,y0,col="red",lty=2,lwd=2)

1/M

avec<-c(1:8)
x0vec<-2*(avec-1)
Mvec<-(2/gamma(avec))*x0vec^(avec-1)*exp(-x0vec/2)
round(1/Mvec,4)

N<-100000
Y<-matrix(rexp(N*8,1/2),nrow=N,ncol=8)
Gmat<-Emat<-Umat<-Amat<-Y*0
for(j in 1:8){
	Gmat[,j]<-dgamma(Y[,j],avec[j])
	Emat[,j]<-dexp(Y[,j],1/2)
	Umat[,j]<-runif(N)
	Amat[,j]<-as.numeric(Umat[,j]<=Gmat[,j]/(Emat[,j]*Mvec[j]))
}
apply(Amat,2,mean)



plot(xv,yv,type="l",lwd=2,xlab="x",ylim=range(0,M/2),ylab="p(x)")
lines(xv,y0,col="red",lty=2,lwd=2)
polygon(c(xv,rev(xv)),c(yv,0*yv),col="gray")
polygon(c(xv,rev(xv)),c(yv,rev(y0)),col="cyan")
legend(6,0.8,c("Accepted","Rejected"),col=c("gray","cyan"),pch=15,cex=1.5)

set.seed(1010)
X1<-rexp(1,1/2)
pX1<-dgamma(X1,al,be)
pX1.0<-dexp(X1,al)

pX1/(M*pX1.0)
U1<-runif(1)
plot(xv,yv,type="l",lwd=2,xlab="x",ylim=range(0,M/2))
lines(xv,y0,col="red",lty=2,lwd=2)
points(X1,(U1*M*pX1.0),pch=19)

plot(xv,log(yv),type="l",lwd=2,xlab="x",ylim=range(-10,1))
lines(xv,log(y0),col="red",lty=2,lwd=2)
points(X1,log(U1*M*pX1.0),pch=19)

###########################################################
#Adaptive Rejection

xv<-seq(0,12,by=0.01)
yv<-dgamma(xv,al,be)
xp<-c(0,0.25,0.75,2,5,10)
yp<-dgamma(xp,al,be)

#Log density
plot(xv,log(yv),type="l",lwd=2,xlab="x",ylim=range(-10,1),ylab=expression(log(p(x))))

#Initial Lower Squeeze
plot(xv,log(yv),type="l",lwd=2,xlab="x",ylim=range(-10,1),ylab=expression(log(p(x))))
lines(xp,log(yp),col="red",lwd=2)

#Upper and lower tangents
xpts<-rep(" ",length(xp))
xpts[1]<-expression(x[0])
xpts[2]<-expression(x[1])
xpts[3]<-expression(x[2])
xpts[4]<-expression(x[3])
xpts[5]<-expression(x[4])
plot(xv,log(yv),type="l",lwd=2,xlab="x",ylim=range(-10,1),ylab=expression(log(p(x))))
abline(v=xp,lty=3,col="red")
for(ip in 2:5){
	x1<-xp[ip]
	x2<-xp[ip+1]
	y1<-log(yp[ip])
	y2<-log(yp[ip+1])
	m<-(y2-y1)/(x2-x1)
	cval<-y1-m*x1
	abline(a=cval,b=m,lty=2,col="red",lwd=2)
}
text(xp+0.05,xp*0-10.1,xpts)


#################################################################
xv<-seq(0,12,by=0.01)
yvl<-log(dgamma(xv,al,be))

set.seed(335)
(xp<-sort(c(0,runif(6,0,6),Inf)))

ypl<-log(dgamma(xp,al,be))
xpts<-rep(" ",length(xp))
xpts[1]<-expression(x[0])
xpts[2]<-expression(x[1])
xpts[3]<-expression(x[2])
xpts[4]<-expression(x[3])
xpts[5]<-expression(x[4])
xpts[6]<-expression(x[5])
xpts[7]<-expression(x[6])
xpts[8]<-expression(x[7])
uv<-uv.l<-uv.r<-rep(NA,length(xv))
lv<-rep(-Inf,length(xv))

nint<-length(xp)-1

max.val<-max(yvl)

slope.vec<-rep(0,nint)
intercept.vec<-rep(0,nint)

for(i in 2:length(xv)){

	xi<-xv[i]
	ip<-as.numeric(cut(xi,xp))
	if(ip == 1){
		#No lower squeeze
		x1l<-xp[ip+1]
		x2l<-xp[ip+2]
		y1l<-ypl[ip+1]
		y2l<-ypl[ip+2]
		m<-(y2l-y1l)/(x2l-x1l)		
		uv[i]<-y1l+m*(xi-x1l)

		slope.vec[ip]<-m
		intercept.vec[ip]<-y1l-m*x1l

	}else if(ip == nint){
		#No lower squeeze
		x1r<-xp[ip-1]
		x2r<-xp[ip]
		y1r<-ypl[ip-1]
		y2r<-ypl[ip]
		m<-(y2r-y1r)/(x2r-x1r)
		uv[i]<-y1r+m*(xi-x1r)

		slope.vec[ip]<-m
		intercept.vec[ip]<-y1r-m*x1r


	}else if(ip == 2){
		#Lower squeeze
		#Upper squeeze from right-hand interval
	
		x1<-xp[ip]
		x2<-xp[ip+1]
		y1<-ypl[ip]
		y2<-ypl[ip+1]
		m<-(y2-y1)/(x2-x1)
		lv[i]<-y1+m*(xi-x1)

		x1l<-xp[ip+1]
		x2l<-xp[ip+2]
		y1l<-ypl[ip+1]
		y2l<-ypl[ip+2]
		m<-(y2l-y1l)/(x2l-x1l)		
		uv[i]<-y1l+m*(xi-x1l)

		slope.vec[ip]<-m
		intercept.vec[ip]<-y1l-m*x1l



	}else if(ip == nint-1){
		#Lower squeeze
		#Upper squeeze from left-hand interval

		x1<-xp[ip]
		x2<-xp[ip+1]
		y1<-ypl[ip]
		y2<-ypl[ip+1]
		m<-(y2-y1)/(x2-x1)
		lv[i]<-y1+m*(xi-x1)

		x1l<-xp[ip-1]
		x2l<-xp[ip]
		y1l<-ypl[ip-1]
		y2l<-ypl[ip]
		m<-(y2l-y1l)/(x2l-x1l)
		uv[i]<-y1l+m*(xi-x1l)

		slope.vec[ip]<-m
		intercept.vec[ip]<-y1l-m*x1l



	}else{

		#Lower squeeze
		#Upper squeeze from right-hand interval when left of mode
		#Upper squeeze from left-hand interval when right of mode

		x1<-xp[ip]
		x2<-xp[ip+1]
		y1<-ypl[ip]
		y2<-ypl[ip+1]
		m<-(y2-y1)/(x2-x1)
		lv[i]<-y1+m*(xi-x1)

		x1l<-xp[ip-1]
		x2l<-xp[ip]
		y1l<-ypl[ip-1]
		y2l<-ypl[ip]
		m<-(y2l-y1l)/(x2l-x1l)
		uv.l[i]<-y1l+m*(xi-x1l)
		x1r<-xp[ip+1]
		x2r<-xp[ip+2]
		y1r<-ypl[ip+1]
		y2r<-ypl[ip+2]
		m<-(y2r-y1r)/(x2r-x1r)
		uv.r[i]<-y1r+m*(xi-x1r)

		if(y2 <= max.val){
			uv[i]<-uv.r[i]
			slope.vec[ip]<-m
			intercept.vec[ip]<-y1r-m*x1r
		}else{
			uv[i]<-uv.l[i]
			slope.vec[ip]<-m
			ntercept.vec[ip]<-y1l-m*x1l
		}
		

	
	} 

}

plot(xv,yvl,type="l",lwd=2,xlab="x",ylim=range(-10,1),ylab=expression(log(p(x))))
abline(v=xp[1:8],lty=3,col="red")
text(xp[1:8]+0.05,xp[1:8]*0-10.1,xpts)
lines(xv,lv,col="blue")
lines(xv,uv.l,col="red",lwd=2)
lines(xv,uv.r,col="green",lwd=2)
lines(xv,uv,col="orange",lwd=2)

plot(xv,exp(yvl),type="l",lwd=2,xlab="x",ylim=range(0,1),ylab=expression(p(x)))
abline(v=xp[1:8],lty=3,col="red")
text(xp[1:8]+0.05,xp[1:7]*0,xpts)
lines(xv,exp(lv),col="blue")
lines(xv,exp(uv),col="orange",lwd=2)


piecewise.func<-function(xv,xpo,mv,cv){
	xcat<-as.numeric(cut(xv,xpo,include.lowest=TRUE))
	yv<-mv[xcat]*xv+cv[xcat]
	return(yv)
}

xvals<-seq(0,12,by=0.001)
yvals<-piecewise.func(xvals,xp,slope.vec,intercept.vec)
plot(xv,yvl,type="l",lwd=2,xlab="x",ylim=range(-10,1),ylab=expression(log(p(x))))
abline(v=xp[1:8],lty=3,col="red")
text(xp[1:8]+0.05,xp[1:7]*0-10.1,xpts)
lines(xvals,yvals,col='orange')

plot(xv,exp(yvl),type="l",lwd=2,xlab="x",ylim=range(0,1),ylab=expression(p(x)))
abline(v=xp[1:8],lty=3,col="red")
text(xp[1:8]+0.05,xp[1:7]*0-10.1,xpts)
lines(xvals,exp(yvals),col='orange')



find.parameters<-function(xvv,ypv,xpv){

	nint<-length(xpv)-1

	uvv<-uvv.l<-uvv.r<-rep(NA,length(xvv))
	lvv<-rep(-Inf,length(xvv))
	svec<-rep(0,nint)
	ivec<-rep(0,nint)

	for(i in 2:length(xvv)){

		xi<-xvv[i]
		ip<-as.numeric(cut(xi,xpv))
		if(ip == 1){
			#No lower squeeze
			x1l<-xpv[ip+1]
			x2l<-xpv[ip+2]
			y1l<-ypv[ip+1]
			y2l<-ypv[ip+2]
			m<-(y2l-y1l)/(x2l-x1l)		
			uvv[i]<-y1l+m*(xi-x1l)
			svec[ip]<-m
			ivec[ip]<-y1l-m*x1l
		}else if(ip == nint){
			#No lower squeeze
			x1r<-xpv[ip-1]
			x2r<-xpv[ip]
			y1r<-ypv[ip-1]
			y2r<-ypv[ip]
			m<-(y2r-y1r)/(x2r-x1r)
			uvv[i]<-y1r+m*(xi-x1r)
			svec[ip]<-m
			ivec[ip]<-y1r-m*x1r
		}else if(ip == 2){
			#Lower squeeze
			#Upper squeeze from right-hand interval
			x1<-xpv[ip]
			x2<-xpv[ip+1]
			y1<-ypv[ip]
			y2<-ypv[ip+1]
			m<-(y2-y1)/(x2-x1)
			lvv[i]<-y1+m*(xi-x1)
			x1l<-xp[ip+1]
			x2l<-xp[ip+2]
			y1l<-ypl[ip+1]
			y2l<-ypl[ip+2]
			m<-(y2l-y1l)/(x2l-x1l)		
			uvv[i]<-y1l+m*(xi-x1l)
			svec[ip]<-m
			ivec[ip]<-y1l-m*x1l
		}else if(ip == nint-1){
			#Lower squeeze
			#Upper squeeze from left-hand interval

			x1<-xpv[ip]
			x2<-xpv[ip+1]
			y1<-ypv[ip]
			y2<-ypv[ip+1]
			m<-(y2-y1)/(x2-x1)
			lvv[i]<-y1+m*(xi-x1)
			x1l<-xpv[ip-1]
			x2l<-xpv[ip]
			y1l<-ypv[ip-1]
			y2l<-ypv[ip]
			m<-(y2l-y1l)/(x2l-x1l)
			uvv[i]<-y1l+m*(xi-x1l)
			svec[ip]<-m
			ivec[ip]<-y1l-m*x1l
		}else{

			#Lower squeeze
			#Upper squeeze from right-hand interval when left of mode
			#Upper squeeze from left-hand interval when right of mode

			x1<-xpv[ip]
			x2<-xpv[ip+1]
			y1<-ypv[ip]
			y2<-ypv[ip+1]
			m<-(y2-y1)/(x2-x1)
			lvv[i]<-y1+m*(xi-x1)
			x1l<-xpv[ip-1]
			x2l<-xpv[ip]
			y1l<-ypv[ip-1]
			y2l<-ypv[ip]
			m<-(y2l-y1l)/(x2l-x1l)
			uvv.l[i]<-y1l+m*(xi-x1l)
			x1r<-xpv[ip+1]
			x2r<-xpv[ip+2]
			y1r<-ypv[ip+1]
			y2r<-ypv[ip+2]
			m<-(y2r-y1r)/(x2r-x1r)
			uvv.r[i]<-y1r+m*(xi-x1r)
			if(y2 <= max.val){
				uvv[i]<-uvv.r[i]
				svec[ip]<-m
				ivec[ip]<-y1r-m*x1r
			}else{
				uvv[i]<-uvv.l[i]
				svec[ip]<-m
				ivec[ip]<-y1l-m*x1l
			}
		} 
	}

	return(list(slope=svec,intercept=ivec,lower=lvv,upper=uvv))

}

xp.new<-c(0,1,2,3,4,5,Inf)
ypl.new<-log(dgamma(xp.new,al,be))

current<-find.parameters(xvv=xv,ypv=ypl.new,xpv=xp.new)

plot(xv,exp(yvl),type="l",lwd=2,xlab="x",ylim=range(0,1),ylab=expression(p(x)))
abline(v=xp.new,lty=3,col="red")
lines(xv,exp(current$upper),col='orange',lwd=2)
lines(xv,exp(current$lower),col='blue',lwd=2)

integrate.piecewise<-function(svec,ivec,xpv){

	cdfv<-xpv*0
	for(j in 2:length(xpv)){
		cdfv[j]<-cdfv[j-1]+exp(ivec[j-1])*(exp(xpv[j]*svec[j-1])-exp(xpv[j-1]*svec[j-1]))/svec[j-1]
	}
	return(list(normalizing.const=cdfv[length(xpv)]))
}

cval<-integrate.piecewise(slope.new,int.new,xp.new)$normalizing.const

plot(xv,exp(yvl),type="l",lwd=2,xlab="x",ylim=range(0,1),ylab=expression(p(x)))
abline(v=xp.new,lty=3,col="red")
lines(xv,exp(current$upper)/cval,col='orange',lwd=2)
lines(xv,exp(current$lower),col='blue',lwd=2)


sample.piecewise<-function(nsamp,svec,ivec,xpv){

	cdfv<-xpv*0

	for(j in 2:length(xpv)){
		cdfv[j]<-cdfv[j-1]+exp(ivec[j-1])*(exp(xpv[j]*svec[j-1])-exp(xpv[j-1]*svec[j-1]))/svec[j-1]
	}
	cdfv<-cdfv/cdfv[length(xpv)]


	

	return(cdf=cdfv)
}


slope.new<-current$slope
int.new<-current$intercept
sample.piecewise(1,slope.new,int.new,xp.new)