suppressMessages({library(lme4); library(Matrix)})
args <- commandArgs(TRUE)
L1 <- as.integer(args[1]); L2 <- as.integer(args[2]); nobs <- as.integer(args[3]); out <- args[4]
set.seed(1)
# Crossed random-effects design: factor f1 (L1 levels, the sparse "field"),
# f2 (L2 levels, crossing -> hubs). Observations randomly cross f1 x f2.
f1 <- factor(sample(L1, nobs, replace=TRUE), levels=1:L1)
f2 <- factor(sample(L2, nobs, replace=TRUE), levels=1:L2)
# Build random-effects design Z = [Z1 | Z2], system A = Z'Z + I  (lme4 factorizes this)
Z1 <- sparse.model.matrix(~ 0 + f1)
Z2 <- sparse.model.matrix(~ 0 + f2)
Z  <- cbind(Z1, Z2)
A  <- crossprod(Z) + Diagonal(L1 + L2) * 1.0   # SPD, crossed-effects structure
A  <- forceSymmetric(A, "L")
At <- as(tril(A), "TsparseMatrix")
n  <- nrow(A)
con <- file(out, "w")
writeLines("%%MatrixMarket matrix coordinate real symmetric", con)
writeLines(sprintf("%d %d %d", n, n, length(At@x)), con)
writeLines(sprintf("%d %d %.15e", At@i + 1L, At@j + 1L, At@x), con)
close(con)
cat(sprintf("wrote %s  n=%d nnz_lower=%d\n", out, n, length(At@x)))