Running 'BinBerTest.R' [8s]
  Comparing 'BinBerTest.Rout' to 'BinBerTest.Rout.save' ... OK
  Running 'binomfamtest.R' [178s]
  Running 'cfamilies.R' [1s]
  Comparing 'cfamilies.Rout' to 'cfamilies.Rout.save' ... OK
  Running 'coreTest.R' [175s]
  Running 'distRandtests.R' [2s]
  Comparing 'distRandtests.Rout' to 'distRandtests.Rout.save' ... OK
  Running 'elTest.R' [5s]
  Comparing 'elTest.Rout' to 'elTest.Rout.save' ... OK
  Running 'families.R' [154s]
  Running 'familiesFiniteDiffs.R' [2s]
  Comparing 'familiesFiniteDiffs.Rout' to 'familiesFiniteDiffs.Rout.save' ... OK
  Running 'matvecmult.R' [1s]
  Comparing 'matvecmult.Rout' to 'matvecmult.Rout.save' ... OK
  Running 'mcseTest.R' [155s]
  Running 'objfunTest.R' [6s]
  Comparing 'objfunTest.Rout' to 'objfunTest.Rout.save' ... OK
  Running 'salamFiniteDiffs.R' [9s]
  Comparing 'salamFiniteDiffs.Rout' to 'salamFiniteDiffs.Rout.save' ... OK
  Running 'testpiecesBH.R' [6s]
  Comparing 'testpiecesBH.Rout' to 'testpiecesBH.Rout.save' ... OK
  Running 'testt.R' [159s]
  Running 'weightsTest.R' [8s]
  Comparing 'weightsTest.Rout' to 'weightsTest.Rout.save' ... OK
Running the tests in 'tests/binomfamtest.R' failed.
Complete output:
  > #the cumulant for binomial coded in R
  > cumR<-function(eta,ntrials){
  +    if(eta <= 0){
  + 	out <- ntrials * log1p(exp(eta))
  + 	}
  + 	
  + 	if(eta>0){
  + 	out <- ntrials * (eta + log1p(exp(-eta)))
  + 	}
  + 	
  + 	out
  + }
  > 
  > 
  > library(glmm)
  Loading required package: trust
  Loading required package: mvtnorm
  Loading required package: Matrix
  Loading required package: parallel
  Loading required package: doParallel
  Loading required package: foreach
  Loading required package: iterators
  > 
  > #choose a positive value for eta and compare to 
  > #calc of cum for bernoulli
  > eta <- 1
  > ntrials <- 1
  > this <- cumR(eta, ntrials=ntrials)
  > that <- bernoulli.glmm()$cum(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > #choose a negative value for eta and compare to 
  > #calc of cum for bernoulli
  > eta <- -1
  > this <- cumR(eta, ntrials=ntrials)
  > that <- bernoulli.glmm()$cum(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > #########
  > #the first deriv of the cumulant for binom in R
  > cpR<-function(eta,ntrials){
  + 
  + 	denom <- 1 + exp(-eta)
  + 	out <- ntrials/denom
  + 	out
  + }
  > 
  > #choose a couple values for eta and compare to 
  > #calc of cp for bernoulli
  > eta <- 1
  > this <- cpR(eta, ntrials=ntrials)
  > that <- bernoulli.glmm()$cp(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta <- -1
  > this <- cpR(eta, ntrials=ntrials)
  > that <- bernoulli.glmm()$cp(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > #########
  > #the second deriv of the cumulant for binom in R
  > cppR<-function(eta,ntrials){
  + 
  + 	first <- 1+exp(-eta)
  + 	second <-  1+exp(eta)
  + 	out <- ntrials/(first*second)
  + 	out
  + }
  > 
  > #choose a couple values for eta and compare to 
  > #calc of cpp for bernoulli
  > eta <- 1
  > this <- cppR(eta, ntrials=ntrials)
  > that <- bernoulli.glmm()$cpp(eta)
  > all.equal(this,that)
  [1] TRUE
  > 
  > eta <- -1
  > this <- cppR(eta, ntrials=ntrials)
  > that <- bernoulli.glmm()$cpp(eta)
  > all.equal(this,that)
  [1] TRUE
  > 
  > #########
  > # check cumR for ntrials = 5 for both pos, neg eta
  > ntrials <- 5
  > eta <- -1
  > this <- cumR(eta, ntrials)
  > that <- ntrials * log(1+exp(eta))
  > all.equal(this, that)
  [1] TRUE
  > eta <- 1
  > this <- cumR(eta,ntrials)
  > that <- ntrials * (eta+log(1+exp(-eta)))
  > all.equal(this, that)
  [1] TRUE
  > 
  > # check cpR for ntrials = 5 for both pos, neg eta
  > ntrials <- 5
  > eta <- -1
  > this <- cpR(eta,ntrials)
  > that <- ntrials/(1+exp(-eta))
  > all.equal(this, that)
  [1] TRUE
  > eta <- 1
  > this <- cpR(eta,ntrials)
  > that <- ntrials/(1+exp(-eta))
  > all.equal(this, that)
  [1] TRUE
  > 
  > # check cppR for ntrials = 5 for both pos, neg eta
  > ntrials <- 5
  > eta <- -1
  > this <- cppR(eta, ntrials)
  > that <- ntrials/((1+exp(-eta))*(1+exp(eta)))
  > all.equal(this, that)
  [1] TRUE
  > eta <- 1
  > this <- cppR(eta, ntrials)
  > that <- ntrials/((1+exp(-eta))*(1+exp(eta)))
  > all.equal(this, that)
  [1] TRUE
  > 
  > #########
  > # check central finite differences for cumR, cpR, cppR
  > ntrials <- 5
  > eta <- 1
  > delta <- .0001
  > 
  > #going from value to deriv
  > this <- cpR(eta, ntrials) * delta
  > that <- cumR(eta + delta/2, ntrials) - cumR(eta - delta/2, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > #going from first to second deriv
  > this <- cppR(eta, ntrials) * delta
  > that <- cpR(eta + delta/2, ntrials) - cpR(eta - delta/2, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > 
  > # ntrials*cumulant for bern = cumulant for binomial in R
  > eta <- .8
  > ntrials <-2
  > this <- cumR(eta, ntrials)
  > that <- ntrials * bernoulli.glmm()$cum(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta <- -.8
  > this <- cumR(eta, ntrials)
  > that <- ntrials * bernoulli.glmm()$cum(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > ntrials <-5
  > this <- cumR(eta, ntrials)
  > that <- ntrials * bernoulli.glmm()$cum(eta)
  > all.equal(this, that)
  [1] TRUE
  > 
  > 
  > 
  > ##################  
  > #I am now confident in the R calcs for binom's cum and its 2 derivs
  > ##################  
  > 
  > 
  > ##################  
  > #Now I'll check the C calcs for bern against C calcs for bin
  > #void cum3(double *etain, int *neta, int *typein, int *ntrials, double *wts, double *cumout)
  > #bin = bern for cumulant, ntrials = 1
  > eta <- 1
  > ntrials <- 1
  > wts <- 1
  > this <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(1.1))
  > that <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(1.1))
  > all.equal(this$out, that$out)
  [1] TRUE
  > 
  > #bin = 2 bern for cumulant
  > eta <- 1
  > ntrials <- 2
  > bincum <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > berncum <- ntrials* .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > all.equal(bincum, berncum)
  [1] TRUE
  > doublecheck <- cumR(eta, ntrials)
  > all.equal(bincum, doublecheck)
  [1] TRUE
  > 
  > eta <- -1
  > ntrials <-2
  > bincum <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > berncum <- ntrials* .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > all.equal(bincum, berncum)
  [1] TRUE
  > doublecheck <- cumR(eta, ntrials)
  > all.equal(bincum, doublecheck)
  [1] TRUE
  > 
  > 
  > #compare cumulant's derivative (R vs C)
  > eta <- 1
  > ntrials <- 4
  > this <- .C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
  > that <- cpR(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > #compare cumulant's 2nd derivative  (R vs C)
  > eta <- 1
  > ntrials <- 4
  > this <- .C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
  > that <- cppR(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > #make sure cumulant still works for eta length greater than 1
  > neta<- 6
  > eta <- 1:neta
  > wts <- rep(1,neta)
  > ntrials <- rep(4, times=neta)
  > length(eta) == length(ntrials)
  [1] TRUE
  > length(eta) == neta
  [1] TRUE
  > length(wts) == neta
  [1] TRUE
  > this <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(neta), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > that <- 0
  > for(i in 1:neta){
  + 	that <- that + cumR(eta[i], ntrials[i])
  + }
  > all.equal(this, that)
  [1] TRUE
  > #no segfault errors and the same answer. Yay!
  > 
  > 
  > # one last thing to check about the cumulant calcs
  > eta<- 1
  > ntrials<- 1
  > wts <- 1
  > length(eta) == length(ntrials)
  [1] TRUE
  > length(wts) == length(eta)
  [1] TRUE
  > that <- bernoulli.glmm()$cum(eta)
  > this <- binomial.glmm()$cum(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > that <-  .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<--2
  > ntrials<-1
  > that <- bernoulli.glmm()$cum(eta)
  > this <- binomial.glmm()$cum(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<--2
  > ntrials<-2
  > that <- 2*bernoulli.glmm()$cum(eta)
  > this <- binomial.glmm()$cum(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > that <-  .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<-1
  > ntrials<-1
  > that <- bernoulli.glmm()$cp(eta)
  > this <- binomial.glmm()$cp(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<-1
  > ntrials<-3
  > that <- 3*bernoulli.glmm()$cp(eta)
  > this <- binomial.glmm()$cp(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > that <-  .C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<-1
  > ntrials<-1
  > that <- bernoulli.glmm()$cpp(eta)
  > this <- binomial.glmm()$cpp(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<-1
  > ntrials<-3
  > that <- 3*bernoulli.glmm()$cpp(eta)
  > this <- binomial.glmm()$cpp(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > that <-  .C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
  > all.equal(this, that)
  [1] TRUE
  > 
  > eta<--1
  > ntrials<-1
  > that <- bernoulli.glmm()$cpp(eta)
  > this <- binomial.glmm()$cpp(eta, ntrials)
  > all.equal(this, that)
  [1] TRUE
  > 
  > 
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.28    0.29    1.56 
Running the tests in 'tests/coreTest.R' failed.
Complete output:
  > library(glmm)
  Loading required package: trust
  Loading required package: mvtnorm
  Loading required package: Matrix
  Loading required package: parallel
  Loading required package: doParallel
  Loading required package: foreach
  Loading required package: iterators
  > data(BoothHobert)
  > clust <- makeCluster(2)
  > set.seed(1234)
  > out<-glmm(y~0+x1,list(y~0+z1),varcomps.names=c("z1"),data=BoothHobert,
  +           family.glmm=bernoulli.glmm,m=50,doPQL=FALSE,debug=TRUE, cluster=clust)
  > 
  > vars <- new.env(parent = emptyenv())
  > debug<-out$debug
  > vars$m1 <- debug$m1
  > m2 <- debug$m2
  > m3 <- debug$m3
  > vars$zeta <- 5
  > vars$cl <- clust
  > registerDoParallel(vars$cl)                   #making cluster usable with foreach
  > vars$no_cores <- length(vars$cl)
  > vars$umat<-debug$umat
  > vars$newm <- nrow(vars$umat)
  > vars$u.star<-debug$u.star
  > vars$mod.mcml<-out$mod.mcml
  > vars$nu.pql <- debug$nu.pql
  > D.star.inv <- Dstarnotsparse <- vars$D.star <- as.matrix(debug$D.star)
  > 
  > getEk<-glmm:::getEk
  > addVecs<-glmm:::addVecs
  > genRand<-glmm:::genRand
  > 
  > vars$family.glmm<-out$family.glmm
  > vars$ntrials<- rep(1, length(out$y))
  > beta.pql <- debug$beta.pql
  > 
  > if(is.null(out$weights)){
  +   wts <- rep(1, length(out$y))
  + } else{
  +   wts <- out$weights
  + }
  > 
  > vars$wts<-wts
  > 
  > simulate <- function(vars, Dstarnotsparse, m2, m3, beta.pql, D.star.inv){
  +   #generate m1 from t(0,D*)
  +   if(vars$m1>0) genData<-rmvt(ceiling(vars$m1/vars$no_cores),sigma=Dstarnotsparse,df=vars$zeta,type=c("shifted"))
  +   if(vars$m1==0) genData<-NULL		
  +   
  +   #generate m2 from N(u*,D*)
  +   if(m2>0) genData2<-genRand(vars$u.star,vars$D.star,ceiling(m2/vars$no_cores))
  +   if(m2==0) genData2<-NULL
  +   
  +   
  +   #generate m3 from N(u*,(Z'c''(Xbeta*+zu*)Z+D*^{-1})^-1)
  +   if(m3>0){
  +     Z=do.call(cbind,vars$mod.mcml$z)
  +     eta.star<-as.vector(vars$mod.mcml$x%*%beta.pql+Z%*%vars$u.star)
  +     if(vars$family.glmm$family.glmm=="bernoulli.glmm") {cdouble<-vars$family.glmm$cpp(eta.star)}
  +     if(vars$family.glmm$family.glmm=="poisson.glmm"){cdouble<-vars$family.glmm$cpp(eta.star)}
  +     if(vars$family.glmm$family.glmm=="binomial.glmm"){cdouble<-vars$family.glmm$cpp(eta.star, vars$ntrials)}
  +     #still a vector
  +     cdouble<-Diagonal(length(cdouble),cdouble)
  +     wtsmat <- diag(vars$wts)
  +     Sigmuh.inv<- t(Z)%*%cdouble%*%wtsmat%*%Z+D.star.inv
  +     Sigmuh<-solve(Sigmuh.inv)
  +     genData3<-genRand(vars$u.star,Sigmuh,ceiling(m3/vars$no_cores))
  +   }
  +   if(m3==0) genData3<-NULL
  +   
  +   #	#these are from distribution based on data
  +   #	if(distrib=="tee")genData<-genRand(sigma.gen,s.pql,mod.mcml$z,m1,distrib="tee",gamm)
  +   #	if(distrib=="normal")genData<-genRand(sigma.pql,s.pql,mod.mcml$z,m1,distrib="normal",gamm)
  +   #	#these are from standard normal
  +   #	ones<-rep(1,length(sigma.pql))
  +   #	zeros<-rep(0,length(s.pql))
  +   #	genData2<-genRand(ones,zeros,mod.mcml$z,m2,distrib="normal",gamm)
  +   
  +   umat<-rbind(genData,genData2,genData3)
  +   m <- nrow(umat)
  +   list(umat=umat, m=m, Sigmuh.inv=Sigmuh.inv)
  + }
  > 
  > clusterSetRNGStream(vars$cl, 1234)
  > 
  > clusterExport(vars$cl, c("vars", "Dstarnotsparse", "m2", "m3", "beta.pql", "D.star.inv", "simulate", "genRand"), envir = environment())     #installing variables on each core
  > noprint <- clusterEvalQ(vars$cl, umatparams <- simulate(vars=vars, Dstarnotsparse=Dstarnotsparse, m2=m2, m3=m3, beta.pql=beta.pql, D.star.inv=D.star.inv))
  > 
  > vars$nbeta <- 1
  > vars$p1=vars$p2=vars$p3=1/3
  > par<-c(6,1.5)
  > del<-rep(10^-8,2)
  > 
  > objfun<-glmm:::objfun
  > 
  > core2<-objfun(par=par, vars=vars)
  > 
  > umats <- clusterEvalQ(vars$cl, umatparams$umat)
  > umat <- Reduce(rbind, umats)
  > 
  > Sigmuh.invs <- clusterEvalQ(vars$cl, umatparams$Sigmuh.inv)
  > Sigmuh.inv <- Sigmuh.invs[[1]]
  > 
  > stopCluster(clust)
  > 
  > vars$cl <- makeCluster(1)
  > registerDoParallel(vars$cl)                   #making cluster usable with foreach
  > vars$no_cores <- length(vars$cl)
  > 
  > clusterExport(vars$cl, c("umat", "vars", "Sigmuh.inv"), envir = environment())
  > noprint <- clusterEvalQ(vars$cl, umatparams <- list(umat=umat, m=vars$newm, Sigmuh.inv=Sigmuh.inv))
  > 
  > core1<-objfun(par=par, vars=vars)
  > 
  > all.equal(core1, core2)
  [1] TRUE
  > 
  > stopCluster(vars$cl)
  > 
  > proc.time()
     user  system elapsed 
     3.07    0.34    6.64 
Running the tests in 'tests/families.R' failed.
Complete output:
  > library(glmm)
  Loading required package: trust
  Loading required package: mvtnorm
  Loading required package: Matrix
  Loading required package: parallel
  Loading required package: doParallel
  Loading required package: foreach
  Loading required package: iterators
  > theta<-1
  > delta<-1E-8
  > 
  > #check derivatives for bernoulli using finite differences
  > this<-bernoulli.glmm()$cp(theta)*delta
  > that<-bernoulli.glmm()$cum(theta+delta)-bernoulli.glmm()$cum(theta)
  > all.equal(this,that)
  [1] TRUE
  > 
  > this<-bernoulli.glmm()$cpp(theta)*delta
  > that<-bernoulli.glmm()$cp(theta+delta)-bernoulli.glmm()$cp(theta)
  > all.equal(this,that)
  [1] TRUE
  > 
  > #check derivatives for poisson using finite differences
  > this<-poisson.glmm()$cp(theta)*delta
  > that<-poisson.glmm()$cum(theta+delta)-poisson.glmm()$cum(theta)
  > all.equal(this,that)
  [1] TRUE
  > 
  > this<-poisson.glmm()$cpp(theta)*delta
  > that<-poisson.glmm()$cp(theta+delta)-poisson.glmm()$cp(theta)
  > all.equal(this,that)
  [1] TRUE
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.15    0.15    1.29 
Running the tests in 'tests/mcseTest.R' failed.
Complete output:
  > #some functions from glmm that we need
  > library(glmm)
  Loading required package: trust
  Loading required package: mvtnorm
  Loading required package: Matrix
  Loading required package: parallel
  Loading required package: doParallel
  Loading required package: foreach
  Loading required package: iterators
  > getFamily <- glmm:::getFamily
  > getEk<-glmm:::getEk
  > addVecs<-glmm:::addVecs
  > tconstant <- glmm:::tconstant
  > addMats<-glmm:::addMats
  > 
  > #other functions we'll need
  > 
  > #here is el, likelihood
  > elR <- function(Y,X,eta,family.mcml,wts){
  + 	family.mcml<-getFamily(family.mcml)
  + 	neta<-length(eta)
  + 	ntrials <- rep(1, neta)
  + 
  + 	if(family.mcml$family.glmm=="bernoulli.glmm"){
  + 		foo<-.C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), cumout=double(1))$cumout
  + 		mu<-.C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(1), ntrials=as.integer(ntrials), cpout=double(neta))$cpout
  + 		cdub<-.C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(1), ntrials=as.integer(ntrials), cppout=double(neta))$cppout
  + 	}
  + 	if(family.mcml$family.glmm=="poisson.glmm"){
  + 		foo<-.C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(2), ntrials=as.integer(ntrials), wts=as.double(wts),cumout=double(1))$cumout
  + 		mu<-.C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(2), ntrials=as.integer(ntrials), cpout=double(neta))$cpout
  + 		cdub<-.C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(neta), type=as.integer(2), ntrials=as.integer(ntrials), cppout=double(neta))$cppout
  + 	}
  + 
  + 	value<-as.numeric(Y%*%eta-foo)
  + 	gradient<-t(X)%*%(Y-mu)	
  + 	cdubmat<-diag(cdub)
  + 	hessian<-t(X)%*%(-cdubmat)%*%X
  + 	
  + 	list(value=value,gradient=gradient,hessian=hessian)
  + }
  > 
  > #for calc t
  > tdist2<-function(tconst,u, Dstarinv,zeta,myq){
  + 	inside<-1+t(u)%*%Dstarinv%*%u/zeta
  + 	logft<-tconst - ((zeta+myq)/2)*log(inside)
  + 	as.vector(logft)
  + }
  > 
  > 
  > distRandGeneral<-function(uvec,mu,Sigma.inv){
  + 	logDetSigmaInv<-sum(log(eigen(Sigma.inv,symmetric=TRUE)$values))
  + 	umu<-uvec-mu
  + 	piece2<-t(umu)%*%Sigma.inv%*%umu
  + 	out<-as.vector(.5*(logDetSigmaInv-piece2))
  + 	const<-length(uvec)*.5*log(2*pi)
  + 	out<-out-const
  + 	out
  + }
  > 
  > 
  > distRand <-
  + function(nu,U,z.list,mu){
  + 	# T=number variance components
  + 	T<-length(z.list)
  + 	
  + 	#nrandom is q_t
  + 	nrand<-lapply(z.list,ncol)
  + 	nrandom<-unlist(nrand)
  + 	totnrandom<-sum(nrandom)
  + 	
  + 	mu.list<-U.list<-NULL
  + 	if(T==1) {
  + 		U.list[[1]]<-U
  + 		mu.list[[1]]<-mu
  + 		}
  + 
  + 	if(T>1){
  + 		U.list[[1]]<-U[1:nrandom[1]] 
  + 		mu.list[[1]]<-mu[1:nrandom[1]]
  + 		for(t in 2:T){
  + 			thing1<-sum(nrandom[1:t-1])+1
  + 			thing2<-sum(nrandom[1:t])
  + 			U.list[[t]]<-U[thing1:thing2]
  + 			mu.list[[t]]<-mu[thing1:thing2]
  + 		}
  + 	}
  + 	
  + 	val<-gradient<-Hessian<-rep(0,T)
  + 	
  + 	#for each variance component
  + 	for(t in 1:T){
  + 		you<-as.vector(U.list[[t]])
  + 		mew<-as.vector(mu.list[[t]])
  + 		Umu<-(you-mew)%*%(you-mew)
  + 		val[t]<- -length(U)*log(2*pi)/2+as.numeric(-.5*nrandom[t]*log(nu[t])-Umu/(2*nu[t]))
  + 		
  + 		gradient[t]<- -nrandom[t]/(2*nu[t])+Umu/(2*(nu[t])^2)
  + 		
  + 		Hessian[t]<- nrandom[t]/(2*(nu[t])^2)- Umu/((nu[t])^3)
  + 		
  + 	}
  + 		
  + 	value<-sum(val)
  + 	if(T>1) hessian<-diag(Hessian)
  + 	if(T==1) hessian<-matrix(Hessian,nrow=1,ncol=1)
  + 	
  + 	list(value=value,gradient=gradient,hessian=hessian)		
  + }
  > 
  > mcseTEST <- function(mod){
  + 
  + 
  + 
  + 	beta <- mod$beta
  + 	nu <-mod$nu
  + 	umat<-mod$umat
  + 	m<-nrow(umat)
  + 
  + 	x <- mod$x
  + 	y <- mod$y
  + 	Z=do.call(cbind,mod$z)
  + 	T<-length(mod$z)
  + 	nrand<-lapply(mod$z,ncol)
  + 	nrandom<-unlist(nrand)
  + 	
  + 	if(is.null(mod$weights)){
  + 	  wts <- rep(1, length(y))
  + 	} else{
  + 	  wts <- mod$weights
  + 	}
  + 
  + 
  + 	family.glmm<-getFamily(mod$family.glmm)
  + 	if(family.glmm$family.glmm=="bernoulli.glmm"){family_glmm=1}	
  + 	if(family.glmm$family.glmm=="poisson.glmm"){family_glmm=2}	
  + 	if(family.glmm$family.glmm=="binomial.glmm"){family_glmm=3}	
  + 
  + 	pvec <-mod$pvec
  + 	p1 <- pvec[1]
  + 	p2 <- pvec[2]
  + 	p3 <- pvec[3]
  + 	zeta <- mod$zeta
  + 	beta.pql <- mod$beta.pql
  + 	nu.pql <- mod$nu.pql
  + 
  + 
  + 
  + 	#need to get D*
  + 	eek<-getEk(mod$z)
  + 	Aks<-Map("*",eek,nu.pql)
  + 	D.star<-addVecs(Aks) 
  + 	D.star<-diag(D.star)
  + 	D.star.inv<-solve(D.star)
  + 
  + 	#need to recreate the variance matrix of  imp sampling distribution
  + 	
  + 	eta.star<-as.vector(x%*%beta.pql+Z%*%mod$u.pql)
  + 	cdouble<-bernoulli.glmm()$cpp(eta.star) #still a vector
  + 	cdouble<-diag(cdouble)
  + 	Sigmuh.inv<- t(Z)%*%cdouble%*%Z+D.star.inv
  + 	Sigmuh<-solve(Sigmuh.inv)
  + 
  + 	Dstinvdiag<-diag(D.star.inv)
  + 
  + 	tconst<-tconstant(zeta,nrow(D.star.inv),Dstinvdiag)
  + 
  + 	eta<-b<-rep(0,m)
  + 	lfu<-lfyu<-list(rep(c(0,0,0),m))
  + 	lfu.twid<-matrix(data=NA,nrow=m,ncol=4)
  + 
  + 	#stuff for the numerator of v-hat
  + 	npar <- length(beta)+length(nu)
  + 	Gpiece<-rep(0,npar)
  + 	squaretop <- b
  + 	numpiece <- matrix(data=NA, nrow=m, ncol=npar*npar)
  + 	
  + 	#for each simulated random effect vector
  + 	for(k in 1:m){
  + 		Uk<-umat[k,]  #use the simulated vector as our random effect vec
  + 		eta<-mod$x%*%beta+Z%*%Uk # calculate eta using it
  + 		zeros<-rep(0,length(Uk))
  + 
  + 		#log f_theta(u_k)
  + 		lfu[[k]]<-distRand(nu,Uk,mod$z,zeros) 
  + 
  + 		#log f_theta(y|u_k)
  + 		lfyu[[k]]<-elR(mod$y,mod$x,eta,family.glmm, wts) 
  + 
  + 		#log f~_theta(u_k)
  + 		lfu.twid[k,1]<-tdist2(tconst,Uk,D.star.inv,zeta=zeta,myq=nrow(D.star.inv))
  + 		lfu.twid[k,2]<-distRandGeneral(Uk,mod$u.pql,D.star.inv)
  + 		lfu.twid[k,3]<-distRandGeneral(Uk,mod$u.pql,Sigmuh.inv)
  + 		
  + 		tempmax<-max(lfu.twid[k,1:3])
  + 		blah<-exp(lfu.twid[k,1:3]-tempmax)
  + 		pea<-c(p1,p2,p3)
  + 		qux<-pea%*%blah
  + 		lfu.twid[k,4]<-tempmax+log(qux)
  + 		
  + 		b[k]<-as.numeric(lfu[[k]]$value)+as.numeric(lfyu[[k]]$value)-lfu.twid[k,4]
  + 
  + 		#things for the numerator now
  + 		Gpiece <- c(lfyu[[k]]$gradient,lfu[[k]]$gradient)
  + 		GGT <- Gpiece %*% t(Gpiece)	
  + 
  + 		squaretop <- exp(	2*as.numeric(lfu[[k]]$value) + 2*as.numeric(lfyu[[k]]$value) - 2*lfu.twid[k,4] )
  + 	
  + 		numpiece[k,] <- as.numeric(GGT*squaretop)
  + 	}
  + 
  + 	#return to denominator of Vhat
  + 	#bk are unnormalized log weights
  + 	a<-max(b) #biggest logftheta(uk,y)-logftwiddle
  + 	thing<-exp(b-a) 
  + 	mygamma<- exp(a)*sum(thing)/m
  + 	Vdenom <- mygamma^2 	#denom of V hat
  + 	#now the denominator is finished
  + 
  + 	
  + 	#return to numerator of Vhat
  + 
  + 	num<-apply( numpiece   ,2,sum)/m
  + 	temp<-num/Vdenom
  + 	Vhat<- matrix(temp,nrow=npar,ncol=npar)
  + 	Vhat 
  + 
  + 	Uhat <- mod$loglike.hessian
  + 	Uhatinv <- qr.solve(Uhat)
  + 
  + 	UinvVUinv <- Uhatinv %*% Vhat %*% Uhatinv
  + 
  + 	MCSE <- sqrt(diag(UinvVUinv)/m)
  + 	MCSE
  + }
  > 
  > 
  > data(BoothHobert)
  > clust <- makeCluster(2)
  > set.seed(123)
  > mod.mcml1<-glmm(y~0+x1, list(y~0+z1), varcomps.names=c("z1"), data=BoothHobert, family.glmm=bernoulli.glmm, m=1000, doPQL=TRUE, cluster=clust)
  > all.equal(mcseTEST(mod.mcml1), as.numeric(mcse(mod.mcml1)))
  [1] TRUE
  > 
  > stopCluster(clust)
  > 
  > proc.time()
     user  system elapsed 
     3.85    0.46    6.81 
Running the tests in 'tests/testt.R' failed.
Complete output:
  > # test the t distribution
  > 
  > library(glmm)
  Loading required package: trust
  Loading required package: mvtnorm
  Loading required package: Matrix
  Loading required package: parallel
  Loading required package: doParallel
  Loading required package: foreach
  Loading required package: iterators
  > 
  > tdist2<-function(tconst,u, Dstarinv,zeta,myq){
  + 	inside<-1+t(u)%*%Dstarinv%*%u/zeta
  + 	logft<-tconst - ((zeta+myq)/2)*log(inside)
  + 	as.vector(logft)
  + }
  > 
  > zeta<-5
  > myq<-10
  > set.seed(1234)
  > u<-rt(10,df=zeta)
  > Dstarinvdiag<-rep(1.1,10)
  > tconstant<-glmm:::tconstant
  > tconst<-tconstant(zeta,myq,Dstarinvdiag)
  > 
  > Dstarinv<-diag(Dstarinvdiag)
  > tR<-tdist2(tconst,u,Dstarinv,zeta,myq)
  > 
  > #compare version in R to version in C
  > tC<-.C(glmm:::C_tdist,as.double(Dstarinv),  as.integer(myq), as.double(u), as.integer(zeta), as.double(tconst), double(1))
  > all.equal(tC[[6]],tR)
  [1] TRUE
  > 
  > #compare my version to the mvtnorm version
  > Dstar<-diag(10)/1.1
  > theirs<-dmvt(u,sigma=Dstar,df=zeta)
  > theirs==tR
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     1.17    0.14    1.29