R Under development (unstable) (2025-06-22 r88341 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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.

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.

> 
> library(INLAtools)
Loading required package: Matrix
> 
> ## first dim
> G1 <- sparseMatrix(
+     i = c(2, 3, 1, 4, 1, 4, 2, 3),
+     j = c(1, 1, 2, 2, 3, 3, 4, 4),
+ )
> R1 <- Diagonal(n = nrow(G1), x = colSums(G1)) - 0.9*G1
> R1
4 x 4 sparse Matrix of class "dgCMatrix"
                        
[1,]  2.0 -0.9 -0.9  .  
[2,] -0.9  2.0  .   -0.9
[3,] -0.9  .    2.0 -0.9
[4,]  .   -0.9 -0.9  2.0
> 
> ## 2nd dim
> R2 <- sparseMatrix(
+     i = c(1L, 1L, 2L, 2L, 2L, 3L, 3L),
+     j = c(1L, 2L, 1L, 2L, 3L, 3L, 2L),
+     x = c(2,-1, -1,3,-1, 4, -1))
> R2
3 x 3 sparse Matrix of class "dgCMatrix"
             
[1,]  2 -1  .
[2,] -1  3 -1
[3,]  . -1  4
> 
> (r <- 0.9)
[1] 0.9
> R3 <- sparseMatrix(
+     i = c(1:3, 1:2, 2:3),
+     j = c(1:3, 2:3, 1:2),
+     x = c(1, 1+r^2, 1, rep(-r, 4))
+ )/(1-r^2)
> R3
3 x 3 sparse Matrix of class "dgCMatrix"
                                  
[1,]  5.263158 -4.736842  .       
[2,] -4.736842  9.526316 -4.736842
[3,]  .        -4.736842  5.263158
> 
> c(n1 <- nrow(R1),
+   n2 <- nrow(R2),
+   n3 <- nrow(R3))
[1] 4 3 3
> 
> ### cgeneric models for each dim
> cg1 <- cgeneric(
+     model = "generic0", R = R1,
+     constr = FALSE, scale = FALSE,
+     param = c(1, 0.5))
> cg2 <- cgeneric(
+     model = "generic0", R = R2,
+     constr = FALSE, scale = FALSE,
+     param = c(1, NA)) ## fixed to 1
> cg3 <- cgeneric(
+     model = "generic0", R = R3,
+     constr = FALSE, scale = FALSE,
+     param = c(1, NA)) ## fixed to 1
> 
> ## kronecker 2x1
> cg21 <- kronecker(cg2, cg1)
> cg21Q <- prec(cg21, theta = 0)
> 
> ## compare
> R21 <- kronecker(R2, R1)
> all.equal(Sparse(R21),
+           Sparse(cg21Q))
[1] TRUE
> 
> ## kronecker 3x(2x1)
> cg321 <- kronecker(cg3, cg21)
> cg321Q <- prec(cg321, theta = 0)
> 
> R321 <- kronecker(R3, R21)
> all.equal(Sparse(R321),
+           Sparse(cg321Q))
[1] TRUE
> 
> if("INLA" %in% loadedNamespaces()) {
+ 
+     dataf <- as.data.frame(
+     expand.grid(i1 = 1:n1,
+                 i2 = 1:n2,
+                 i3 = 1:n3)
+     )
+     nrow(dataf)
+ 
+     ## additional indexes (combined i1 and i2), and 1:(n1*n2*n3)
+     dataf$i12 <- rep(1:(n1*n2), n3)
+     dataf$iii <- 1:nrow(dataf)
+ 
+     ## no data
+     dataf$y <- rep(NA, nrow(dataf))
+ 
+     ## use the group
+     m1f <- y ~ 0 +
+         f(i12, model = 'generic0', Cmatrix = R21,
+           group = i3,
+           control.group = list(
+               model = 'ar1',
+               hyper = list(rho = list(
+                                initial = log((1+r)/(1-r)),
+                                fixed = TRUE))))
+ 
+     hfix <- list(prec = list(initial = 10, fixed = TRUE))
+ 
+     fit1 <- INLA::inla(
+         formula = m1f,
+         data = dataf,
+         control.mode = list(theta = 0, fixed = TRUE),
+         control.family = list(hyper = hfix),
+         control.compute = list(config = TRUE)
+     )
+ 
+     Q1 <- prec(fit1)
+ 
+     all.equal(Sparse(R321),
+               Sparse(Q1))
+ 
+     ## 'fit' cg321
+     m321 <- y ~ 0 + f(iii, model=cg321)
+     fit2 <- INLA::inla(
+         formula = m321,
+         data = dataf,
+         control.mode = list(theta = 0, fixed = TRUE),
+         control.family = list(hyper = hfix),
+         control.compute = list(config = TRUE)
+     )
+ 
+ 
+     Q2 <- prec(fit2)
+ 
+     all.equal(Sparse(R321),
+               Sparse(Q2))
+ 
+ }
> 
> proc.time()
   user  system elapsed 
   1.39    0.21    1.59