fsr.fast.pv<-function(pv.orig,m,gam0=.05,digits=4){
# fast fsr from a sequence of forward addition p-values
# m is original number of variables
m1<-length(pv.orig)
vm<-1:m1
pvm<-rep(0,m1)                       # to create pvm below
for (i in 1:m1){pvm[i]<-max(pv.orig[1:i])}  # sequential max of pvalues
alpha<-c(0,pvm)
ng<-length(alpha)
S<-rep(0,ng)                         # will contain num. of true entering in orig.
for (ia in 2:ng){                    # loop through alpha values for S=size
S[ia] <- sum(pvm<=alpha[ia])         # size of models at alpha[ia], S[1]=0
}
ghat<-(m-S)*alpha/(1+S)              # gammahat_ER
# add additional points to make jumps
alpha2<-alpha[2:ng]-.0000001
ghat2<-(m-S[1:(ng-1)])*alpha2/(1+S[1:(ng-1)])
zp<-data.frame(a=c(alpha,alpha2),g=c(ghat,ghat2))
zp<-zp[order(zp$a),]
gmax<-max(zp$g)
index.max<-which.max(zp$g)           # index of largest ghat
alphamax<-zp$a[index.max]            # alpha with largest ghat
ind<-(ghat<= gam0 & alpha<=alphamax)*1
Sind<-S[max(which(ind > 0))]          # model size with ghat just below gam0
alphahat.fast<-(1+Sind)*gam0/(m-Sind)  # ER est.
size1<-sum(pvm<=alphahat.fast)+1       # size of model including intercept
plot(zp$a,zp$g,type="b",xlab="Alpha",ylab="Estimated Gamma",xlim=c(0,alphamax))
points(alphahat.fast,gam0,pch=19)
lines(c(-1,alphahat.fast),c(gam0,gam0))
lines(c(alphahat.fast,alphahat.fast),c(-1,gam0))
print(round(data.frame(pval=pv.orig,pvmax=pvm,ghigh=ghat2,glow=ghat[2:ng]),digits))
print(data.frame(m1,m,gam0,size=size1-1,alphamax,alphahat.fast))
return(list(zp=zp,alphahat.fast=alphahat.fast,alphamax=alphamax,gam0=gam0))
}