#######################################################################
# R function: This function fits a logisitic model to group testing 
#             data according to the Bayesian methods discussed in McMahan et al. 
#
# Input:
# Z = A matrix of testing responses whose jth row is of the form (col1=Z_j, col2=cj, col3=assay used, col4:col(4+cj-1)=indices of the individuals assigned to the jth pool) 
# X = Covariate matrix in which the ith row contains the covariate information pertaining to the ith individual
# Y = matrix whose ith row is of the form (col1=0, col2=number of pools ith individual involved in, col3:col(3+col2-1)=pools ith individual was assigned to)
# GI = The number of Gibbs iterates 
# b0 = Initialization value for the regression parameters
# na = Number of different assays being used
# ase(bse) = prior parameters for the testing sensitivities, if these values are to be estimated, note these should be of length na (Assuming a beta prior)
# asp(bsp) = prior parameters for the testing specificities, if these values are to be estimated, note these should be of length na (Assuming a beta prior)
# Se = vector of sensitivity values, if known
# Sp = vector of specificity values, if known
# est.error = logical if FALSE function uses the known Se and Sp to fit model, if TRUE Se and Sp are estimated along with the other parameters
# (a,R)     = The prior mean vector (a) and covariance matrix (R)
#
#
# Required Packages: mnormt
#
#
#

Bayes.reg.GT.ip<-function(Z,Y,X,GI,b0,na, ase=rep(1,na), bse=rep(1,na), asp=rep(1,na), bsp=rep(1,na), Se=NULL, Sp=NULL, est.error=FALSE, visual=FALSE, step=50, a, R){


J<-dim(Z)[1]
N<-dim(X)[1]

beta<-matrix(-99,nrow=GI,ncol=length(b0))

if(est.error==TRUE){
Se.mat<-matrix(-99,nrow=GI,ncol=na)
Sp.mat<-matrix(-99,nrow=GI,ncol=na)
}


###############################
# Initializing latent variables
Y[,1]<-rbinom(N,1,exp(X%*%b0)/(1+exp(X%*%b0)))

#########################
# Load Fortran DLL
dyn.load("samplatentstatus.dll")
dyn.load("errorupdate.dll")

##########################
#GIBBS SAMPLER STARTS HERE

for(g in 1:GI){


######################################################################################################## 
#Sampling the beta parameters (This assumes a diffuse prior easily changeable see Gamerman (1997)
beta[g,]<-WLS.Metropolis.Hastings.ip(b0=b0,X=X,y=Y[,1],a,R)
b0<-beta[g,]


if(est.error==TRUE){
###################
#Sampling Se and Sp

PSe<-matrix(0,nrow=na,ncol=2)
PSp<-matrix(0,nrow=na,ncol=2)
zm<-dim(Z)[2]
ym<-dim(Y)[2]

res<-.C("errorupdate_",
            as.integer(N),
            as.integer(J),
            as.integer(zm),
            as.integer(ym),
            as.integer(Y),
            as.integer(Z),
            as.integer(PSe),
            as.integer(PSe),
            as.integer(na))
Se.up<-matrix(res[[7]],na,2)
Sp.up<-matrix(res[[8]],na,2)

for(t in 1:na){
Se.mat[g,t]<-rbeta(1,ase[t]+Se.up[t,1],bse[t]+Se.up[t,2])
Sp.mat[g,t]<-rbeta(1,asp[t]+Sp.up[t,1],bsp[t]+Sp.up[t,2])
}

Se<-Se.mat[g,]
Sp<-Sp.mat[g,]
}



######################
# Sample true status Y

p<-exp(X%*%b0)/(1+exp(X%*%b0))

zm<-dim(Z)[2]
ym<-dim(Y)[2]
w<-rep(0,N)
u<-runif(N)

res<-.C("samplatentstatus_",
            as.integer(N),
            as.integer(J),
            as.integer(zm),
            as.integer(ym),
            p,
            as.integer(Y),
            as.integer(Z),
            as.integer(w),
            u,
            Se,
            Sp,
            as.integer(na))
Y[,1]<-res[[8]]

#print(g)

if(visual==TRUE){
if(sum(seq(100,GI,step)==g)==1){
par(mfrow=c(length(b0),1))
for(h in 1:length(b0)){
plot(1:g,beta[1:g,h],type="l")
lines(1:g,cumsum(beta[1:g,h])/(1:g),col="red")
}}
}


} # End Gibbs sampling 


if(est.error==TRUE){
return(list("beta"=beta, "Se"=Se.mat, "Sp"=Sp.mat))
}
if(est.error==FALSE){
return(list("beta"=beta))
}

} # End of function



##########################################################################
# Note to run code need package mnormt
library(mnormt)
library(matrixcalc)


#############################################################################################
# Carries out Metropolis hasting sampling as discussed in Gamerman (1997)
WLS.Metropolis.Hastings.ip<-function(b0,X,y,a,R){

res<-mt_ct.ip(b=b0, y=y, X=X, a=a, R=R, sig=0)
Mb0<-res$M
Cb0<-res$C
Lb0<-res$L
sig<-res$s


b.star<-as.vector(rmnorm(1,Mb0,Cb0))
qb.star<-dmnorm(b.star,as.vector(Mb0),Cb0)

res<-mt_ct.ip(b=b.star, y=y, X=X, a=a, R=R, sig=sig)
Mb.star<-res$M
Cb.star<-res$C
Lb.star<-res$L

qb0<-dmnorm(b0,as.vector(Mb.star),Cb.star)

ratio.pi<-exp(Lb.star-Lb0)

alpha<-min(1,max(0,((ratio.pi*qb0)/(qb.star))))
accept<-rbinom(1,1,alpha)
if(accept==1){b1<-b.star}
if(accept==0){b1<-b0}
return("b1"=b1)
}

############################################
# Finds the Mean and Covariance of proposal 
# density according to Gammerman (1997)
mt_ct.ip<-function(b,y,X,a,R,sig){
P<-length(b)
R.inv<-solve(R)
xb<-X%*%b

p<-exp(xb)/(1+exp(xb))

W<-p*(1-p)
Yb<-xb*(p*(1-p))+(y-p)
XWX<-t(X*matrix(W,ncol=P, nrow=length(W),byrow=FALSE))%*%X
if(is.non.singular.matrix(XWX)==FALSE || sig==1){
XWX<-XWX+0.1*diag(rep(1,P))
sig<-1
}

C<-solve(R.inv+XWX)
C<-(C+t(C))/2
M<-C%*%(R.inv%*%a+t(X)%*%(Yb))
L<- (-1/2)*((b-a)%*%R.inv%*%(b-a)) + sum(y*(xb)-log(1+exp(xb)))  

return(list("M"=M, "C"=C, "L"=L, "s"=sig))
}