==1845835== Memcheck, a memory error detector
==1845835== Copyright (C) 2002-2017, and GNU GPL'd, by Julian Seward et al.
==1845835== Using Valgrind-3.18.1 and LibVEX; rerun with -h for copyright info
==1845835== Command: /data/blackswan/ripley/R/R-devel-vg/bin/exec/R --vanilla
==1845835== 

R Under development (unstable) (2023-01-26 r83691) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> pkgname <- "glmmrMCML"
> source(file.path(R.home("share"), "R", "examples-header.R"))
> options(warn = 1)
> library('glmmrMCML')
Loading required package: glmmrBase
Loading required package: Matrix
> 
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
> cleanEx()
> nameEx("ModelMCML")
> ### * ModelMCML
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ModelMCML
> ### Title: Extension to the Model class to use Markov Chain Monte Carlo
> ###   Maximum Likelihood
> ### Aliases: ModelMCML
> 
> ### ** Examples
> 
> 
> ## ------------------------------------------------
> ## Method `ModelMCML$MCML`
> ## ------------------------------------------------
> 
> ## Not run: 
> ##D df <- nelder(~(cl(10)*t(5)) > ind(10))
> ##D df$int <- 0
> ##D df[df$cl > 5, 'int'] <- 1
> ##D # specify parameter values in the call for the data simulation below
> ##D des <- ModelMCML$new(
> ##D   covariance = list( formula = ~ (1|gr(cl)*ar1(t)),
> ##D                      parameters = c(0.25,0.7)), 
> ##D   mean = list(formula = ~ factor(t) + int - 1,
> ##D                         parameters = c(rep(0,5),0.2)),
> ##D   data = df,
> ##D   family = gaussian()
> ##D )
> ##D ysim <- des$sim_data() # simulate some data from the model
> ##D fit1 <- des$MCML(y = ysim,method="mcnr",usestan=FALSE)
> ##D #fits the models using Stan 
> ##D fit2 <- des$MCML(y = ysim,method="mcnr")
> ##D  #adds a simulated likelihood step after the MCEM algorithm
> ##D fit3 <- des$MCML(y = ysim, sim.lik.step = TRUE)  
> ##D 
> ##D  # we could use LA to find better starting values
> ##D fit4 <- des$LA(y=ysim)
> ##D # to set better starting values we can use the update_parameters function 
> ##D des$update_parameters(mean = fit4$coefficients$est[1:6], 
> ##D                      cov = fit4$coefficients$est[7:8])
> ##D fit5 <- des$MCML(y = ysim,method="mcnr") # it should converge much more quickly                    
> ## End(Not run)
> 
> ## ------------------------------------------------
> ## Method `ModelMCML$LA`
> ## ------------------------------------------------
> 
> ## Not run: 
> ##D df <- nelder(~(cl(10)*t(5)) > ind(10))
> ##D df$int <- 0
> ##D df[df$cl > 5, 'int'] <- 1
> ##D # specify parameter values in the call for the data simulation below
> ##D des <- ModelMCML$new(
> ##D   covariance = list( formula = ~ (1|gr(cl)*ar1(t)),
> ##D                      parameters = c(0.25,0.7)), 
> ##D   mean = list(formula = ~ factor(t) + int - 1,
> ##D                         parameters = c(rep(0,5),-0.2)),
> ##D   data = df,
> ##D   family = stats::binomial()
> ##D )
> ##D ysim <- des$sim_data() # simulate some data from the model
> ##D fit1 <- des$LA(y = ysim)
> ## End(Not run)
> 
> 
> 
> cleanEx()
> nameEx("aic_mcml")
> ### * aic_mcml
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: aic_mcml
> ### Title: Calculates the conditional Akaike Information Criterion for the
> ###   GLMM
> ### Aliases: aic_mcml
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx()
> nameEx("gen_u_samples")
> ### * gen_u_samples
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: gen_u_samples
> ### Title: Generate samples of random effects using MCMC
> ### Aliases: gen_u_samples
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ## small example with simulated data
> ##D df <- nelder(~(j(10) * t(3)) > i(5))
> ##D des <- ModelMCML$new(
> ##D  covariance = list(
> ##D   formula =  ~(1|gr(j)*ar1(t)),
> ##D   parameters = c(0.25,0.7)
> ##D ),
> ##D mean = list(
> ##D   formula = ~factor(t)-1,
> ##D   parameters = rnorm(3)
> ##D ),
> ##D data=df,
> ##D family=gaussian()
> ##D )
> ##D ## simulate data
> ##D y <- des$sim_data()
> ##D ## get covariance definition matrix
> ##D ddata <- des$covariance$get_D_data()
> ##D ## simulate some values of the random effects
> ##D ## first, we need to extract the Cholesky decomposition of the covariance matrix D
> ##D L <- des$covariance$get_chol_D()
> ##D samp <- gen_u_samples(y=y,
> ##D   Z = as.matrix(des$covariance$Z),
> ##D   L = as.matrix(L),
> ##D   X = as.matrix(des$mean_function$X),
> ##D   beta = des$mean_function$parameters,,
> ##D   family = des$mean_function$family
> ##D  )
> ## End(Not run)
> 
> 
> 
> cleanEx()
> nameEx("glmmrMCML-package")
> ### * glmmrMCML-package
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: glmmrMCML-package
> ### Title: Markov Chain Monte Carlo Maximum Likelihood for Generalised
> ###   Linear Mixed Models
> ### Aliases: glmmrMCML-package glmmrMCML
> ### Keywords: package
> 
> ### ** Examples
> 
> ## Not run
> 
> 
> 
> cleanEx()
> nameEx("mcmc_sample")
> ### * mcmc_sample
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcmc_sample
> ### Title: Hamiltonian Monte Carlo Sampler for Model Random Effects
> ### Aliases: mcmc_sample
> 
> ### ** Examples
> 
> ## small example with simulated data
> df <- nelder(~(j(10) * t(3)) > i(5))
> des <- ModelMCML$new(
+  covariance = list(
+   formula =  ~(1|gr(j)*ar1(t)),
+   parameters = c(0.25,0.7)
+ ),
+ mean = list(
+   formula = ~factor(t)-1,
+   parameters = rnorm(3)
+ ),
+ data=df,
+ family=gaussian()
+ )
changes found, updating Z
changes found, updating
> ## simulate data
> y <- des$sim_data()
> ## get covariance definition matrix
> ddata <- des$covariance$get_D_data()
> ## simulate some values of the random effects
> ## first, we need to extract the Cholesky decomposition of the covariance matrix D
> L <- des$covariance$get_chol_D()
> ## generate samples using HMC
> mat <- mcmc_sample(Z = as.matrix(des$covariance$Z),
+    L = as.matrix(L),
+    X = as.matrix(des$mean_function$X),
+    y = as.vector(y),
+    beta = des$mean_function$parameters,
+    var_par = 1,
+    family = des$mean_function$family[[1]],
+    link = des$mean_function$family[[2]],
+    warmup = 250, 
+    nsamp = 250,
+    lambda = 5,
+    maxsteps = 100,
+    trace=1,
+    target_accept = 0.95)
 
 STARTING SAMPLING

MCMC Sampling
Warmup: Iter 0 of 500
Sampling: Iter 250 of 500
Accept rate: 1 steps: 94 step size: 0.0530929
----------------------------------------> 
> 
> 
> cleanEx()
> nameEx("mcml_full")
> ### * mcml_full
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_full
> ### Title: Markov Chain Monte Carlo Maximum Likelihood Algorithm
> ### Aliases: mcml_full
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ## small example with simulated data
> ##D ## create data and model object with 
> ##D ## parameters to simulate data and 
> ##D ## act as starting values
> ##D df <- nelder(~(j(10) * t(3)) > i(5))
> ##D des <- ModelMCML$new(
> ##D  covariance = list(
> ##D   formula =  ~(1|gr(j)*ar1(t)),
> ##D   parameters = c(0.25,0.7)
> ##D ),
> ##D mean = list(
> ##D   formula = ~factor(t)-1,
> ##D   parameters = rnorm(3)
> ##D ),
> ##D data=df,
> ##D family=gaussian()
> ##D )
> ##D ## simulate data
> ##D y <- des$sim_data()
> ##D ## get covariance definition matrix
> ##D ddata <- des$covariance$get_D_data()
> ##D ## fit model using default arguments for HMC sampler
> ##D out <- mcml_full(cov=ddata$cov,
> ##D    data=ddata$data,
> ##D    eff_range = rep(0,30),
> ##D    Z = as.matrix(des$covariance$Z),
> ##D    X = as.matrix(des$mean_function$X),
> ##D    y = y,
> ##D    family = des$mean_function$family[[1]],
> ##D    link=des$mean_function$family[[2]],
> ##D    start = c(rep(0.01,3),0.25,0.7,1.0),
> ##D    mcnr = TRUE)
> ## End(Not run)
> 
> 
> 
> cleanEx()
> nameEx("mcml_hess")
> ### * mcml_hess
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_hess
> ### Title: Generate Hessian matrix of GLMM
> ### Aliases: mcml_hess
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx()
> nameEx("mcml_hess_sparse")
> ### * mcml_hess_sparse
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_hess_sparse
> ### Title: Generate Hessian matrix of GLMM using sparse matrix methods
> ### Aliases: mcml_hess_sparse
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx()
> nameEx("mcml_la")
> ### * mcml_la
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_la
> ### Title: Maximum Likelihood with Laplace Approximation and Derivative
> ###   Free Optimisation
> ### Aliases: mcml_la
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx()
> nameEx("mcml_la_nr")
> ### * mcml_la_nr
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_la_nr
> ### Title: Maximum Likelihood with Laplace Approximation and Newton-Raphson
> ### Aliases: mcml_la_nr
> 
> ### ** Examples
> 
> ## small example with simulated data
> ## create data and model object with 
> ## parameters to simulate data and 
> ## act as starting values
> df <- nelder(~(j(10) * t(3)) > i(5))
> des <- ModelMCML$new(
+  covariance = list(
+   formula =  ~(1|gr(j)*ar1(t)),
+   parameters = c(0.25,0.7)
+ ),
+ mean = list(
+   formula = ~factor(t)-1,
+   parameters = rnorm(3)
+ ),
+ data=df,
+ family=gaussian()
+ )
changes found, updating Z
changes found, updating
> ## simulate data
> y <- des$sim_data()
> ## get covariance definition matrix
> ddata <- des$covariance$get_D_data()
> out3 <- mcml_la_nr(cov=ddata$cov,
+    data=ddata$data,
+    eff_range = rep(0,30),
+    Z = as.matrix(des$covariance$Z),
+    X = as.matrix(des$mean_function$X),
+    y = y,
+    family = des$mean_function$family[[1]],
+    link=des$mean_function$family[[2]],
+    start = c(rep(0.01,5),0.25,0.7,1.0),
+    usehess = FALSE,
+    tol=1e-2,verbose=FALSE,trace=0)
==1845835== Conditional jump or move depends on uninitialised value(s)
==1845835==    at 0x1FD9E918: mcml_la_nr(Eigen::Array<int, -1, -1, 0, -1, -1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Array<double, -1, 1, 0, -1, 1>, bool, double, bool, int, int) (packages/tests-vg/glmmrMCML/src/mcml_la.cpp:315)
==1845835==    by 0x1FD52C7B: _glmmrMCML_mcml_la_nr (packages/tests-vg/glmmrMCML/src/RcppExports.cpp:114)
==1845835==    by 0x4A2CF9: R_doDotCall (svn/R-devel/src/main/dotcode.c:927)
==1845835==    by 0x4E1439: bcEval (svn/R-devel/src/main/eval.c:7947)
==1845835==    by 0x4F4C2F: Rf_eval (svn/R-devel/src/main/eval.c:1013)
==1845835==    by 0x4F658D: R_execClosure (svn/R-devel/src/main/eval.c:2183)
==1845835==    by 0x4F72E2: Rf_applyClosure (svn/R-devel/src/main/eval.c:2109)
==1845835==    by 0x4F4D3B: Rf_eval (svn/R-devel/src/main/eval.c:1136)
==1845835==    by 0x4F8E39: do_set (svn/R-devel/src/main/eval.c:3246)
==1845835==    by 0x4F4FB4: Rf_eval (svn/R-devel/src/main/eval.c:1088)
==1845835==    by 0x52972C: Rf_ReplIteration (svn/R-devel/src/main/main.c:262)
==1845835==    by 0x529A67: R_ReplConsole (svn/R-devel/src/main/main.c:314)
==1845835==  Uninitialised value was created by a heap allocation
==1845835==    at 0x484086F: malloc (/builddir/build/BUILD/valgrind-3.18.1/coregrind/m_replacemalloc/vg_replace_malloc.c:381)
==1845835==    by 0x1FD9FB5E: aligned_malloc (R-devel/site-library/RcppEigen/include/Eigen/src/Core/util/Memory.h:159)
==1845835==    by 0x1FD9FB5E: conditional_aligned_malloc<true> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/util/Memory.h:214)
==1845835==    by 0x1FD9FB5E: conditional_aligned_new_auto<double, true> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/util/Memory.h:374)
==1845835==    by 0x1FD9FB5E: resize (R-devel/site-library/RcppEigen/include/Eigen/src/Core/DenseStorage.h:557)
==1845835==    by 0x1FD9FB5E: resize (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:319)
==1845835==    by 0x1FD9FB5E: _init1<long int> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:782)
==1845835==    by 0x1FD9FB5E: Array<long int> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/Array.h:171)
==1845835==    by 0x1FD9FB5E: mcml_la_nr(Eigen::Array<int, -1, -1, 0, -1, -1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Array<double, -1, 1, 0, -1, 1>, bool, double, bool, int, int) (packages/tests-vg/glmmrMCML/src/mcml_la.cpp:280)
==1845835==    by 0x1FD52C7B: _glmmrMCML_mcml_la_nr (packages/tests-vg/glmmrMCML/src/RcppExports.cpp:114)
==1845835==    by 0x4A2CF9: R_doDotCall (svn/R-devel/src/main/dotcode.c:927)
==1845835==    by 0x4E1439: bcEval (svn/R-devel/src/main/eval.c:7947)
==1845835==    by 0x4F4C2F: Rf_eval (svn/R-devel/src/main/eval.c:1013)
==1845835==    by 0x4F658D: R_execClosure (svn/R-devel/src/main/eval.c:2183)
==1845835==    by 0x4F72E2: Rf_applyClosure (svn/R-devel/src/main/eval.c:2109)
==1845835==    by 0x4F4D3B: Rf_eval (svn/R-devel/src/main/eval.c:1136)
==1845835==    by 0x4F8E39: do_set (svn/R-devel/src/main/eval.c:3246)
==1845835==    by 0x4F4FB4: Rf_eval (svn/R-devel/src/main/eval.c:1088)
==1845835==    by 0x52972C: Rf_ReplIteration (svn/R-devel/src/main/main.c:262)
==1845835== 
==1845835== Conditional jump or move depends on uninitialised value(s)
==1845835==    at 0x1FD9E942: mcml_la_nr(Eigen::Array<int, -1, -1, 0, -1, -1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Array<double, -1, 1, 0, -1, 1>, bool, double, bool, int, int) (packages/tests-vg/glmmrMCML/src/mcml_la.cpp:293)
==1845835==    by 0x1FD52C7B: _glmmrMCML_mcml_la_nr (packages/tests-vg/glmmrMCML/src/RcppExports.cpp:114)
==1845835==    by 0x4A2CF9: R_doDotCall (svn/R-devel/src/main/dotcode.c:927)
==1845835==    by 0x4E1439: bcEval (svn/R-devel/src/main/eval.c:7947)
==1845835==    by 0x4F4C2F: Rf_eval (svn/R-devel/src/main/eval.c:1013)
==1845835==    by 0x4F658D: R_execClosure (svn/R-devel/src/main/eval.c:2183)
==1845835==    by 0x4F72E2: Rf_applyClosure (svn/R-devel/src/main/eval.c:2109)
==1845835==    by 0x4F4D3B: Rf_eval (svn/R-devel/src/main/eval.c:1136)
==1845835==    by 0x4F8E39: do_set (svn/R-devel/src/main/eval.c:3246)
==1845835==    by 0x4F4FB4: Rf_eval (svn/R-devel/src/main/eval.c:1088)
==1845835==    by 0x52972C: Rf_ReplIteration (svn/R-devel/src/main/main.c:262)
==1845835==    by 0x529A67: R_ReplConsole (svn/R-devel/src/main/main.c:314)
==1845835==  Uninitialised value was created by a heap allocation
==1845835==    at 0x484086F: malloc (/builddir/build/BUILD/valgrind-3.18.1/coregrind/m_replacemalloc/vg_replace_malloc.c:381)
==1845835==    by 0x1FD9FB5E: aligned_malloc (R-devel/site-library/RcppEigen/include/Eigen/src/Core/util/Memory.h:159)
==1845835==    by 0x1FD9FB5E: conditional_aligned_malloc<true> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/util/Memory.h:214)
==1845835==    by 0x1FD9FB5E: conditional_aligned_new_auto<double, true> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/util/Memory.h:374)
==1845835==    by 0x1FD9FB5E: resize (R-devel/site-library/RcppEigen/include/Eigen/src/Core/DenseStorage.h:557)
==1845835==    by 0x1FD9FB5E: resize (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:319)
==1845835==    by 0x1FD9FB5E: _init1<long int> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:782)
==1845835==    by 0x1FD9FB5E: Array<long int> (R-devel/site-library/RcppEigen/include/Eigen/src/Core/Array.h:171)
==1845835==    by 0x1FD9FB5E: mcml_la_nr(Eigen::Array<int, -1, -1, 0, -1, -1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Array<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Array<double, -1, 1, 0, -1, 1>, bool, double, bool, int, int) (packages/tests-vg/glmmrMCML/src/mcml_la.cpp:280)
==1845835==    by 0x1FD52C7B: _glmmrMCML_mcml_la_nr (packages/tests-vg/glmmrMCML/src/RcppExports.cpp:114)
==1845835==    by 0x4A2CF9: R_doDotCall (svn/R-devel/src/main/dotcode.c:927)
==1845835==    by 0x4E1439: bcEval (svn/R-devel/src/main/eval.c:7947)
==1845835==    by 0x4F4C2F: Rf_eval (svn/R-devel/src/main/eval.c:1013)
==1845835==    by 0x4F658D: R_execClosure (svn/R-devel/src/main/eval.c:2183)
==1845835==    by 0x4F72E2: Rf_applyClosure (svn/R-devel/src/main/eval.c:2109)
==1845835==    by 0x4F4D3B: Rf_eval (svn/R-devel/src/main/eval.c:1136)
==1845835==    by 0x4F8E39: do_set (svn/R-devel/src/main/eval.c:3246)
==1845835==    by 0x4F4FB4: Rf_eval (svn/R-devel/src/main/eval.c:1088)
==1845835==    by 0x52972C: Rf_ReplIteration (svn/R-devel/src/main/main.c:262)
==1845835== 
> 
> 
> 
> cleanEx()
> nameEx("mcml_optim")
> ### * mcml_optim
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_optim
> ### Title: Likelihood maximisation for the GLMM
> ### Aliases: mcml_optim
> 
> ### ** Examples
> 
> ## small example with simulated data
> ## create data and model object with 
> ## parameters to simulate data and 
> ## act as starting values
> df <- nelder(~(j(8) * t(3)) > i(4))
> des <- ModelMCML$new(
+  covariance = list(
+   formula =  ~(1|gr(j)*ar1(t)),
+   parameters = c(0.25,0.7)
+ ),
+ mean = list(
+   formula = ~factor(t)-1,
+   parameters = rnorm(3)
+ ),
+ data=df,
+ family=gaussian()
+ )
changes found, updating Z
changes found, updating
> ## simulate data
> y <- des$sim_data()
> ## get covariance definition matrix
> ddata <- des$covariance$get_D_data()
> ## simulate some values of the random effects
> ## first, we need to extract the Cholesky decomposition of the covariance matrix D
> L <- des$covariance$get_chol_D()
> ## generate samples using HMC
> mat <- mcmc_sample(Z = as.matrix(des$covariance$Z),
+    L = as.matrix(L),
+    X = as.matrix(des$mean_function$X),
+    y = as.vector(y),
+    beta = des$mean_function$parameters,
+    var_par = 1,
+    family = des$mean_function$family[[1]],
+    link = des$mean_function$family[[2]],
+    warmup = 250, 
+    nsamp = 250,
+    lambda = 5,
+    maxsteps = 100,
+    trace=1,
+    target_accept = 0.95)
 
 STARTING SAMPLING

MCMC Sampling
Warmup: Iter 0 of 500
Sampling: Iter 250 of 500
Accept rate: 0.998 steps: 94 step size: 0.053138
----------------------------------------> out <- mcml_optim(cov=ddata$cov,
+    data=ddata$data,
+    eff_range = rep(0,30),
+    Z = as.matrix(des$covariance$Z),
+    X = as.matrix(des$mean_function$X),
+    y = y,
+    u = mat,
+    family = des$mean_function$family[[1]],
+    link=des$mean_function$family[[2]],
+    start = c(rep(0.01,3),0.25,0.7,1.0),
+    trace=0,
+    mcnr = TRUE)
> 
> 
> 
> cleanEx()
> nameEx("mcml_optim_sparse")
> ### * mcml_optim_sparse
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_optim_sparse
> ### Title: Likelihood maximisation for the GLMM using sparse matrix methods
> ### Aliases: mcml_optim_sparse
> 
> ### ** Examples
> 
> ## small example with simulated data
> ## create data and model object with 
> ## parameters to simulate data and 
> ## act as starting values
> df <- nelder(~(j(8) * t(3)) > i(4))
> des <- ModelMCML$new(
+  covariance = list(
+   formula =  ~(1|gr(j)*ar1(t)),
+   parameters = c(0.25,0.7)
+ ),
+ mean = list(
+   formula = ~factor(t)-1,
+   parameters = rnorm(3)
+ ),
+ data=df,
+ family=gaussian()
+ )
changes found, updating Z
changes found, updating
> ## simulate data
> y <- des$sim_data()
> ## get covariance definition matrix
> ddata <- des$covariance$get_D_data()
> ## simulate some values of the random effects
> ## first, we need to extract the Cholesky decomposition of the covariance matrix D
> L <- des$covariance$get_chol_D()
> ## generate samples using HMC - note this step does not use the sparse matrix
> mat <- mcmc_sample(Z = as.matrix(des$covariance$Z),
+    L = as.matrix(L),
+    X = as.matrix(des$mean_function$X),
+    y = as.vector(y),
+    beta = des$mean_function$parameters,
+    var_par = 1,
+    family = des$mean_function$family[[1]],
+    link = des$mean_function$family[[2]],
+    warmup = 250, 
+    nsamp = 250,
+    lambda = 5,
+    maxsteps = 100,
+    trace=1,
+    target_accept = 0.95)
 
 STARTING SAMPLING

MCMC Sampling
Warmup: Iter 0 of 500
Sampling: Iter 250 of 500
Accept rate: 1 steps: 94 step size: 0.053138
----------------------------------------> ## the specification of the covariance above results in a sparse covariance matrix,
> ## so we can just extract the components as below.
> out <- mcml_optim_sparse(cov=ddata$cov,
+    data=ddata$data,
+    eff_range = rep(0,30),
+    Ap = des$covariance$D@p,
+    Ai = des$covariance$D@i,
+    Z = as.matrix(des$covariance$Z),
+    X = as.matrix(des$mean_function$X),
+    y = y,
+    u = mat,
+    family = des$mean_function$family[[1]],
+    link=des$mean_function$family[[2]],
+    start = c(rep(0.01,3),0.25,0.7,1.0),
+    trace=0,
+    mcnr = TRUE)
> 
> 
> 
> cleanEx()
> nameEx("mcml_simlik")
> ### * mcml_simlik
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_simlik
> ### Title: Simulated likelihood optimisation step for MCML
> ### Aliases: mcml_simlik
> 
> ### ** Examples
> 
> 
> 
> 
> cleanEx()
> nameEx("mcml_simlik_sparse")
> ### * mcml_simlik_sparse
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mcml_simlik_sparse
> ### Title: Simulated likelihood optimisation step for MCML using sparse
> ###   matrix methods
> ### Aliases: mcml_simlik_sparse
> 
> ### ** Examples
> 
> 
> 
> 
> ### * <FOOTER>
> ###
> cleanEx()
> options(digits = 7L)
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed:  3265.997 54.39 3305.055 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')
==1845835== 
==1845835== HEAP SUMMARY:
==1845835==     in use at exit: 201,379,200 bytes in 44,449 blocks
==1845835==   total heap usage: 2,184,101 allocs, 2,139,652 frees, 1,997,584,071 bytes allocated
==1845835== 
==1845835== 352 bytes in 1 blocks are possibly lost in loss record 108 of 2,677
==1845835==    at 0x4845464: calloc (/builddir/build/BUILD/valgrind-3.18.1/coregrind/m_replacemalloc/vg_replace_malloc.c:1328)
==1845835==    by 0x401364B: UnknownInlinedFun (/usr/src/debug/glibc-2.33-21.fc34.x86_64/elf/../include/rtld-malloc.h:44)
==1845835==    by 0x401364B: allocate_dtv (/usr/src/debug/glibc-2.33-21.fc34.x86_64/elf/../elf/dl-tls.c:366)
==1845835==    by 0x401364B: _dl_allocate_tls (/usr/src/debug/glibc-2.33-21.fc34.x86_64/elf/../elf/dl-tls.c:612)
==1845835==    by 0x53C9008: pthread_create@@GLIBC_2.2.5 (in /usr/lib64/libpthread-2.33.so)
==1845835==    by 0x5396283: ??? (in /usr/lib64/libgomp.so.1.0.0)
==1845835==    by 0x538D2D0: GOMP_parallel (in /usr/lib64/libgomp.so.1.0.0)
==1845835==    by 0x1FD831B3: void Eigen::internal::parallelize_gemm<true, Eigen::internal::gemm_functor<double, long, Eigen::internal::general_matrix_matrix_product<long, double, 0, false, double, 0, false, 0, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::internal::gemm_blocking_space<0, double, double, -1, -1, -1, 1, false> >, long>(Eigen::internal::gemm_functor<double, long, Eigen::internal::general_matrix_matrix_product<long, double, 0, false, double, 0, false, 0, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::internal::gemm_blocking_space<0, double, double, -1, -1, -1, 1, false> > const&, long, long, long, bool) (R-devel/site-library/RcppEigen/include/Eigen/src/Core/products/Parallelizer.h:135)
==1845835==    by 0x1FD83BB0: scaleAndAddTo<Eigen::Matrix<double, -1, -1> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/products/GeneralMatrixMatrix.h:487)
==1845835==    by 0x1FD83BB0: void Eigen::internal::generic_product_impl<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::DenseShape, Eigen::DenseShape, 8>::evalTo<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&) (R-devel/site-library/RcppEigen/include/Eigen/src/Core/products/GeneralMatrixMatrix.h:437)
==1845835==    by 0x1FD6753D: run (R-devel/site-library/RcppEigen/include/Eigen/src/Core/ProductEvaluators.h:148)
==1845835==    by 0x1FD6753D: call_assignment_no_alias<Eigen::Matrix<double, -1, -1>, Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0>, Eigen::internal::assign_op<double, double> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:836)
==1845835==    by 0x1FD6753D: _set_noalias<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:732)
==1845835==    by 0x1FD6753D: _init1<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0>, Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:818)
==1845835==    by 0x1FD6753D: Matrix<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/Matrix.h:294)
==1845835==    by 0x1FD6753D: void Eigen::internal::call_assignment<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0>, Eigen::internal::assign_op<double, double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1>&, Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0> const&, Eigen::internal::assign_op<double, double> const&, Eigen::internal::enable_if<Eigen::internal::evaluator_assume_aliasing<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0>, Eigen::internal::evaluator_traits<Eigen::Product<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 0> >::Shape>::value, void*>::type) [clone .isra.0] (R-devel/site-library/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:796)
==1845835==    by 0x1FD8CE21: call_assignment<Eigen::Matrix<double, -1, -1>, Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/AssignEvaluator.h:782)
==1845835==    by 0x1FD8CE21: _set<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/PlainObjectBase.h:714)
==1845835==    by 0x1FD8CE21: operator=<Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0> > (R-devel/site-library/RcppEigen/include/Eigen/src/Core/Matrix.h:225)
==1845835==    by 0x1FD8CE21: glmmr::mcmlModel::mcmlModel(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>*, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>*, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >) (packages/tests-vg/glmmrMCML/src/../inst/include/glmmrmcml/mcmlmodel.h:67)
==1845835==    by 0x1FD6A034: mcmc_sample(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, int, int, double, double, int, int, int, double) (packages/tests-vg/glmmrMCML/src/mcml_full.cpp:407)
==1845835==    by 0x1FD512B9: _glmmrMCML_mcmc_sample (packages/tests-vg/glmmrMCML/src/RcppExports.cpp:66)
==1845835==    by 0x4A2C9F: R_doDotCall (svn/R-devel/src/main/dotcode.c:939)
==1845835== 
==1845835== LEAK SUMMARY:
==1845835==    definitely lost: 0 bytes in 0 blocks
==1845835==    indirectly lost: 0 bytes in 0 blocks
==1845835==      possibly lost: 352 bytes in 1 blocks
==1845835==    still reachable: 201,378,848 bytes in 44,448 blocks
==1845835==         suppressed: 0 bytes in 0 blocks
==1845835== Reachable blocks (those to which a pointer was found) are not shown.
==1845835== To see them, rerun with: --leak-check=full --show-leak-kinds=all
==1845835== 
==1845835== For lists of detected and suppressed errors, rerun with: -s
==1845835== ERROR SUMMARY: 19 errors from 3 contexts (suppressed: 0 from 0)