require(mvtnorm)
library(mnormt)

  matrix.sqrt <- function(A)
{
   sva <- svd(A)
    if (min(sva$d)>=0)
       Asqrt <- sva$u%*%diag(sqrt(sva$d))%*%t(sva$v) # svd e decomposição espectral
    else
       stop("Matrix square root is not defined")
    return(Asqrt)
}



nlf<-function(x,betas,u){
resp<- exp(betas[1]+u[1]) + exp(betas[2])/(1+exp((x-exp(betas[3]))/exp(betas[4]))) + (betas[5]+u[2])*(x-50)
return(resp)
}


funcH<-function(x,u,betas){
ni<-length(x)
resp<-matrix(0,ni,2)
resp[,1]<- rep(-exp(betas[1]+u[1]),ni)
resp[,2]<- (-x+50)
return(-resp)
}


Ymed<-function(x,betas,u){
resp<- nlf(x,betas,u)- funcH(x,u,betas)%*%u
return(resp)
}


grafT<-function(y0,mu0,sigma0,Qc,muc,sigmac,nu){
 f<-function(u) dmnorm(y0,mu0,u^(-1)*sigma0)*(1-pmnorm(Qc,muc,u^(-1)*sigmac))*dgamma(u,nu/2,nu/2)
 resp <- integrate(Vectorize(f),0,Inf)$value
 return(resp)
}


dSL<-function(y,mu,sigma,nu){
 f<-function(u) nu*u^(nu-1)*dmnorm(y,mu,u^(-1)*sigma)
 resp <- integrate(Vectorize(f),0,1)$value
 return(resp)
}


acumSL<-function(y,mu,sigma,nu){
 f<-function(u) pmnorm(y,mu,u^(-1)*sigma)*nu*u^(nu-1)
 resp <- integrate(Vectorize(f),0,1)$value
 return(resp)
}


grafSL<-function(y0,mu0,sigma0,Qc,muc,sigmac,nu){
 f<-function(u) dmnorm(y0,mu0,u^(-1)*sigma0)*(1-pmnorm(Qc,muc,u^(-1)*sigmac))*nu*u^(nu-1)
 resp <- integrate(Vectorize(f),0,1)$value
 return(resp)
}


dCN<-function(y,mu,sigma,nu){
 resp <-  nu[1]*dmnorm(y,mu,nu[2]^(-1)*sigma)+(1-nu[1])*dmnorm(y,mu,sigma)
 return(resp)
}

acumCN<-function(y,mu,sigma,nu){
 resp <- nu[1]*pmnorm(y,mu,nu[2]^(-1)*sigma)+(1-nu[1])*pmnorm(y,mu,sigma)
 return(resp)
}

grafCN<-function(y0,mu0,sigma0,Qc,muc,sigmac,nu){
 resp <- nu[1]*dmnorm(y0,mu0,nu[2]^(-1)*sigma0)*(1-pmnorm(Qc,muc,nu[2]^(-1)*sigmac))+(1-nu[1])*dmnorm(y0,mu0,sigma0)*(1-pmnorm(Qc,muc,sigmac))
 return(resp)
}




################################################################################
#### Log Veross Censored-LMM
################################################################################


LogVerosCens<-function(bb,cc,y,x,nj,beta1,sigmae,D1,nu,type="Normal"){

## type="Normal"
## type="T"


  m<-length(nj)
  
  p<-dim(x)[2]
  q1=2

ver<-matrix(0,m,1)

    D1<-as.matrix(D1)
    beta1<-as.matrix(beta1)

if(type=="T") {
nu=ceiling(nu)
 for (j in 1:m ){
            
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          bi=bb[j,]
         
          
          if(sum(cc1)==0){
          
          y1<-y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1<-matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1<-funcH(x1,bi,beta1)    
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- Ymed(x1,beta1,bi)
          ver[j]<- (dmt(y1,gammai,Psi,df=nu))
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1=funcH(x1,bi,beta1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-Ymed(x1,beta1,bi)
          if(sum(cc1)==nj[j]){
          muc=Ymed(x1,beta1,bi)
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
          ver[j]<-1-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)
                            }
          else {
          muc<- gammai[cc1==1,]+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-gammai[cc1==0,])
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          ver[j]<- grafT(y1[cc1==0],gammai[cc1==0],Psi[cc1==0,cc1==0],y1[cc1==1],muc,Sc,nu)
          #Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          #Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          #ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
           # ver[j]<-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)*dmt(y1[cc1==0],gammai[cc1==0,],Psi[cc1==0,cc1==0],df=nu)
                }
                        }

} }

if(type=="SL") {
 for (j in 1:m ){
            nu=nu
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          bi=bb[j,]
         
          
          if(sum(cc1)==0){
          
          y1<-y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1<-matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1<-funcH(x1,bi,beta1)    
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- Ymed(x1,beta1,bi)
          ver[j]<- (dSL(y1,gammai,Psi,nu))
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1=funcH(x1,bi,beta1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-Ymed(x1,beta1,bi)
          if(sum(cc1)==nj[j]){
          muc=Ymed(x1,beta1,bi)
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
          ver[j]<-1-acumSL(ai,rep(0,sum(cc1)),Ri,nu)
                            }
          else {
          muc<- gammai[cc1==1,]+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-gammai[cc1==0,])
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          ver[j]<- grafSL(y1[cc1==0],gammai[cc1==0],Psi[cc1==0,cc1==0],y1[cc1==1],muc,Sc,nu)
          #Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          #Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          #ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
           # ver[j]<-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)*dmt(y1[cc1==0],gammai[cc1==0,],Psi[cc1==0,cc1==0],df=nu)
                }
                        }

} }

if(type=="CN") {
 for (j in 1:m ){
            
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          bi=bb[j,]
         
          
          if(sum(cc1)==0){
          
          y1<-y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1<-matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1<-funcH(x1,bi,beta1)    
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- Ymed(x1,beta1,bi)
          ver[j]<- (dCN(y1,gammai,Psi,nu))
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1=funcH(x1,bi,beta1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-Ymed(x1,beta1,bi)
          if(sum(cc1)==nj[j]){
          muc=Ymed(x1,beta1,bi)
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
          ver[j]<-1-acumCN(ai,rep(0,sum(cc1)),Ri,nu)
                            }
          else {
          muc<- gammai[cc1==1,]+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-gammai[cc1==0,])
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          ver[j]<- grafCN(y1[cc1==0],gammai[cc1==0],Psi[cc1==0,cc1==0],y1[cc1==1],muc,Sc,nu)
          #Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          #Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          #ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
           # ver[j]<-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)*dmt(y1[cc1==0],gammai[cc1==0,],Psi[cc1==0,cc1==0],df=nu)
                }
                        }

} }



if(type=="Normal") {
 for (j in 1:m ){
 
                    bi=bb[j,]
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          if(sum(cc1)==0){

          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1=funcH(x1,bi,beta1)    
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- Ymed(x1,beta1,bi)
          ver[j]<-dmnorm(as.vector(y1),gammai,Psi)
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ])
          z1=funcH(x1,bi,beta1)    
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-Ymed(x1,beta1,bi)
          if(sum(cc1)==nj[j]){
          muc= gammai
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
           ver[j]<- 1-pmnorm(as.vector(ai),rep(0,sum(cc1)),Ri)
                            }
          else {
          muc<- gammai[cc1==1,]+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-gammai[cc1==0,])
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
#          ver[j]<-pmnorm(as.vector(ai),rep(0,sum(cc1)),Ri)*dmnorm(as.vector(y1[cc1==0]),as.vector(gammai[cc1==0,]),Psi[cc1==0,cc1==0])
          ver[j]<-dmnorm(y1[cc1==0],gammai[cc1==0,],Psi[cc1==0,cc1==0])*(1-pmnorm(y1[cc1==1],muc,Sc))
                }
                        }

} }


return(ver)
}

#LogVerosCens(resu.simT$mean$bi,cens,y,x,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,resu.simT$mean$nu,"T")

# LogVerosCens(cens,y,x,z,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,3,"Normal")
 #  sum(log(LogVerosCens(cens,y,x,z,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,2500,"T")))
#sum(log(LogVerosCens(cens,y,x,z,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,3,"Normal")))


################################################################################
###                                 CRITERIOS
################################################################################



################################################################################
# CPO e DIC
################################################################################


criterios<-function(cc,y,x,nj,espac=50,cadeia,type="T"){
m<-length(nj)
if (type == "T"){
dd<-attach.bugs(cadeia,overwrite = FALSE)
ver<-LogVerosCens(cadeia$mean$bi,cc,y,x,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,cadeia$mean$nu,type=type)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i,,]
sigma2es<-dd$sigmae[i]
bis<-dd$bi[i,,]
nu1<-dd$nu[i]
Loglikaux[,k]<- LogVerosCens(bis,cc,y,x,nj,betass,sigma2es,Ds,nu1,type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}

Np<-length(c(betass,sigma2es,Ds[upper.tri(matrix.sqrt(Ds), diag = T)],nu1))
}


if (type == "SL"){
dd<-attach.bugs(cadeia,overwrite = FALSE)
ver<-LogVerosCens(cadeia$mean$bi,cc,y,x,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,cadeia$mean$nu,type=type)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i,,]
sigma2es<-dd$sigmae[i]
bis<-dd$bi[i,,]
nu1<-dd$nu[i]
Loglikaux[,k]<- LogVerosCens(bis,cc,y,x,nj,betass,sigma2es,Ds,nu1,type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}

Np<-length(c(betass,sigma2es,Ds[upper.tri(matrix.sqrt(Ds), diag = T)],nu1))
}

if (type == "CN"){
dd<-attach.bugs(cadeia,overwrite = FALSE)
ver<-LogVerosCens(cadeia$mean$bi,cc,y,x,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,c(cadeia$mean$nu,cadeia$mean$gamma),type=type)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i,,]
sigma2es<-dd$sigmae[i]
bis<-dd$bi[i,,]
nu1<-dd$nu[i]
gamma1<-dd$gamma[i]
Loglikaux[,k]<- LogVerosCens(bis,cc,y,x,nj,betass,sigma2es,Ds,c(nu1,gamma1),type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}

Np<-length(c(betass,sigma2es,Ds[upper.tri(matrix.sqrt(Ds), diag = T)],nu1,gamma1))
}

if (type == "Normal"){
dd<-attach.bugs(cadeia,overwrite = FALSE)
ver<-LogVerosCens(cadeia$mean$bi,cc,y,x,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,1,type=type)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i,,]
sigma2es<-dd$sigmae[i]
bis<-dd$bi[i,,]
Loglikaux[,k]<-LogVerosCens(bis,cc,y,x,nj,betass,sigma2es,Ds,1,type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}
Np<-length(c(betass,sigma2es,Ds[upper.tri(matrix.sqrt(Ds), diag = T)]))
}

CPO<-sum(log(1/(apply(CPOaux,1,mean))))
DIC<-mean(apply(-2*log(Loglikaux),2,sum))-sum(2*ver)
EAIC<- -2*sum(log(ver))+2*Np
EBIC<- -2*sum(log(ver))+log(length(nj))*Np
IL<- -log(1/(apply(CPOaux,1,mean)))+ apply(log(Loglikaux),1,mean)
return(list(CPO=CPO,DIC=DIC,EAIC=EAIC,EBIC=EBIC,IL=IL))
}               


