diff --git a/nimbleModel/R/indexRuleAll.R b/nimbleModel/R/indexRuleAll.R index f36c7c5..8ce54f1 100644 --- a/nimbleModel/R/indexRuleAll.R +++ b/nimbleModel/R/indexRuleAll.R @@ -97,7 +97,14 @@ indexRuleAll_setup <- function(toIndexExprList, unrolledSize <- unrolledIndicesEnv$outputSize if (length(unrolledResults) > 1) { - allIndices <- do.call("cbind", unrolledResults) + lists <- sapply(unrolledResults, is.list) + if(any(lists)) { # Cases like `y[i,n1[i]:n2[i]]`. We unroll these. + allIndices <- do.call("rbind", lapply(seq_along(unrolledResults[[1]]), \(i) + as.matrix(do.call("expand.grid", + lapply(unrolledResults, \(x) x[[i]]) + )) + )) + } else allIndices <- do.call("cbind", unrolledResults) } else { allIndices <- matrix(unrolledResults[[1]], nrow = unrolledSize) } diff --git a/nimbleModel/R/modelDecl.R b/nimbleModel/R/modelDecl.R index 78b092c..3875dab 100644 --- a/nimbleModel/R/modelDecl.R +++ b/nimbleModel/R/modelDecl.R @@ -34,7 +34,8 @@ modelDeclClass <- R6Class( context, sourceLineNumber, truncated = FALSE, - boundExprs = NULL) { + boundExprs = NULL, + warnNonConstantBlocks = FALSE) { context <<- context sourceLineNumber <<- sourceLineNumber code <<- code @@ -69,8 +70,8 @@ modelDeclClass <- R6Class( targetVarExpr <- targetExpr[[2]] targetNodeExpr <<- targetExpr indexedBlocks <- sapply(indexExpr, checkForIndexedIntervals, context) - if (any(indexedBlocks)) { - stop("Non-constant indexing found in `", safeDeparse(indexExpr[[which(indexedBlocks)[1]]]), "` in `", safeDeparse(targetNodeExpr), "`. Block indices must be constant.") + if (any(indexedBlocks) && warnNonConstantBlocks) { + messageIfVerbose("Non-constant indexing found in `", safeDeparse(indexExpr[[which(indexedBlocks)[1]]]), "` in `", safeDeparse(targetNodeExpr), "`. Graph queries will incorrectly show elements of multivariate nodes as individual nodes but other functionality should operate correctly.") } } else { # There is a transformation, possibly with a subscript. diff --git a/nimbleModel/R/modelDef.R b/nimbleModel/R/modelDef.R index 068caa0..9b051f5 100644 --- a/nimbleModel/R/modelDef.R +++ b/nimbleModel/R/modelDef.R @@ -192,7 +192,8 @@ modelDefClass <- R6Class( declInfo[[iAns]] <<- modelDeclClass$new( code[[i]], contexts[[contextID]], - lineNumber + lineNumber, + warnNonConstantBlocks = TRUE ) } if (code[[i]][[1]] == "for") { diff --git a/nimbleModel/R/nodeRules.R b/nimbleModel/R/nodeRules.R index 230acba..aeb778e 100644 --- a/nimbleModel/R/nodeRules.R +++ b/nimbleModel/R/nodeRules.R @@ -497,7 +497,7 @@ nodeRangeClass <- R6Class( indexVars <- paste0("idx", seq_len(sum(boolExternalIndexRanges & !multiSlot))) exprs[boolExternalIndexRanges & !multiSlot] <- lapply(indexVars, as.name) - forText <- paste(indexVars, "in", sv) + forText <- paste0("`", indexVars, "`", " in ", sv) result <- safeDeparse(do.call("call", c( list("[", nm), @@ -512,7 +512,7 @@ nodeRangeClass <- R6Class( forText <- "`" } result <- paste0("`", result, forText) - } + } else result <- paste0("`", result, "`") return(result) }, diff --git a/nimbleModel/R/originalIndexingRules.R b/nimbleModel/R/originalIndexingRules.R index d726403..051735a 100644 --- a/nimbleModel/R/originalIndexingRules.R +++ b/nimbleModel/R/originalIndexingRules.R @@ -38,11 +38,19 @@ originalIndexingRuleClass <- R6Class( # Produces a varRange, though it's not really a range for a variable # but rather a range for the indices. + # (2023-06-10, commit 30ede6) # Do not remove duplicates because in generation of `calcRange`s there # can be cases where we need duplicated values in order to have correct # number of logProbs. + # (2026-07-09) actually we need to remove duplicates for y[i,n1[i]:n2[i]] case + # as that has one 'node' per scalar element. + # It appears that the 2023-06-10 thinking was incorrect as it seems to have + # been based on wanting `calculate()` to return as many elements as were + # computed in any given deterministic calculation, e.g., y[1:2] <- foo() + # returning 2 elements, but we actually don't return deterministic results anyway. + # See test in line 1302 of test-nodeRules.R in commit 30ede6. apply = function(fromVarRange) { - graphRule$apply(fromVarRange, removeDuplicates = FALSE) + graphRule$apply(fromVarRange, removeDuplicates = TRUE) } ) ) diff --git a/nimbleModel/tests/testthat/test-modelDefClass.R b/nimbleModel/tests/testthat/test-modelDefClass.R index 3bfd2fc..fb8d16e 100644 --- a/nimbleModel/tests/testthat/test-modelDefClass.R +++ b/nimbleModel/tests/testthat/test-modelDefClass.R @@ -282,19 +282,3 @@ test_that("detection of duplicated declarations", { }) -test_that("detection of non-constant block indices", { - code <- quote({ - for(i in 1:2) { - x[i:5] ~ dmnorm(z[i:5], Q[i:5,i:5]) - } - }) - expect_error(m <- modelClass$new(code), "Non-constant indexing") - - code <- quote({ - for(i in 1:2) { - x[i,(i+1):5] ~ dmnorm(z[(i+1):5], Q[(i+1):5,(i+1):5]) - } - }) - expect_error(m <- modelClass$new(code), "Non-constant indexing") - -}) diff --git a/nimbleModel/tests/testthat/test-modelGraph.R b/nimbleModel/tests/testthat/test-modelGraph.R index e94ce40..63aeba2 100644 --- a/nimbleModel/tests/testthat/test-modelGraph.R +++ b/nimbleModel/tests/testthat/test-modelGraph.R @@ -1548,16 +1548,14 @@ test_that("basic hierarchical models", { }) test_that("mixed-length block dependences", { + # See issue 31 regarding getting scalar elements of nodes in such cases. code <- quote({ for(i in 1:3) y[i, n1[i]:n2[i]] ~ dmulti(p[n1[i]:n2[i]], 10) }) - # Issue #31 - expect_error({ - modelDef <- modelDefClass$new(code, constants = list(n1 = c(3,1,2), n2 = c(6,3,2))); - result <- getDependencies(modelDef, 'p[2]'); + modelDef <- modelDefClass$new(code, constants = list(n1 = c(3,1,2), n2 = c(6,3,2))) + result <- getDependencies(modelDef, 'p[2]') expect_equal(result[[1]]$indexRanges, list(newIndexRange(matrix(c(2,2,2,3,1,2,3,2), ncol = 2)))) - }) }) diff --git a/nimbleModel/tests/testthat/test-nimbleModel.R b/nimbleModel/tests/testthat/test-nimbleModel.R index 39b29d0..25edf2c 100644 --- a/nimbleModel/tests/testthat/test-nimbleModel.R +++ b/nimbleModel/tests/testthat/test-nimbleModel.R @@ -659,7 +659,7 @@ test_that("non-sequential indexing cases", { nr <- m$getNodes()[[1]] expect_true(inherits(nr$indexRanges[[1]], "indexRangeMatrixClass")) expect_identical(nr$numExternalIndexRanges, 1L) - expect_identical(nr$toChar(), "`y[idx1]`, for idx1 in c(2, 3, 5)") + expect_identical(nr$toChar(), "`y[idx1]`, for `idx1` in c(2, 3, 5)") code <- nimbleCode({ y[c(2,3,5)] ~ dmnorm(mu[1:3], pr[1:3,1:3]) @@ -686,22 +686,23 @@ test_that("non-sequential indexing cases", { nr <- m$getNodes()[[2]] expect_true(inherits(nr$indexRanges[[1]], "indexRangeMatrixClass")) expect_identical(nr$numExternalIndexRanges, 0L) - expect_identical(nr$toChar(), "y[c(2, 3, 5)]") + expect_identical(nr$toChar(), "`y[c(2, 3, 5)]`") - if(FALSE) { # No operator def for nimC: This was part of the call: y[i = nimC(2, 3, 5)] - code <- nimbleCode({ + # Cannot compile. No operator def for nimC: This was part of the call: y[i = nimC(2, 3, 5)] + code <- nimbleCode({ y[c(2,3,5)] <- x[1:3] + 1 - }) - mclass <- nimbleModel(code, inits = list(x = 1:3), - returnClass = TRUE) - m <- mclass$new() - cmclass <- nCompile(mclass) - cm <- cmclass$new() - m$calculate() - cm$calculate() - expect_identical(m$y, c(NA,2,3,NA,4)) - expect_identical(cm$y, c(NA,2,3,NA,4)) - } + }) + mclass <- nimbleModel(code, inits = list(x = 1:3), + returnClass = TRUE) + m <- mclass$new() + m$calculate() + expect_identical(m$y, c(NA,2,3,NA,4)) + expect_error({ + cmclass <- nCompile(mclass) + cm <- cmclass$new() + cm$calculate() + expect_identical(cm$y, c(NA,2,3,NA,4)) + }) }) @@ -1406,6 +1407,152 @@ test_that("isData", { }) +test_that("non-constant block indexing", { + # Changing dimensions. + code <- quote({ + for(i in 1:3) + y[i, n1[i]:n2[i]] ~ dmnorm(mu[n1[i]:n2[i]], pr[n1[i]:n2[i],n1[i]:n2[i]]) + }) + model <- nimbleModel(code, data=list(y=matrix(rnorm(7*6),7)), constants = list(mu=rep(0,6), pr=diag(6),n1 = c(3,1,2), n2 = c(6,3,3))) + expect_identical(model$getNodes('y',nodesAsChars=TRUE), + c(paste0("y[1, ", 3:6, "]"), + paste0("y[2, ", 1:3, "]"), + paste0("y[3, ", 2:3, "]"))) + + truth <- dmnorm_chol(model$y[1,3:6],rep(0,4),diag(4), log=TRUE)+ + dmnorm_chol(model$y[2,1:3],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[3,2:3],rep(0,2),diag(2), log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) # Formerly incorrect as discussed in issue 31. + + # All same dimensions, constant precision. + code <- quote({ + for(i in 1:3) + y[i, n1[i]:n2[i]] ~ dmnorm(mu[n1[i]:n2[i]], pr[1:3,1:3]) + }) + model <- nimbleModel(code, data=list(y=matrix(rnorm(7*6),7)), constants = list(mu=rep(0,6), pr=diag(6),n1 = c(3,1,2), n2 = c(5,3,4))) + expect_identical(model$getNodes('y',nodesAsChars=TRUE), + c(paste0("y[1, ", 3:5, "]"), + paste0("y[2, ", 1:3, "]"), + paste0("y[3, ", 2:4, "]"))) + truth <- dmnorm_chol(model$y[1,3:5],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[2,1:3],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[3,2:4],rep(0,3),diag(3), log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) + + # Various crazy cases. + code <- quote({ + for(i in 1:3) + y[n1[i]:n2[i], k[i], i] ~ dmnorm(mu[n1[i]:n2[i]], pr[1:3,1:3]) + }) + model <- nimbleModel(code, data=list(y=array(rnorm(5*5*3),c(5,5,3))), constants = list(k = c(2,4,5), mu=rep(0,6), pr=diag(6), n1 = c(3,1,2), n2 = c(5,3,4))) + expect_identical(sort(model$getNodes('y',nodesAsChars=TRUE)), + sort(c(paste0("y[", 3:5, ", 2, 1]"), + paste0("y[", 1:3, ", 4, 2]"), + paste0("y[", 2:4, ", 5, 3]")))) + truth <- dmnorm_chol(model$y[3:5,2,1],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[1:3,4,2],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[2:4,5,3],rep(0,3),diag(3), log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) + + code <- quote({ + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i]] ~ dmnorm(mu[n1[i]:n2[i]], pr[n1[i]:n2[i],n1[i]:n2[i]]) + }) + model <- nimbleModel(code, data=list(y=array(rnorm(7*3*2),c(3,2,7))), constants = list(mu=rep(0,6), pr=diag(6),n1 = c(3,1,2), n2 = c(6,3,3))) + expect_identical(sort(model$getNodes('y',nodesAsChars=TRUE)), + sort(c(paste0("y[1, 1, ", 3:6, "]"), + paste0("y[1, 2, ", 3:6, "]"), + paste0("y[2, 1, ", 1:3, "]"), + paste0("y[2, 2, ", 1:3, "]"), + paste0("y[3, 1, ", 2:3, "]"), + paste0("y[3, 2, ", 2:3, "]")))) + truth <- dmnorm_chol(model$y[1,1,3:6],rep(0,3),diag(4), log=TRUE)+ + dmnorm_chol(model$y[1,2,3:6],rep(0,3),diag(4), log=TRUE)+ + dmnorm_chol(model$y[2,1,1:3],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[2,2,1:3],rep(0,3),diag(3), log=TRUE)+ + dmnorm_chol(model$y[3,1,2:3],rep(0,3),diag(2), log=TRUE)+ + dmnorm_chol(model$y[3,2,2:3],rep(0,3),diag(2), log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) + + code <- quote({ + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i],n1[i]:n2[i]] ~ dwish(pr[1:4,1:4], df = 10) + }) + y <- array(0, c(3,2,6,6)) + n1 <- c(3,1,2); n2 <- c(6,4,5) + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i],n1[i]:n2[i]] <- rwish_chol(1, diag(4), 10) + model <- nimbleModel(code, data=list(y=y), constants =list( pr=diag(6),n1 = n1,n2=n2)) + expect_length(model$getNodes('y',nodesAsChars=TRUE), 6*4*4) + truth <- dwish_chol(model$y[1,1,3:6,3:6],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[1,2,3:6,3:6],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[2,1,1:4,1:4],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[2,2,1:4,1:4],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[3,1,2:5,2:5],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[3,2,2:5,2:5],diag(4),10, log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) + + code <- quote({ + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i],n1[i]:n2[i]] ~ dwish(pr[n1[i]:n2[i],n1[i]:n2[i]], df = 10) + }) + y <- array(0, c(3,2,6,6)) + n1 <- c(3,1,2); n2 <- c(6,4,5) + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i],n1[i]:n2[i]] <- rwish_chol(1, diag(4), 10, 10) + model <- nimbleModel(code, data=list(y=y), constants =list( pr=diag(6),n1 = n1,n2=n2)) + expect_length(model$getNodes('y',nodesAsChars=TRUE), 6*4*4) + truth <- dwish_chol(model$y[1,1,3:6,3:6],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[1,2,3:6,3:6],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[2,1,1:4,1:4],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[2,2,1:4,1:4],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[3,1,2:5,2:5],diag(4), 10, log=TRUE)+ + dwish_chol(model$y[3,2,2:5,2:5],diag(4),10, log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) + + code <- quote({ + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i],n3[j]:n4[j]] ~ dwish(pr[n1[i]:n2[i],n3[j]:n4[j]], df = 10) + }) + y <- array(0, c(3,2,6,7)) + n1 <- c(3,1,2); n2 <- c(6,4,5) + n3 <- c(2,4); n4 <- c(5,7) + for(i in 1:3) + for(j in 1:2) + y[i,j,n1[i]:n2[i],n3[j]:n4[j]] <- rwish_chol(1, diag(4), 10, 10) + model <- nimbleModel(code, data=list(y=y), constants =list( pr=diag(7),n1 = n1,n2=n2,n3=n3,n4=n4)) + expect_length(model$getNodes('y',nodesAsChars=TRUE), 6*4*4) + + ## SS case + code <- quote({ + for(i in 2:4) + y[i, n1[i]:n2[i]] ~ dmnorm(y[i-1,n1[i]:n2[i]], pr[1:3,1:3]) + }) + model <- nimbleModel(code, data=list(y=matrix(rnorm(4*5),4)), constants = list(pr=diag(3), n1 = c(1,3,1,2), n2 = c(1,5,3,4))) + expect_identical(sort(model$getNodes('y',nodesAsChars=TRUE)), + sort(c(paste0("y[2, ", 3:5, "]"), + paste0("y[3, ", 1:3, "]"), + paste0("y[4, ", 2:4, "]")))) + + truth <- dmnorm_chol(model$y[2,3:5],model$y[1,3:5],diag(3), log=TRUE)+ + dmnorm_chol(model$y[3,1:3],model$y[2,1:3],diag(3), log=TRUE)+ + dmnorm_chol(model$y[4,2:4],model$y[3,2:4],diag(3), log=TRUE) + model$calculate() + expect_identical(model$getLogProb('y'), truth) +}) + diff --git a/nimbleModel/tests/testthat/test-nodeRules.R b/nimbleModel/tests/testthat/test-nodeRules.R index dcc3c86..c1c8fa5 100644 --- a/nimbleModel/tests/testthat/test-nodeRules.R +++ b/nimbleModel/tests/testthat/test-nodeRules.R @@ -1358,17 +1358,17 @@ test_that("nodeRange::print works correctly", { LHSrule <- nodeRuleClass$new(LHS, 1, context_i) expect_identical(LHSrule$apply()$toChar(), - "`mu[1:3, idx1]`, for idx1 in 2:8") + "`mu[1:3, idx1]`, for `idx1` in 2:8") LHS <- quote(mu[1:3,j,i]) LHSrule <- nodeRuleClass$new(LHS, 1, context_ij) expect_identical(LHSrule$apply()$toChar(), - "`mu[1:3, idx1, idx2]`, for idx1 in 3:5, idx2 in 2:8") + "`mu[1:3, idx1, idx2]`, for `idx1` in 3:5, `idx2` in 2:8") LHS <- quote(mu[1:3,j,i-1]) LHSrule <- nodeRuleClass$new(LHS, 1, context_ij) expect_identical(LHSrule$apply()$toChar(), - "`mu[1:3, idx1, idx2]`, for idx1 in 3:5, idx2 in 1:7") + "`mu[1:3, idx1, idx2]`, for `idx1` in 3:5, `idx2` in 1:7") })