diff --git a/packages/nimble/DESCRIPTION b/packages/nimble/DESCRIPTION
index 751574f43..94abd6c3e 100644
--- a/packages/nimble/DESCRIPTION
+++ b/packages/nimble/DESCRIPTION
@@ -50,6 +50,7 @@ Collate:
     genCpp_buildIntermediates.R
     genCpp_processSpecificCalls.R
     genCpp_sizeProcessing.R
+    genCpp_toEigenize.R
     genCpp_insertAssertions.R
     genCpp_maps.R
     genCpp_liftMaps.R
diff --git a/packages/nimble/R/RCfunction_compile.R b/packages/nimble/R/RCfunction_compile.R
index 81fbd8e30..85cdd9bad 100644
--- a/packages/nimble/R/RCfunction_compile.R
+++ b/packages/nimble/R/RCfunction_compile.R
@@ -265,10 +265,15 @@ RCfunProcessing <- setRefClass('RCfunProcessing',
                                            browser()
                                        }
 
+                                       if(nimbleOptions('experimentalNewSizeProcessing')) {
+                                           exprClasses_setToEigenize(compileInfo$nimExpr, compileInfo$newLocalSymTab, compileInfo$typeEnv)
+                                       }
+
                                        tryResult <- try(exprClasses_insertAssertions(compileInfo$nimExpr))
                                        if(inherits(tryResult, 'try-error')) {
                                            stop(paste('There is some problem at the insertAdditions processing step for this code:\n', paste(deparse(compileInfo$origRcode), collapse = '\n'), collapse = '\n'), call. = FALSE)
                                        }
+                                       
                                        if(debug) {
                                            print('compileInfo$nimExpr$show(showAssertions = TRUE)')
                                            compileInfo$nimExpr$show(showAssertions = TRUE)
diff --git a/packages/nimble/R/all_utils.R b/packages/nimble/R/all_utils.R
index 82dac2e1a..1f6333240 100644
--- a/packages/nimble/R/all_utils.R
+++ b/packages/nimble/R/all_utils.R
@@ -2,20 +2,48 @@
 
 double <- function(ndim, dims) {}
 
-labelFunctionCreator <- function(lead, start = 1) {
-  nextIndex <- start
-  force(lead)
-  labelGenerator <- function(reset = FALSE, count = 1, envName = "") {
-    if(reset) {
-      nextIndex <<- 1
-      return(invisible(NULL))
+## Sequential label generation system:
+## labelFunctionMetaCreator returns a function that returns a function.
+## labelFunctionMetaCreator is only called once, immediately below, to create labelFunctionCreator
+## The outer layer allows allLabelFunctionCreators to be in the closure of every function returned
+## by labelFunctionCreator.  Each of those functions is registered as an element of allLableFunctionCreators.
+## 
+## This scheme allows the function resetLabelFunctionCreators below to work simply,
+## resetting the count to 1 for all of the label generators.
+##
+## The motivation for resetLabelFunctionCreators is for testing: If we want to check
+## that two pathways to code generation (one existing, one experimental) create identical
+## code, it is helpful to have identical generated labels.  Resetting all label generators
+## supports this goal.
+labelFunctionMetaCreator <- function() {
+    allLabelFunctionCreators <- list()
+
+    creatorFun <- function(lead, start = 1) {
+        nextIndex <- start
+        force(lead)
+        labelGenerator <- function(reset = FALSE, count = 1, envName = "") {
+            if(reset) {
+                nextIndex <<- 1
+                return(invisible(NULL))
+            }
+            lead <- paste(lead, envName , sep = '_')
+            ans <- paste0(lead, nextIndex - 1 + (1:count))
+            nextIndex <<- nextIndex + count
+            ans
+        }
+        allLabelFunctionCreators[[ length(allLabelFunctionCreators) + 1 ]] <<- labelGenerator
+        labelGenerator
+    }
+    creatorFun
+}
+
+labelFunctionCreator <- labelFunctionMetaCreator()
+
+resetLabelFunctionCreators <- function() {
+    allLabelFunctionCreators <- environment(labelFunctionCreator)$allLabelFunctionCreators
+    for(i in allLabelFunctionCreators) {
+        i(reset = TRUE)
     }
-    lead <- paste(lead, envName , sep = '_')
-    ans <- paste0(lead, nextIndex - 1 + (1:count))
-    nextIndex <<- nextIndex + count
-    ans
-  }
-  labelGenerator
 }
 
 nimbleUniqueID <- labelFunctionCreator("UID")
diff --git a/packages/nimble/R/cppDefs_RCfunction.R b/packages/nimble/R/cppDefs_RCfunction.R
index 8852a2045..bb0abf232 100644
--- a/packages/nimble/R/cppDefs_RCfunction.R
+++ b/packages/nimble/R/cppDefs_RCfunction.R
@@ -96,7 +96,6 @@ RCfunctionDef <- setRefClass('RCfunctionDef',
                                    
                                    if(is.character(SEXPinterfaceCname) && is.null(dll) && eval) {
                                      warning("creating a .Call() expression with no DLL information")
-                                     browser()                                     
                                    }
                                                                      
                                    # avoid R CMD check problem with registration
diff --git a/packages/nimble/R/genCpp_exprClass.R b/packages/nimble/R/genCpp_exprClass.R
index c9ce337c4..14990692a 100644
--- a/packages/nimble/R/genCpp_exprClass.R
+++ b/packages/nimble/R/genCpp_exprClass.R
@@ -20,12 +20,13 @@ exprClass <- setRefClass('exprClass',
                              sizeExprs =  'ANY',	#'list', ## a list of size expressions (using R parse trees for each non-numeric expression
                              type =  'ANY',		#'character', ## type label
                              args =  'ANY',		#'list', ## list of exprClass objects for the arguments
-                             eigMatrix =  'ANY',		#'logical', ## vs. Array.  Only used for Eigenized expressions
-                             toEigenize =  'ANY',		#'character', ##'yes', 'no', or 'maybe'
-                             caller = 'ANY', ## exprClass object for the call to which this is an argument (if any)
-                             callerArgID =  'ANY',		#'numeric', ## index in the calling object's args list for this object.
-                             assertions =  'ANY',		#'list'
-                             cppADCode = 'ANY'    #'logical ## is expr in code generated for cppad?
+                             eigMatrix =  'ANY',	#'logical', ## vs. Array.  Only used for Eigenized expressions
+                             toEigenize =  'ANY',	#'character', ##'yes', 'no', or 'maybe'
+                             caller = 'ANY',            # exprClass object for the call to which this is an argument (if any)
+                             callerArgID =  'ANY',	#'numeric', ## index in the calling object's args list for this object.
+                             assertions =  'ANY',	#'list'
+                             cppADCode = 'ANY',         #'logical' ## is expr in code generated for cppad?
+                             aux = 'ANY'                # anything needed for specific operators
                              ),
                          methods = list(
                              initialize = function(...) {sizeExprs <<- list(); args <<- list();toEigenize <<- 'unknown';assertions <<- list(); eigMatrix <<- logical(); cppADCode <<- FALSE; callSuper(...)},
diff --git a/packages/nimble/R/genCpp_insertAssertions.R b/packages/nimble/R/genCpp_insertAssertions.R
index f2d691467..db5d8b2d9 100644
--- a/packages/nimble/R/genCpp_insertAssertions.R
+++ b/packages/nimble/R/genCpp_insertAssertions.R
@@ -34,7 +34,7 @@ exprClasses_insertAssertions <- function(code) {
             if(length(code$args[[i]]$assertions) > 0) {
                 toInsert <- lapply(code$args[[i]]$assertions, function(x) if(inherits(x, 'exprClass')) x else RparseTree2ExprClasses(x))
                 before <- unlist(lapply(toInsert, function(x) {if(x$name == 'after') FALSE else TRUE}))
-                
+                if(nimbleOptions('experimentalNewSizeProcessing')) code$args[[i]]$assertions <- list()  ## Clear assertions field so it can be used by later compiler stages
                 newExpr <- newBracketExpr(args = c(lapply(toInsert[before], ## assertions will be inserted recursively IF they are in {}
                                                           function(z) {exprClasses_insertAssertions(z); z}),
                                                    code$args[i],
diff --git a/packages/nimble/R/genCpp_sizeProcessing.R b/packages/nimble/R/genCpp_sizeProcessing.R
index 2215ba1f7..c37726a87 100644
--- a/packages/nimble/R/genCpp_sizeProcessing.R
+++ b/packages/nimble/R/genCpp_sizeProcessing.R
@@ -252,6 +252,8 @@ sizeDim <- function(code, symTab, typeEnv) {
 }
 
 sizeDiagonal <- function(code, symTab, typeEnv) {
+    ## experimentalNewSizeProcessing: code$name change step stays here
+    ## experimentalNewSizeProcessing: because the 3 cases are not implementation-specific
     asserts <- recurseSetSizes(code, symTab, typeEnv)
     argIsExprClass <- inherits(code$args[[1]], 'exprClass')
     nDimArg <- if(argIsExprClass) code$args[[1]]$nDim else 0
@@ -302,8 +304,10 @@ sizeWhich <- function(code, symTab, typeEnv) {
     code$toEigenize <- 'yes'
     code$name <- 'setWhich'
 
-    if(!(code$caller$name %in% assignmentOperators)) {
-        asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+    if(!nimbleOptions('experimentalSelfLiftStage')) {
+        if(!(code$caller$name %in% assignmentOperators)) {
+            asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+        }
     }
     if(length(asserts) == 0) NULL else asserts
 }
@@ -332,7 +336,7 @@ sizeRecyclingRule <- function(code, symTab, typeEnv) { ## also need an entry in
     code$sizeExprs <- newSizeExprs
     code$type <- 'double' ## will need to look up from a list
     code$nDim <- 1
-    code$toEigenize <- TRUE
+    code$toEigenize <- 'yes' ## toEigen: N.B. This had TRUE
     return(asserts)
 }
 
@@ -346,6 +350,7 @@ sizeRecyclingRuleRfunction <- function(code, symTab, typeEnv) {
     ## If scalar, that gives size
     ## If vector, size is length of first argument.
     ## Problem is vector of length 1, where size should be value of first element, not length of 1.
+    ## toEigen: keep this lift here for now, since it sets up sizes.
     if(inherits(code$args[[1]], 'exprClass')) {
         if(!code$args[[1]]$isName) {
             asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
@@ -895,7 +900,6 @@ sizeGetBound <- function(code, symTab, typeEnv) {
         asserts <- list()
     }
  
-    
     boundInfoSym <- symTab$getSymbolObject(code$args[[3]]$name, inherits = TRUE)
     code$type <- boundInfoSym$boundInfo$type
     code$nDim <- boundInfoSym$boundInfo$nDim
@@ -951,6 +955,7 @@ sizeAsRowOrCol <- function(code, symTab, typeEnv) {
 
 ## a$b becomes nfVar(a, 'b')
 sizeNFvar <- function(code, symTab, typeEnv) {
+    ## toEigen: Is it correct that this does not mark toEigen?
     asserts <- list()
     if(!inherits(code$args[[1]], 'exprClass'))
         stop(exprClassProcessingErrorMsg(code, 'Problem using $: no name on the right?'), call. = FALSE)
@@ -1068,10 +1073,12 @@ sizeNimbleListReturningFunction <- function(code, symTab, typeEnv) {
   code$sizeExprs <- symbolObject
   code$toEigenize <- "yes"  # This is specialized for nimSvd and nimEigen.
   if(code$name == 'getDerivs')
-    code$toEigenize <- 'no'  ## Temp. solution to ensure that derivsOrders argument is a nimArray and not an eigen type.
+      code$toEigenize <- 'no'  ## Temp. solution to ensure that derivsOrders argument is a nimArray and not an eigen type.
   code$nDim <- 0
-  if(!(code$caller$name %in% assignmentOperators))
-      asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+  if(!nimbleOptions('experimentalSelfLiftStage')) {
+      if(!(code$caller$name %in% assignmentOperators))
+          asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+  }
   if(length(asserts) == 0) NULL else asserts
 }
 
@@ -1357,6 +1364,7 @@ assignmentTypeWarn <- function(LHS, RHS) {
 }
 
 ## used for setAll
+## toEigen: N.B. This may be deprecated.
 sizeOneEigenCommand <- function(code, symTab, typeEnv) {
     if(!code$args[[1]]$isName) stop(exprClassProcessingErrorMsg(code, 'In sizeOneEigenCommand:  First arg should be a name.'), call. = FALSE)
     recurseSetSizes(code, symTab, typeEnv)
@@ -1385,6 +1393,7 @@ sizeforceEigenize <- function(code, symTab, typeEnv) {
     if(length(asserts) == 0) NULL else asserts
 }
 
+## toEigen: N.B. Deprecated (see external-c-calls in development)
 sizecallC <- function(code, symTab, typeEnv) {
     asserts <- recurseSetSizes(code$args[[1]], symTab, typeEnv)
     asserts
@@ -1479,6 +1488,7 @@ sizeSetSize <- function(code, symTab, typeEnv) {
 
 
 ## This was redundant and we should eventually be able to remove it
+## toEigen: N.B. omitting this
 sizeResizeNoPtr <- function(code, symTab, typeEnv){
     sym <- symTab$getSymbolObject(code$args[[1]]$name, inherits = TRUE)
     if(length(code$args[[2]]) != 1)  stop(exprClassProcessingErrorMsg(code, 'In sizeResizeNoPtr: Problem with number of dimensions provided in resize.'), call. = FALSE)
@@ -1525,7 +1535,9 @@ sizeInsertIntermediate <- function(code, argID, symTab, typeEnv, forceAssign = F
         newArgExpr <- RparseTree2ExprClasses(as.name(newName))
         newArgExpr$type <- code$args[[argID]]$type
         newArgExpr$sizeExprs <- code$args[[argID]]$sizeExprs
-        newArgExpr$toEigenize <- 'maybe'
+        if(!nimbleOptions('experimentalNewSizeProcessing')) {
+            newArgExpr$toEigenize <- 'maybe'
+        }
         newArgExpr$nDim <- code$args[[argID]]$nDim
     } else {
 
@@ -1542,7 +1554,9 @@ sizeInsertIntermediate <- function(code, argID, symTab, typeEnv, forceAssign = F
         newArgExpr <- RparseTree2ExprClasses(as.name(newName))
         newArgExpr$type <- newExpr$args[[1]]$type
         newArgExpr$sizeExprs <- newExpr$args[[1]]$sizeExprs
-        newArgExpr$toEigenize <- 'maybe'
+        if(!nimbleOptions('experimentalNewSizeProcessing')) {
+            newArgExpr$toEigenize <- 'maybe'
+        }
         newArgExpr$nDim <- newExpr$args[[1]]$nDim
     }
     setArg(code, argID, newArgExpr)
@@ -1564,8 +1578,6 @@ sizeAssign <- function(code, symTab, typeEnv) {
     if(length(asserts) == 0) NULL else asserts
 }
 
-
-
 ## Handler for assignment
 sizeAssignAfterRecursing <- function(code, symTab, typeEnv, NoEigenizeMap = FALSE) {
     LHS <- code$args[[1]]
@@ -1628,7 +1640,9 @@ sizeAssignAfterRecursing <- function(code, symTab, typeEnv, NoEigenizeMap = FALS
                 }
             } else { ## yes in symTab
                 ## This case is ok.  It is in the symbol table but not the typeEnv.  So it is something like ptr <- getPtr(A)
+                 if(!nimbleOptions('experimentalNewSizeProcessing')) {
                 code$toEigenize <- 'no'
+                 } ##experimentalNewSizeProcessing
                 code$nDim <- 0
                 code$type <- 'unknown'
                 code$sizeExprs <- list()
@@ -1666,6 +1680,8 @@ sizeAssignAfterRecursing <- function(code, symTab, typeEnv, NoEigenizeMap = FALS
     ## Note this can use LHS$name for RHSsizeExprs when returning from a nimbleFunction on RHS.  But this is probably not needed any more.
     if(any(unlist(lapply(RHSsizeExprs, is.null)))) RHSsizeExprs <- makeSizeExpressions(rep(NA, RHSnDim), LHS$name) ## reset sizeExprs for the LHS var. re-using RHSsizeExprs for LHS.  This would only be valid if it is a nimbleFunction returning something on the RHS.  For assignment to be executed in Eigen, the RHS sizes MUST be known
 
+     if(!nimbleOptions('experimentalNewSizeProcessing')) {
+                
     if(LHS$toEigenize == 'yes') {
         code$toEigenize <- 'yes'
 ##        message('Warning from sizeAssign: not expecting LHS to have toEigenize == yes')
@@ -1684,7 +1700,8 @@ sizeAssignAfterRecursing <- function(code, symTab, typeEnv, NoEigenizeMap = FALS
             else 'no'
         }
     }
-
+    
+                
     if(code$toEigenize == 'yes') { ## this would make more sense in eigenize_assign
     ## generate setSize(LHS, ...) where ... are dimension expressions
         if(length(RHSnDim) == 0) {
@@ -1737,7 +1754,7 @@ sizeAssignAfterRecursing <- function(code, symTab, typeEnv, NoEigenizeMap = FALS
             if(LHS$name == 'map') assert <- c(assert, sizeInsertIntermediate(code, 1, symTab, typeEnv) )
         }
     }
-
+     } ##experimentalNewSizeProcessing
     if(!(LHS$name %in% c('eigenBlock', 'diagonal', 'coeffSetter', 'nimNonseqIndexedd', 'nimNonseqIndexedi','nimNonseqIndexedb'))) {
         ## should already be annotated if it is an indexed assignment.
         ## It should be harmless to re-annotated EXCEPT in case like out[1:5] <- scalar
@@ -2392,9 +2409,11 @@ sizeUnaryCwise <- function(code, symTab, typeEnv) {
     ## lift intermediates
     a1 <- code$args[[1]]
     if(inherits(a1, 'exprClass')) {
-        if(a1$toEigenize == 'no') {
-            asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
-            a1 <- code$args[[1]]
+        if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+            if(a1$toEigenize == 'no') {
+                asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
+                a1 <- code$args[[1]]
+            }
         }
         code$nDim <- a1$nDim
         code$sizeExprs <- a1$sizeExprs
@@ -2404,7 +2423,7 @@ sizeUnaryCwise <- function(code, symTab, typeEnv) {
     }
     code$type <- setReturnType(code$name, getArgumentType(a1))
     if(length(code$nDim) != 1) stop(exprClassProcessingErrorMsg(code, 'In sizeUnaryCwise: nDim is not set.'), call. = FALSE)
-    code$toEigenize <- if(code$nDim > 0) 'yes' else 'maybe'
+    if(!nimbleOptions('experimentalNewSizeProcessing') ) code$toEigenize <- if(code$nDim > 0) 'yes' else 'maybe'
     return(asserts)
 }
 
@@ -2521,9 +2540,11 @@ sizeUnaryReduction <- function(code, symTab, typeEnv) {
                 stop(exprClassProcessingErrorMsg(code, 'NIMBLE compiler does not support var with a matrix (or higher dimensional) argument.'), call. = FALSE) 
             }
         }
-        if(!code$args[[1]]$isName) {
-            if(code$args[[1]]$toEigenize == 'no') {
-                asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
+        if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+            if(!code$args[[1]]$isName) {
+                if(code$args[[1]]$toEigenize == 'no') {
+                    asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
+                }
             }
         }
     }
@@ -2531,10 +2552,12 @@ sizeUnaryReduction <- function(code, symTab, typeEnv) {
     code$nDim <- 0
     code$sizeExprs <- list()
     code$type <- setReturnType(code$name, code$args[[1]]$type)
-    code$toEigenize <- 'yes'
+    if(!nimbleOptions('experimentalNewSizeProcessing') ) code$toEigenize <- 'yes'
 
-    if(!(code$caller$name %in% c('{','<-','<<-','='))) {
-        asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+    if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+        if(!(code$caller$name %in% c('{','<-','<<-','='))) {
+            asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+        }
     }
 
     if(length(asserts) == 0) NULL else asserts
@@ -2603,14 +2626,15 @@ sizeMatrixMult <- function(code, symTab, typeEnv) {
     a2 <- code$args[[2]]
     
     if(a1$nDim == 0 | a2$nDim == 0) stop(exprClassProcessingErrorMsg(code, 'In sizeMatrixMult: Cannot do matrix multiplication with a scalar.'), call. = FALSE) 
-
-    if(a1$toEigenize == 'no') {
-        asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
-        a1 <- code$args[[1]]
-    }
-    if(a2$toEigenize == 'no') {
-        asserts <- c(asserts, sizeInsertIntermediate(code, 2, symTab, typeEnv))
-        a2 <- code$args[[2]]
+    if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+        if(a1$toEigenize == 'no') {
+            asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
+            a1 <- code$args[[1]]
+        }
+        if(a2$toEigenize == 'no') {
+            asserts <- c(asserts, sizeInsertIntermediate(code, 2, symTab, typeEnv))
+            a2 <- code$args[[2]]
+        }
     }
     
     ## Note that we could insert RUN-TIME adaptation of mat %*% vec and vec %*% mat
@@ -2655,7 +2679,7 @@ sizeMatrixMult <- function(code, symTab, typeEnv) {
     code$nDim <- 2
     code$sizeExprs <- list(a1$sizeExprs[[1]], a2$sizeExprs[[2]])
     code$type <- setReturnType(code$name, arithmeticOutputType(a1$type, a2$type))
-    code$toEigenize <- 'yes'
+    if(!nimbleOptions('experimentalNewSizeProcessing') ) code$toEigenize <- 'yes'
     assertMessage <- paste0("Run-time size error: expected ", deparse(a1$sizeExprs[[2]]), " == ", deparse(a2$sizeExprs[[1]]))
     newAssert <- identityAssert(a1$sizeExprs[[2]], a2$sizeExprs[[1]], assertMessage)
     if(is.null(newAssert))
@@ -2766,40 +2790,44 @@ sizeBinaryCwise <- function(code, symTab, typeEnv) {
     a2 <- code$args[[2]]
     ## pull out aXDropNdim, aXnDim, aXsizeExprs, and aXtype (X = 1 or 2)
     if(inherits(a1, 'exprClass')) {
-        if(a1$toEigenize == 'no') {
-            asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
-            a1 <- code$args[[1]]
+        if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+            if(a1$toEigenize == 'no') {
+                asserts <- c(asserts, sizeInsertIntermediate(code, 1, symTab, typeEnv))
+                a1 <- code$args[[1]]
+            }
         }
         a1Drop <- dropSingleSizes(a1$sizeExprs)
         a1DropNdim <- length(a1Drop$sizeExprs)
         a1nDim <- a1$nDim
         a1sizeExprs <- a1$sizeExprs
         a1type <- a1$type
-        a1toEigenize <- a1$toEigenize
+        if(!nimbleOptions('experimentalNewSizeProcessing') ) a1toEigenize <- a1$toEigenize
     } else {
         a1DropNdim <- 0
         a1nDim <- 0
         a1sizeExprs <- list()
         a1type <- storage.mode(a1)
-        a1toEigenize <- 'maybe'
+        if(!nimbleOptions('experimentalNewSizeProcessing') ) a1toEigenize <- 'maybe'
     }
     if(inherits(a2, 'exprClass')) {
-        if(a2$toEigenize == 'no') {
-            asserts <- c(asserts, sizeInsertIntermediate(code, 2, symTab, typeEnv))
-            a2 <- code$args[[2]]
+        if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+            if(a2$toEigenize == 'no') {
+                asserts <- c(asserts, sizeInsertIntermediate(code, 2, symTab, typeEnv))
+                a2 <- code$args[[2]]
+            }
         }
         a2Drop <- dropSingleSizes(a2$sizeExprs)
         a2DropNdim <- length(a2Drop$sizeExprs)
         a2nDim <- a2$nDim
         a2sizeExprs <- a2$sizeExprs
         a2type <- a2$type
-        a2toEigenize <- a2$toEigenize
+        if(!nimbleOptions('experimentalNewSizeProcessing') )  a2toEigenize <- a2$toEigenize
     } else {
         a2DropNdim <- 0
         a2nDim <- 0
         a2sizeExprs <- list()
         a2type <- storage.mode(a2)
-        a2toEigenize <- 'maybe'
+        if(!nimbleOptions('experimentalNewSizeProcessing') )  a2toEigenize <- 'maybe'
     }
     
     ## Choose the output type by type promotion
@@ -2807,15 +2835,18 @@ sizeBinaryCwise <- function(code, symTab, typeEnv) {
     if(length(a2type) == 0) {warning('Problem with type of arg2 in sizeBinaryCwise', call. = FALSE); browser()}
     code$type <- setReturnType(code$name, arithmeticOutputType(a1type, a2type))
 
-    forceYesEigenize <- identical(a1toEigenize, 'yes') | identical(a2toEigenize, 'yes')
-    code$toEigenize <- if(a1DropNdim == 0 & a2DropNdim == 0)
-                           if(forceYesEigenize)
-                               'yes'
-                           else
-                               'maybe'
-                       else 'yes'
-     
+    if(!nimbleOptions('experimentalNewSizeProcessing') ) {
+        forceYesEigenize <- identical(a1toEigenize, 'yes') | identical(a2toEigenize, 'yes')
+        code$toEigenize <- if(a1DropNdim == 0 & a2DropNdim == 0)
+                               if(forceYesEigenize)
+                                   'yes'
+                               else
+                                   'maybe'
+                           else 'yes'
+    }
+    
     ## Catch the case that there is at least one scalar-equivalent (all lengths == 1)
+    ## experimentalNewSizeProcessing: The 3 'code$toEigenize <- ' should be redundant with above and could be removed during refactor
     if(a1DropNdim == 0 | a2DropNdim == 0) { 
         ## Here we will process effective scalar additions
         ## and not do any other type of size promotion/dropping
@@ -2823,17 +2854,17 @@ sizeBinaryCwise <- function(code, symTab, typeEnv) {
             if(a2DropNdim == 0) { ##both are scalar-equiv
                 code$nDim <- max(a1nDim, a2nDim) ## use the larger nDims
                 code$sizeExprs <- rep(list(1), code$nDim) ## set sizeExprs to all 1
-                code$toEigenize <- if(forceYesEigenize) 'yes' else 'maybe'
+                if(!nimbleOptions('experimentalNewSizeProcessing') ) code$toEigenize <- if(forceYesEigenize) 'yes' else 'maybe'
             } else {
                 ## a2 is not scalar equiv, so take nDim and sizeExprs from it
                 code$nDim <- a2nDim
                 code$sizeExprs <- a2sizeExprs
-                code$toEigenize <- 'yes'
+                if(!nimbleOptions('experimentalNewSizeProcessing') ) code$toEigenize <- 'yes'
             }
         } else { ## a2 is scalar-equiv, and a1 is not
             code$nDim <- a1nDim
             code$sizeExprs <- a1sizeExprs
-            code$toEigenize <- 'yes'
+            if(!nimbleOptions('experimentalNewSizeProcessing') ) code$toEigenize <- 'yes'
         }
         return(if(length(asserts) == 0) NULL else asserts)
     }
diff --git a/packages/nimble/R/genCpp_toEigenize.R b/packages/nimble/R/genCpp_toEigenize.R
new file mode 100644
index 000000000..344a505d4
--- /dev/null
+++ b/packages/nimble/R/genCpp_toEigenize.R
@@ -0,0 +1,472 @@
+## It is hoped substantial sets of these calls can be combined
+## or implemented by a rule system.
+
+toEigenizeNoCalls <- c('dim', 'run.time','nimOptimDefaultControl')
+
+toEigenizeYesCalls <- c(paste0('nimDiagonal', c('D','I','B')),
+                        'diagonal',
+                        'which',
+                        scalar_distribution_dFuns,
+                        scalar_distribution_pFuns,
+                        scalar_distribution_qFuns,
+                        paste0(c('d','q','p'), 't'),
+                        paste0(c('d','q','p'), 'exp'),
+                        scalar_distribution_rFuns,
+                        'rt',
+                        'rexp',
+                        'asRow','asCol',
+                        nimbleListReturningOperators, ## should these be in 'maybe'?
+                        'setAll' ## deprecate?
+                        )
+
+toEigenizeMaybeCalls <- c('map',
+                          c('decide', 'size', 'getsize','getNodeFunctionIndexedInfo', 'endNimbleTimer'),
+                          c('blank', 'nfMethod', 'getPtr', 'startNimbleTimer'))
+
+toEigenizeUseRuleCalls <- c('nimPrint')
+
+toEigenCalls <- c(makeCallList(binaryOperators, 'toEigenBinaryCwise'),             #ok
+                  makeCallList(binaryMidLogicalOperators, 'sizeBinaryCwiseLogical'),
+                  makeCallList(binaryOrUnaryOperators, 'toEigenBinaryUnaryCwise'), #ok
+                  makeCallList(unaryOperators, 'toEigenUnaryCwise'),               #ok
+                  makeCallList(unaryOrNonaryOperators, 'sizeUnaryNonaryCwise'),
+                  makeCallList(assignmentOperators, 'toEigenAssign'),              #ok
+                  makeCallList(reductionUnaryOperators, 'sizeUnaryReduction'),     # drafted
+                  makeCallList(matrixSquareReductionOperators, 'sizeMatrixSquareReduction'),
+                  makeCallList(reductionBinaryOperators, 'sizeBinaryReduction'),
+                  makeCallList(matrixMultOperators, 'toEigenMatrixMult'),          # drafted
+                  makeCallList(matrixFlipOperators, 'sizeTranspose'),
+                  makeCallList(matrixSolveOperators, 'sizeSolveOp'), 
+                  makeCallList(matrixSquareOperators, 'sizeUnaryCwiseSquare'),
+                  ##               makeCallList(nimbleListReturningOperators, 'sizeNimbleListReturningFunction'),
+                  nimOptim = 'toEigenOptim',
+                  ##nimOptimDefaultControl = 'sizeOptimDefaultControl',
+                  
+                  makeCallList(paste0('nimC',c('d','i','b')), 'toEigenConcatenate'),
+                  makeCallList(paste0('nimRep',c('d','i','b')), 'toEigenRep'),
+                  list('debugSizeProcessing' = 'sizeProxyForDebugging',
+                       ##      diag = 'toEigenYes',
+                       ##      dim = 'toEigenNo',
+                       ##       RRtest_add = 'toEigenYes', ## testing!
+                       ##      which = 'toEigenYes',
+                       
+                       nimRep = 'sizeRep',
+                       nimSeqBy = 'sizeSeq',
+                       nimSeqLen = 'sizeSeq',
+                       nimSeqByLen = 'sizeSeq',
+                       'return' = 'sizeReturn',
+                       ##  'asRow' = 'sizeAsRowOrCol',
+                       ##  'asCol' = 'sizeAsRowOrCol',
+                       makeNewNimbleListObject = 'toEigenNewNimbleList',
+                       getParam = 'toEigenGetParam',
+                       getBound = 'toEigenGetBound',
+                       ## nimSwitch = 'sizeSwitch', ## NEVER GETS MARKED WITH TOEIGENIZE
+                       asDoublePtr = 'toEigenIgnoreForNow',
+                       '[' = 'sizeIndexingBracket',
+                       ##'[[' = 'sizeDoubleBracket', ## for nimbleFunctionList, this will always  go through chainedCall(nfList[[i]], 'foo')(arg1, arg2) ## NEVER GETS MARKED WITH TOEIGENIZE
+                       chainedCall = 'toEigenChainedCall',
+                       ## nfVar = 'sizeNFvar', ## NEVER GETS MARKED WITH TOEIGENIZE
+                       ## map = 'sizemap',
+                       ':' = 'sizeColonOperator',
+                       ##dim = 'sizeDimOperator',
+                       'if' = 'recurseSetSizes', ##OK
+                       'while' = 'recurseSetSizes',
+                       callC = 'sizecallC', 
+                       'for' = 'toEigenFor', 
+                       cppPointerDereference = 'toEigenCppPointerDereference',
+                       values = 'toEigenValues',
+                       '(' = 'sizeUnaryCwise',
+                       setSize = 'toEigenSetSize', ## OK but not done for numericLists
+                       ## resizeNoPtr = 'sizeResizeNoPtr', ## may not be used any more 
+                       nimArr_rcat = 'toEigenScalarRecurse',
+                       nimArr_rinterval = 'toEigenScalarRecurse',
+                       ##nimPrint = 'sizeforceEigenize',
+                       ##nimCat = 'sizeforceEigenize',
+                       as.integer = 'sizeUnaryCwise', ## Note as.integer and as.numeric will not work on a non-scalar yet
+                       as.numeric = 'sizeUnaryCwise',
+                       nimArrayGeneral = 'toEigenNimArrayGeneral',
+                       ##setAll = 'sizeOneEigenCommand',
+                       voidPtr = 'sizeVoidPtr'
+                       ),
+                    ##run.time = 'sizeRunTime'),
+                  ##               makeCallList(scalar_distribution_dFuns, 'toEigenYes'),
+                  ##               makeCallList(scalar_distribution_pFuns, 'toEigenYes'),
+                  ##               makeCallList(scalar_distribution_qFuns, 'sizeRecyclingRule'),
+                  ##               makeCallList(scalar_distribution_rFuns, 'sizeRecyclingRuleRfunction'),
+                  makeCallList(distributionFuns[!(distributionFuns %in% c(scalar_distribution_dFuns, scalar_distribution_pFuns, scalar_distribution_qFuns, scalar_distribution_rFuns))], 'toEigenScalarRecurse'),
+                                        # R dist functions that are not used by NIMBLE but we allow in DSL
+                  ## makeCallList(paste0(c('d','q','p'), 't'), 'sizeRecyclingRule'),
+                  ## rt = 'sizeRecyclingRuleRfunction',
+                  ## makeCallList(paste0(c('d','q','p'), 'exp'), 'sizeRecyclingRule'),
+                  ## rexp = 'sizeRecyclingRuleRfunction',
+                  makeCallList(c('isnan','ISNAN','ISNA'), 'toEigenScalarRecurse'),
+                  makeCallList(c('nimArr_dmnorm_chol', 'nimArr_dmvt_chol', 'nimArr_dwish_chol', 'nimArr_dinvwish_chol', 'nimArr_dmulti', 'nimArr_dcat', 'nimArr_dinterval', 'nimArr_ddirch'), 'toEigenScalarRecurse'),
+                  makeCallList(c('nimArr_rmnorm_chol', 'nimArr_rmvt_chol', 'nimArr_rwish_chol', 'nimArr_rinvwish_chol', 'nimArr_rmulti', 'nimArr_rdirch'), 'sizeRmultivarFirstArg'),
+                  ##      makeCallList(c('decide', 'size', 'getsize','getNodeFunctionIndexedInfo', 'endNimbleTimer'), 'sizeScalar'),
+                  ## makeCallList(c('calculate','calculateDiff', 'getLogProb'), 'toEigenScalarModelOp'),
+                  simulate = 'toEigenIgnoreForNow'
+                  ##               makeCallList(c('blank', 'nfMethod', 'getPtr', 'startNimbleTimer'), 'sizeUndefined')
+                  )
+
+
+exprClasses_setToEigenize <- function(code, symTab, typeEnv) { ## input code is exprClass
+    ## name:
+    if(code$isName) {
+        return(list()) ## names don't have toEigen annotation
+    }
+    if(code$isCall) {
+        if(code$name == '{') {
+            ## recurse over lines
+            for(i in seq_along(code$args)) {
+                if(inherits(code$args[[i]], 'exprClass')) {
+                    newAsserts <- exprClasses_setToEigenize(code$args[[i]], symTab, typeEnv)
+                    code$args[[i]]$assertions <- c(code$args[[i]]$assertions,
+                                                   if(is.null(newAsserts)) list() else newAsserts)
+                }
+            }
+            return(invisible(NULL))
+        }
+        thisCall <- toEigenCalls[[code$name]]
+        if(!is.null(thisCall)) {
+            test0 <- eval(call(thisCall, code, symTab, typeEnv))
+            return(test0)
+        }
+        if(symTab$symbolExists(code$name, TRUE)) { ## could be a nimbleFunction object
+            return(toEigenNimbleFunction(code, symTab, typeEnv) )
+        }
+        ## Finally, it could be an RCfunction (a nimbleFunction with no setup == a simple function) {
+
+        if(exists(code$name)) {
+            obj <- get(code$name)
+            if(is.rcf(obj)) { ## it is an RC function
+                message('figure out how to extract the RCfunProc that should already exist')
+                browser()
+                RCfunProc <- typeEnv$.nimbleProject$compileRCfun(obj, initialTypeInference = TRUE)
+                
+                return(toEigenRCfunction(code, symTab, typeEnv, nfmObj, RCfunProc))
+            }
+        }
+    }
+    invisible(NULL)
+}
+
+
+recurseSetToEigenize <- function(code, symTab, typeEnv, useArgs = rep(TRUE, length(code$args))) {
+    asserts <- list()
+    for(i in seq_along(code$args)) {
+        if(useArgs[i]) {
+            if(inherits(code$args[[i]], 'exprClass')) {
+                asserts <- c(asserts, exprClasses_setToEigenize(code$args[[i]], symTab, typeEnv))
+            }
+        }
+    }
+    if(length(asserts)==0) NULL else asserts
+}
+
+toEigenInsertIntermediate <- function(code, argID, symTab, typeEnv) {
+    newLine <- sizeInsertIntermediate(code, argID, symTab, typeEnv)
+    exprClasses_setToEigenize(newLine, symTab, typeEnv)
+    newLine
+}
+
+toEigenNo <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'no'
+    return(list())
+}
+
+toEigenYes <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'yes'
+    return(list())
+}
+
+toEigenNo <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'no'
+    return(list())
+}
+
+toEigenMaybe <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'maybe'
+    return(list())
+}
+
+toEigenUseRule <- function(code, symTab, typeEnv) {
+    message('need to implement rules for toEigenUseRule')
+    return(list())
+}
+
+toEigenAssign <- function(code, symTab, typeEnv) {
+    asserts <- recurseSetToEigenize(code, symTab, typeEnv, useArgs = c(FALSE, TRUE))
+    asserts <- c(asserts, toEigenAssignLHS(code, symTab, typeEnv))
+    asserts
+}
+
+toEigenAssignLHS <- function(code, symTab, typeEnv, NoEigenizeMap = FALSE) {
+    ## In original size processing, NoEigenizeMap exists to be set as TRUE from sizeInsertIntermediate.
+    if(identical(code$type, 'unknown')) {
+        code$toEigenize <- 'no'
+        return(list())
+    }
+    assert <- list()
+    LHS <- code$args[[1]]
+    RHS <- code$args[[2]]
+    
+    if(LHS$toEigenize == 'yes') { ## unexpected?
+        code$toEigenize <- 'yes'
+    } else {
+        code$toEigenize <-
+            if(inherits(RHS, 'exprClass')) {
+                if(RHS$toEigenize == 'no') 'no'
+                else {
+                    if(RHS$toEigenize == 'unknown') 'no'
+                    else {
+                        if(RHS$toEigenize != 'yes' &
+                           (!(LHS$name %in% c('eigenBlock', 'diagonal', 'coeffSetter'))) &
+                           (RHS$nDim == 0 | RHS$isName | (RHS$name == 'map' & NoEigenizeMap)))
+                            'no' ## if it is scalar or is just a name or a map, we will do it via NimArr operator= .  Used to have "| RHS$name == 'map'", but this allowed X[1:3] <- X[2:4], which requires eigen, with eval triggered, to get right
+                        else 'yes' ## if it is 'maybe' and non-scalar and not just a name, default to 'yes'
+                    }
+                }
+            } else {
+                if(is.numeric(LHS$nDim))
+                    if(LHS$nDim > 0) 'yes' ## This is for cases like out[1:4] <- scalar
+                    else 'no'
+                else 'no'
+            }
+    }
+    ## if LHS is set to eigenize 'yes set code to Eigenize 'yes'
+    
+    if(code$toEigenize == 'yes') { ## this would make more sense in eigenize_assign
+        ## generate setSize(LHS, ...) where ... are dimension expressions
+        if(length(RHS$nDim) == 0) {
+            message("confused about trying to eigenize something with nDim = 0")
+            browser()
+        }
+        if(RHS$nDim > 0) {
+            if(!(RHS$name %in% setSizeNotNeededOperators)) {
+                if(LHS$isName | LHS$name == "nfVar") {
+                    assert <- substitute(setSize(LHS), list(LHS = nimbleGeneralParseDeparse(LHS)))
+                    for(i in seq_along(RHS$sizeExprs)) { ## N.B. RHS$sizeExprs may be modified by sizeAssignAfterRecursing
+                        test <- try(assert[[i + 2]] <- RHS$sizeExprs[[i]])
+                        if(inherits(test, 'try-error')) browser()
+                    }
+                    assert[[ length(assert) + 1]] <- 0 ## copyValues = false
+                    assert[[ length(assert) + 1]] <- 0 ## fillZeros  = false
+                    assert <- list(assert)
+                } else { ## We have an indexed LHS of an eigenizable expression
+                    ## need special handling if it is a row assignment like x[i,] <- ...
+                    ## also need to generate size assertions                    
+                    if(LHS$nDim == 1) {
+                        if(RHS$nDim == 2) {
+                            if(is.numeric(RHS$sizeExprs[[1]])) {
+                                if(RHS$sizeExprs[[1]] == 1) {
+                                    newExpr <- insertExprClassLayer(code, 1, 'asRow', type = LHS$type)
+                                    newExpr$sizeExprs <- RHS$sizeExprs 
+                                    newExpr$type <- LHS$type
+                                    newExpr$nDim <- RHS$nDim
+                                    if(!is.numeric(LHS$sizeExprs[[1]]) | !is.numeric(RHS$sizeExprs[[2]])) {
+                                        assertMessage <- paste0("Run-time size error: expected ", deparse(LHS$sizeExprs[[1]]), " == ", deparse(RHS$sizeExprs[[2]]))
+                                        thisAssert <- identityAssert(LHS$sizeExprs[[1]], RHS$sizeExprs[[2]], assertMessage)
+                                        if(!is.null(thisAssert)) assert[[length(assert) + 1]] <- thisAssert 
+                                    } else {
+                                        if(LHS$sizeExprs[[1]] != RHS$sizeExprs[[2]]) stop(exprClassProcessingErrorMsg(code, paste0('In sizeAssignAfterRecursing: Fixed size mismatch.')), call. = FALSE)
+                                    }
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    } else {
+        if(inherits(RHS, 'exprClass')) {
+            ## If we have A <- map(B, ...), we need to generate a setMap for the RHS, which will be done by sizeInsertIntermediate 
+            if(RHS$name == 'map') assert <- c(assert, toEigenInsertIntermediate(code, 2, symTab, typeEnv) )
+        }
+        if(inherits(LHS, 'exprClass')) {
+                                        # ditto
+            if(LHS$name == 'map') assert <- c(assert, toEigenInsertIntermediate(code, 1, symTab, typeEnv) )
+        }
+    }
+    assert
+    
+}
+        
+toEigenConcatenate <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'yes'
+    return(list())
+}
+
+toEigenRep <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'yes'
+    return(list())
+}
+
+toEigenNewNimbleList <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'maybe'
+    return(list())
+}
+
+toEigenNimArrayGeneral <- function(code, symTab, typeEnv) {
+    ## won't even have the same name at this point
+    stop('setting toEigenize for nimArrayGeneral must be handled')
+}
+
+toEigenGetParam <- function(code, symTab, typeEnv) {
+    stop('setting toEigenize for getParam must be handled')
+    return(list())
+}
+
+toEigenGetBound <- function(code, symTab, typeEnv) {
+    stop('setting toEigenize for getBound must be handled')
+    return(list())
+}
+
+toEigenOptim <- function(code, symTab, typeEnv) {
+    code$toEigenize <- 'no'
+    return(list())
+}
+
+toEigenCppPointerDereference <- function(code, symTab, typeEnv) {
+    code$toEigenize <- code$args[[1]]$toEigenize
+    return(list())
+}
+
+toEigenChainedCall <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for chained call must be handled')
+    return(list())
+}
+
+toEigenValues <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for values must be handled')
+    return(list())
+}
+
+toEigenRCfunction <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for RCfunction must be handled')
+    return(list())
+}
+
+toEigenNimbleFunction <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for RCfunction must be handled')
+    return(list())
+}
+
+toEigenSetSize <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for toEigenSetSize must be handled')
+    return(list())
+}
+
+toEigenFor <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for toEigenFor must be handled')
+    return(list())
+}
+
+toEigenIgnoreForNow <- function(code, symTab, typeEnv) {
+    message('Can not set toEigenize: to ignore for now')
+    return(list())
+}
+
+toEigenScalarModelOp <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for toEigenScalarModelOp must be handled')
+    return(list())
+}
+
+toEigenScalarRecurse <- function(code, symTab, typeEnv) {
+    message('setting toEigenize for toEigenScalarRecurse must be handled')
+    return(list())
+}
+
+toEigenBinaryUnaryCwise <- function(code, symTab, typeEnv) {
+    if(length(code$args) == 1) return(toEigenUnaryCwise(code, symTab, typeEnv))
+    if(length(code$args) == 2) return(toEigenBinaryCwise(code, symTab, typeEnv))
+    stop(exprClassProcessingErrorMsg(code, paste0('In toEigenBinaryUnarycWise: Length of arguments is not 1 or 2.')), call. = FALSE)
+}
+
+toEigenBinaryCwise <- function(code, symTab, typeEnv) {
+    asserts <- recurseSetToEigenize(code, symTab, typeEnv)
+    a1 <- code$args[[1]] 
+    a2 <- code$args[[2]]
+    if(inherits(a1, 'exprClass')) {
+        if(a1$toEigenize == 'no') {
+            asserts <- c(asserts, toEigenInsertIntermediate(code, 1, symTab, typeEnv))
+            a1 <- code$args[[1]]
+        }
+        a1DropNdim <- length(dropSingleSizes(a1$sizeExprs)$sizeExprs)
+        a1toEigenize <- a1$toEigenize
+    } else {
+        a1DropNdim <- 0
+        a1toEigenize <- 'maybe'
+    }
+    if(inherits(a2, 'exprClass')) {
+        if(a2$toEigenize == 'no') {
+            asserts <- c(asserts, toEigenInsertIntermediate(code, 2, symTab, typeEnv))
+            a2 <- code$args[[2]]
+        }
+        a2DropNdim <- length(dropSingleSizes(a2$sizeExprs)$sizeExprs)
+        a2toEigenize <- a2$toEigenize
+    } else {
+        a2DropNdim <- 0
+        a2toEigenize <- 'maybe'
+    }
+    forceYesEigenize <- identical(a1toEigenize, 'yes') | identical(a2toEigenize, 'yes')
+    code$toEigenize <- if(a1DropNdim == 0 & a2DropNdim == 0)
+                           if(forceYesEigenize)
+                               'yes'
+                           else
+                               'maybe'
+                       else 'yes'
+    asserts
+}
+
+toEigenUnaryCwise <- function(code, symTab, typeEnv) {
+    asserts <- recurseSetSizes(code, symTab, typeEnv)
+    a1 <- code$args[[1]]
+    ## lifting rule: lift if not eigenizable
+    if(inherits(a1, 'exprClass')) {
+        if(a1$toEigenize == 'no') {
+            asserts <- c(asserts, toEigenInsertIntermediate(code, 1, symTab, typeEnv))
+            a1 <- code$args[[1]]
+        }
+    }
+    ## propagation rule: yes or maybe
+    code$toEigenize <- if(code$nDim > 0) 'yes' else 'maybe'
+    asserts
+}
+
+toEigenMatrixMult <- function(code, symTab, typeEnv) {
+    asserts <- recurseSetToEigenize(code, symTab, typeEnv)
+    a1 <- code$args[[1]]
+    a2 <- code$args[[2]]
+    ## POSSIBLE DIFFERENT BEHAVIOR:
+    ## sizeMatrixMult may insert an asRow() or asCol()
+    ## It would do so *after* sizeInsertIntermediate
+    ## But now that would happen *before* sizeInsertIntermediate
+    ## Options: (1) The new code made be just as valid.  (2) We could make toEigenInsertIntermediate be smart about asRow or asCol
+    if(a1$toEigenize == 'no') {
+        asserts <- c(asserts, toEigenInsertIntermediate(code, 1, symTab, typeEnv))
+    }
+    if(a2$toEigenize == 'no') {
+        asserts <- c(asserts, toEigenInsertIntermediate(code, 2, symTab, typeEnv))
+    }
+    code$toEigenize <- 'yes'
+    asserts
+}
+
+toEigenUnaryReduction <- function(code, symTab, typeEnv) {
+    asserts <- recurseSetToEigenize(code, symTab, typeEnv)
+    if(inherits(code$args[[1]], 'exprClass')) {
+        if(!code$args[[1]]$isName) {
+            if(code$args[[1]]$toEigenize == 'no') {
+                asserts <- c(asserts, toEigenInsertIntermediate(code, 1, symTab, typeEnv))
+            }
+        }
+    }
+
+    code$toEigenize <- 'yes'
+
+    if(!(code$caller$name %in% c('{','<-','<<-','='))) {
+        asserts <- c(asserts, toEigenInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv))
+    }
+    
+    if(length(asserts) == 0) NULL else asserts
+}
diff --git a/packages/nimble/R/options.R b/packages/nimble/R/options.R
index 1892c64bf..e9f9dd244 100644
--- a/packages/nimble/R/options.R
+++ b/packages/nimble/R/options.R
@@ -7,6 +7,8 @@ nimbleUserNamespace <- as.environment(list(sessionSpecificDll = NULL))
 # These options are for development use at this point.
 .nimbleOptions <- as.environment(
     list(
+        experimentalNewSizeProcessing = FALSE,
+        experimentalSelfLiftStage = FALSE,
         enableSpecialHandling = FALSE,
         pauseAfterWritingFiles = FALSE,
         CppAD_directory = NA,
diff --git a/packages/nimble/inst/tests/mathTestLists.R b/packages/nimble/inst/tests/mathTestLists.R
new file mode 100644
index 000000000..51ca794ca
--- /dev/null
+++ b/packages/nimble/inst/tests/mathTestLists.R
@@ -0,0 +1,197 @@
+1### INSTRUCTIONS:
+## enter each test as a list, with an informative name, NIMBLE expression to evaluate, vector of input dimensions, value of output dimension, and (if NIMBLE expression cannot be directly evaluated in R) the equivalent pure R expression whose result should match the NIMBLE result
+
+testsVaried = list(
+  list(name = "matrix direct product", expr = quote(out <- arg1 * arg2), inputDim = c(2,2), outputDim = 2),
+  list(name = "matrix direct product with scalar addition", expr = quote(out <- (arg1+1) * (arg2+1)), inputDim = c(2,2), outputDim = 2),
+  list(name = "matrix absolute value", expr = quote(out <- abs(arg1)), inputDim = c(2), outputDim = 2),
+  list(name = "matrix absolute value with scalar addition", expr = quote(out <- abs(arg1 - 2)), inputDim = c(2), outputDim = 2),
+  list(name = "vector pmin", expr = quote(out <- pmin(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
+  list(name = "vector pmax", expr = quote(out <- pmax(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
+  list(name = "sd with addition", expr = quote(out <- sd(arg1) + 3), inputDim = c(1), outputDim = 0),
+  list(name = "sd of vector with addition", expr = quote(out <- sd(arg1 + 3)), inputDim = c(1), outputDim = 0),
+  list(name = "sd of matrix-vector multiply", expr = quote(out <- sd(arg1 %*% arg2)), inputDim = c(2,1), outputDim = 0),
+  list(name = "var of vector", expr = quote(out <- var(arg1)), inputDim = c(1), outputDim = 0),
+  list(name = "log determinant", expr = quote(out <- logdet(arg1)), inputDim = c(2), outputDim = 0)
+  )
+
+testsBasicMath = list(
+  list(name = 'exp of scalar', expr = quote(out <- exp(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'log of scalar', expr = quote(out <- log(abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'sqrt of scalar', expr = quote(out <- sqrt(abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'abs of scalar', expr = quote(out <- abs(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'step of scalar', expr = quote(out <- step(arg1)), inputDim = 0, outputDim = 0, Rcode = quote( out <- as.numeric(arg1 > 0))),
+  list(name = 'cube of scalar', expr = quote(out <- cube(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'cos of scalar', expr = quote(out <- cos(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'acos of cos of scalar', expr = quote(out <- acos(cos(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'sin of scalar', expr = quote(out <- sin(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'asin of sin of scalar', expr = quote(out <- asin(sin(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'tan of scalar', expr = quote(out <- tan(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'atan of tan of scalar', expr = quote(out <- atan(tan(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'cosh of scalar', expr = quote(out <- cosh(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'sinh of scalar', expr = quote(out <- sinh(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'tanh of scalar', expr = quote(out <- tanh(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'acosh of scalar', expr = quote(out <- acosh(1 + abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'asinh of scalar', expr = quote(out <- asinh(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'atanh of scalar', expr = quote(out <- atanh(abs(arg1)%%1)), inputDim = 0, outputDim = 0),
+  ###
+  list(name = 'exp of vector', expr = quote(out <- exp(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'log of vector', expr = quote(out <- log(abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'sqrt of vector', expr = quote(out <- sqrt(abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'abs of vector', expr = quote(out <- abs(arg1)), inputDim = 1, outputDim = 1),
+##  list(name = 'step of vector', expr = quote(out <- step(arg1)), inputDim = 1, outputDim = 1, Rcode = quote(out <- as.numeric(arg1 > 0))),   ## FAILS on compileNimble(nfR) with Eigen error
+  list(name = 'cube of vector', expr = quote(out <- cube(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'cos of vector', expr = quote(out <- cos(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'acos of cos of vector', expr = quote(out <- acos(cos(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'sin of vector', expr = quote(out <- sin(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'asin of sin of vector', expr = quote(out <- asin(sin(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'tan of vector', expr = quote(out <- tan(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'atan of tan of vector', expr = quote(out <- atan(tan(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'cosh of vector', expr = quote(out <- cosh(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'sinh of vector', expr = quote(out <- sinh(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'tanh of vector', expr = quote(out <- tanh(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'acosh of vector', expr = quote(out <- acosh(1 + abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'asinh of vector', expr = quote(out <- asinh(arg1)), inputDim = 1, outputDim = 1),
+##  list(name = 'atanh of vector', expr = quote(out <- atanh(arg1%%1)), inputDim = 1, outputDim = 1), ## FAILS - issue here is probably that modulo on vecs doesn't work but need to restrict domain for atanh
+  ###
+  list(name = 'scalar + scalar', expr = quote(out <- arg1 + arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'diff of scalars', expr = quote(out <- arg1 - arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'product of scalars', expr = quote(out <- arg1 * arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'ratio of scalars', expr = quote(out <- arg1 / arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'power of scalars via ^', expr = quote(out <- arg1 ^ arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'power of scalars via pow', expr = quote(out <- pow(arg1, arg2)), inputDim = c(0,0), outputDim = 0),
+  list(name = 'power of scalars via ^ with positive first arg', expr = quote(out <- exp(arg1) ^ arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'power of scalars via pow with positive first arg', expr = quote(out <- pow(exp(arg1), arg2)), inputDim = c(0,0), outputDim = 0),
+  list(name = 'modulo of scalars', expr = quote(out <- arg1 %% arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'min of scalars', expr = quote(out <- min(arg1, arg2)), inputDim = c(0,0), outputDim = 0),
+  list(name = 'max of scalars', expr = quote(out <- max(arg1, arg2)), inputDim = c(0,0), outputDim = 0),
+  ###
+  list(name = 'vector + vector', expr = quote(out <- arg1 + arg2), inputDim = c(1,1), outputDim = 1),
+  list(name = 'diff of vectors', expr = quote(out <- arg1 - arg2), inputDim = c(1,1), outputDim = 1),
+  list(name = 'product of vectors', expr = quote(out <- arg1 * arg2), inputDim = c(1,1), outputDim = 1),
+  list(name = 'ratio of vectors', expr = quote(out <- arg1 / arg2), inputDim = c(1,1), outputDim = 1),
+##  list(name = 'power of vectors via ^', expr = quote(out <- arg1 ^ arg2), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen casting
+  ## list(name = 'power of vectors via pow', expr = quote(out <- pow(arg1, arg2)), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen casting
+  ## list(name = 'modulo of vectors', expr = quote(out <- arg1 %% arg2), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen casting
+  list(name = 'pmin of vectors', expr = quote(out <- pmin(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
+  list(name = 'pmax of vectors', expr = quote(out <- pmax(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
+  ###
+  list(name = 'vector + scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'diff of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'product of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'ratio of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'power of vector and scalar via ^', expr = quote(out <- arg1 ^ arg2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'power of vector and scalar via pow', expr = quote(out <- pow(arg1, arg2)), inputDim = c(1,0), outputDim = 1),
+  list(name = 'power of vector and constant via ^', expr = quote(out <- arg1 ^ 2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'power of vector and constant via pow', expr = quote(out <- pow(arg1, 2)), inputDim = c(1,0), outputDim = 1),
+  list(name = 'power of vector and scalar via ^ with positive first arg', expr = quote(out <- exp(arg1) ^ arg2), inputDim = c(1,0), outputDim = 1),
+  list(name = 'power of vector and scalar via pow with positive first arg', expr = quote(out <- pow(exp(arg1), arg2)), inputDim = c(1,0), outputDim = 1)
+  ## list(name = 'modulo of vector and scalar', expr = quote(out <- arg1 %% arg2), inputDim = c(1,0), outputDim = 1) ## FAILS with Eigen casting
+  )
+
+testsMoreMath = list(
+  list(name = 'inverse cloglog of scalar', expr = quote(out <- icloglog(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'cloglog/inverse cloglog of scalar', expr = quote(out <- cloglog(icloglog(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'inverse logit of scalar', expr = quote(out <- ilogit(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'expit of scalar', expr = quote(out <- expit(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'logit/expit of scalar', expr = quote(out <- logit(expit(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'inverse probit of scalar', expr = quote(out <- iprobit(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'inverse probit of scalar via phi', expr = quote(out <- phi(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'probit/iprobit of scalar', expr = quote(out <- probit(iprobit(arg1))), inputDim = 0, outputDim = 0),
+  ###
+  list(name = 'ceiling of scalar', expr = quote(out <- ceiling(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'floor of scalar', expr = quote(out <- floor(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'round of scalar', expr = quote(out <- round(arg1)), inputDim = 0, outputDim = 0),
+  list(name = 'trunc of scalar', expr = quote(out <- trunc(arg1)), inputDim = 0, outputDim = 0),
+  ###
+  list(name = 'gamma of scalar', expr = quote(out <- gamma(abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'lgamma of scalar', expr = quote(out <- lgamma(abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'loggam of scalar', expr = quote(out <- loggam(abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'log1p of scalar', expr = quote(out <- log1p(abs(arg1))), inputDim = 0, outputDim = 0),
+  list(name = 'factorial of scalar', expr = quote(out <- factorial(ceiling(abs(arg1)))), inputDim = 0, outputDim = 0),
+  list(name = 'lfactorial of scalar', expr = quote(out <- lfactorial(ceiling(abs(arg1)))), inputDim = 0, outputDim = 0),
+  ###
+  list(name = 'inverse cloglog of vector', expr = quote(out <- icloglog(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'cloglog/inverse cloglog of vector', expr = quote(out <- cloglog(icloglog(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'inverse logit of vector', expr = quote(out <- ilogit(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'expit of vector', expr = quote(out <- expit(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'logit/expit of vector', expr = quote(out <- logit(expit(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'inverse probit of vector', expr = quote(out <- iprobit(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'inverse probit of vector via phi', expr = quote(out <- phi(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'probit/iprobit of vector', expr = quote(out <- probit(iprobit(arg1))), inputDim = 1, outputDim = 1),
+  ###
+  list(name = 'ceiling of vector', expr = quote(out <- ceiling(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'floor of vector', expr = quote(out <- floor(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'round of vector', expr = quote(out <- round(arg1)), inputDim = 1, outputDim = 1),
+  list(name = 'trunc of vector', expr = quote(out <- trunc(arg1)), inputDim = 1, outputDim = 1),
+  ###
+  list(name = 'gamma of vector', expr = quote(out <- gamma(abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'lgamma of vector', expr = quote(out <- lgamma(abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'loggam of vector', expr = quote(out <- loggam(abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'log1p of vector', expr = quote(out <- log1p(abs(arg1))), inputDim = 1, outputDim = 1),
+  list(name = 'factorial of vector', expr = quote(out <- factorial(ceiling(abs(arg1)))), inputDim = 1, outputDim = 1),
+  list(name = 'lfactorial of vector', expr = quote(out <- lfactorial(ceiling(abs(arg1)))), inputDim = 1, outputDim = 1)
+  )
+
+testsReduction = list(
+  ### vector
+  list(name = 'min of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  list(name = 'max of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  list(name = 'sum of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  list(name = 'mean of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  list(name = 'sd of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  list(name = 'var of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  list(name = 'prod of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
+  ## list(name = 'norm of vector', expr = quote(out <- norm(arg1)), inputDim = 1, outputDim = 0),  ## norm doesn't work on vector in R
+  ### matrix
+  list(name = 'min of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
+  list(name = 'max of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
+  list(name = 'sum of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
+  list(name = 'mean of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
+##  list(name = 'sd of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
+  list(name = 'var of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
+  list(name = 'prod of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0)
+##  list(name = 'norm of matrix', expr = quote(out <- norm(arg1)), inputDim = 2, outputDim = 0, Rcode = quote(out <- norm(arg1, "F"))) ## NIMBLE's C norm is apparently Frobenius, so R and C nimble functions differ => FAILS
+  )
+
+testsComparison = list(
+  ## scalar
+  list(name = 'greater than, scalar', expr = quote(out <- arg1 > arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'equals, scalar', expr = quote(out <- arg1 == arg2), inputDim = c(0,0), outputDim = 0),
+  list(name = 'not equals, scalar', expr = quote(out <- arg1 != arg2), inputDim = c(0,0), outputDim = 0),
+  ## vector
+  ## list(name = 'greater than, vector', expr = quote(out <- arg1 > arg2), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen issue
+  ## list(name = 'equals, vector', expr = quote(out <- arg1 == arg2), inputDim = c(1,1), outputDim = 1),  ## FAILS with Eigen issue
+  ## list(name = 'not equals, vector', expr = quote(out <- arg1 != arg2), inputDim = c(1,1), outputDim = 1),  ## FAILS with Eigen issue
+  ## logical
+  list(name = 'and operator, scalar', expr = quote(out <- arg1 & arg2), inputDim = c(0,0), outputDim = 0, logicalArgs = c(TRUE, TRUE)),
+  list(name = 'or operator, scalar', expr = quote(out <- arg1 | arg2), inputDim = c(0,0), outputDim = 0, logicalArgs = c(TRUE, TRUE)),
+  list(name = 'not operator, scalar', expr = quote(out <- !arg1), inputDim = c(0), outputDim = 0, logicalArgs = c(TRUE))
+)
+
+
+testsMatrix = list(
+    list(name = 'forwardsolve matrix-vector', expr = quote(out <- forwardsolve(arg1, arg2)), inputDim = c(2, 1), outputDim = 1),
+    list(name = 'forwardsolve matrix-matrix', expr = quote(out <- forwardsolve(arg1, arg2)), inputDim = c(2, 2), outputDim = 2),
+    list(name = 'backsolve matrix-vector', expr = quote(out <- backsolve(arg1, arg2)), inputDim = c(2, 1), outputDim = 1),
+    list(name = 'backsolve matrix-matrix', expr = quote(out <- backsolve(arg1, arg2)), inputDim = c(2, 2), outputDim = 2),
+
+    list(name = 'forwardsolve matrix-vector with indices', expr = quote(out <- forwardsolve(arg1[1:2,1:2], arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
+    list(name = 'forwardsolve matrix-matrix with indices', expr = quote(out <- forwardsolve(arg1[1:2,1:2], arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
+    list(name = 'backsolve matrix-vector with indices', expr = quote(out <- backsolve(arg1[1:2,1:2], arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
+    list(name = 'backsolve matrix-matrix with indices', expr = quote(out <- backsolve(arg1[1:2,1:2], arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
+
+    list(name = 'forwardsolve matrix-vector amid expr', expr = quote(out <- arg2[1:2] + forwardsolve(arg1[1:2,1:2], arg2[1:2] + arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
+    list(name = 'forwardsolve matrix-matrix amid expr', expr = quote(out <- arg2[1:2,1:2] + forwardsolve(arg1[1:2,1:2], arg2[1:2,1:2] + arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
+    list(name = 'backsolve matrix-vector amid expr', expr = quote(out <- arg2[1:2] + backsolve(arg1[1:2,1:2], arg2[1:2] + arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
+    list(name = 'backsolve matrix-matrix amid expr', expr = quote(out <- arg2[1:2,1:2] + backsolve(arg1[1:2,1:2], arg2[1:2,1:2] + arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
+
+    list(name = 'chol', expr = quote({ A <- arg1; for(i in 1:dim(A)[1]) A[i,i] <- A[i,i] + 10; out <- chol(A) }), inputDim = c(2), outputDim = 2),
+    list(name = 'matrix-vector multiply', expr = quote(out <- arg1 %*% arg2), inputDim = c(2, 1), outputDim = 2),
+    list(name = 'vector-matrix multiply', expr = quote(out <- t(arg1) %*% arg2), inputDim = c(1, 2), outputDim = 2),
+    list(name = 'matrix-matrix multiply', expr = quote(out <- arg1 %*% arg2), inputDim = c(2, 2), outputDim = 2)
+)
+
+
+
+
diff --git a/packages/nimble/inst/tests/test-math.R b/packages/nimble/inst/tests/test-math.R
index 13f9036cf..9541252e9 100644
--- a/packages/nimble/inst/tests/test-math.R
+++ b/packages/nimble/inst/tests/test-math.R
@@ -1,201 +1,12 @@
 ### INSTRUCTIONS:
+## got to mathTestLists.R:
 ## enter each test as a list, with an informative name, NIMBLE expression to evaluate, vector of input dimensions, value of output dimension, and (if NIMBLE expression cannot be directly evaluated in R) the equivalent pure R expression whose result should match the NIMBLE result
 
 source(system.file(file.path('tests', 'test_utils.R'), package = 'nimble'))
 
 context("Testing of math functions in NIMBLE code")
 
-testsVaried = list(
-  list(name = "matrix direct product", expr = quote(out <- arg1 * arg2), inputDim = c(2,2), outputDim = 2),
-  list(name = "matrix direct product with scalar addition", expr = quote(out <- (arg1+1) * (arg2+1)), inputDim = c(2,2), outputDim = 2),
-  list(name = "matrix absolute value", expr = quote(out <- abs(arg1)), inputDim = c(2), outputDim = 2),
-  list(name = "matrix absolute value with scalar addition", expr = quote(out <- abs(arg1 - 2)), inputDim = c(2), outputDim = 2),
-  list(name = "vector pmin", expr = quote(out <- pmin(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
-  list(name = "vector pmax", expr = quote(out <- pmax(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
-  list(name = "sd with addition", expr = quote(out <- sd(arg1) + 3), inputDim = c(1), outputDim = 0),
-  list(name = "sd of vector with addition", expr = quote(out <- sd(arg1 + 3)), inputDim = c(1), outputDim = 0),
-  list(name = "sd of matrix-vector multiply", expr = quote(out <- sd(arg1 %*% arg2)), inputDim = c(2,1), outputDim = 0),
-  list(name = "var of vector", expr = quote(out <- var(arg1)), inputDim = c(1), outputDim = 0),
-  list(name = "log determinant", expr = quote(out <- logdet(arg1)), inputDim = c(2), outputDim = 0)
-  )
-
-testsBasicMath = list(
-  list(name = 'exp of scalar', expr = quote(out <- exp(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'log of scalar', expr = quote(out <- log(abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'sqrt of scalar', expr = quote(out <- sqrt(abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'abs of scalar', expr = quote(out <- abs(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'step of scalar', expr = quote(out <- step(arg1)), inputDim = 0, outputDim = 0, Rcode = quote( out <- as.numeric(arg1 > 0))),
-  list(name = 'cube of scalar', expr = quote(out <- cube(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'cos of scalar', expr = quote(out <- cos(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'acos of cos of scalar', expr = quote(out <- acos(cos(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'sin of scalar', expr = quote(out <- sin(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'asin of sin of scalar', expr = quote(out <- asin(sin(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'tan of scalar', expr = quote(out <- tan(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'atan of tan of scalar', expr = quote(out <- atan(tan(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'cosh of scalar', expr = quote(out <- cosh(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'sinh of scalar', expr = quote(out <- sinh(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'tanh of scalar', expr = quote(out <- tanh(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'acosh of scalar', expr = quote(out <- acosh(1 + abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'asinh of scalar', expr = quote(out <- asinh(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'atanh of scalar', expr = quote(out <- atanh(abs(arg1)%%1)), inputDim = 0, outputDim = 0),
-  ###
-  list(name = 'exp of vector', expr = quote(out <- exp(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'log of vector', expr = quote(out <- log(abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'sqrt of vector', expr = quote(out <- sqrt(abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'abs of vector', expr = quote(out <- abs(arg1)), inputDim = 1, outputDim = 1),
-##  list(name = 'step of vector', expr = quote(out <- step(arg1)), inputDim = 1, outputDim = 1, Rcode = quote(out <- as.numeric(arg1 > 0))),   ## FAILS on compileNimble(nfR) with Eigen error
-  list(name = 'cube of vector', expr = quote(out <- cube(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'cos of vector', expr = quote(out <- cos(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'acos of cos of vector', expr = quote(out <- acos(cos(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'sin of vector', expr = quote(out <- sin(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'asin of sin of vector', expr = quote(out <- asin(sin(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'tan of vector', expr = quote(out <- tan(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'atan of tan of vector', expr = quote(out <- atan(tan(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'cosh of vector', expr = quote(out <- cosh(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'sinh of vector', expr = quote(out <- sinh(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'tanh of vector', expr = quote(out <- tanh(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'acosh of vector', expr = quote(out <- acosh(1 + abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'asinh of vector', expr = quote(out <- asinh(arg1)), inputDim = 1, outputDim = 1),
-##  list(name = 'atanh of vector', expr = quote(out <- atanh(arg1%%1)), inputDim = 1, outputDim = 1), ## FAILS - issue here is probably that modulo on vecs doesn't work but need to restrict domain for atanh
-  ###
-  list(name = 'sum of scalars', expr = quote(out <- arg1 + arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'diff of scalars', expr = quote(out <- arg1 - arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'product of scalars', expr = quote(out <- arg1 * arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'ratio of scalars', expr = quote(out <- arg1 / arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'power of scalars via ^', expr = quote(out <- arg1 ^ arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'power of scalars via pow', expr = quote(out <- pow(arg1, arg2)), inputDim = c(0,0), outputDim = 0),
-  list(name = 'power of scalars via ^ with positive first arg', expr = quote(out <- exp(arg1) ^ arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'power of scalars via pow with positive first arg', expr = quote(out <- pow(exp(arg1), arg2)), inputDim = c(0,0), outputDim = 0),
-  list(name = 'modulo of scalars', expr = quote(out <- arg1 %% arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'min of scalars', expr = quote(out <- min(arg1, arg2)), inputDim = c(0,0), outputDim = 0),
-  list(name = 'max of scalars', expr = quote(out <- max(arg1, arg2)), inputDim = c(0,0), outputDim = 0),
-  ###
-  list(name = 'sum of vectors', expr = quote(out <- arg1 + arg2), inputDim = c(1,1), outputDim = 1),
-  list(name = 'diff of vectors', expr = quote(out <- arg1 - arg2), inputDim = c(1,1), outputDim = 1),
-  list(name = 'product of vectors', expr = quote(out <- arg1 * arg2), inputDim = c(1,1), outputDim = 1),
-  list(name = 'ratio of vectors', expr = quote(out <- arg1 / arg2), inputDim = c(1,1), outputDim = 1),
-##  list(name = 'power of vectors via ^', expr = quote(out <- arg1 ^ arg2), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen casting
-  ## list(name = 'power of vectors via pow', expr = quote(out <- pow(arg1, arg2)), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen casting
-  ## list(name = 'modulo of vectors', expr = quote(out <- arg1 %% arg2), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen casting
-  list(name = 'pmin of vectors', expr = quote(out <- pmin(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
-  list(name = 'pmax of vectors', expr = quote(out <- pmax(arg1, arg2)), inputDim = c(1,1), outputDim = 1),
-  ###
-  list(name = 'sum of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'diff of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'product of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'ratio of vector and scalar', expr = quote(out <- arg1 + arg2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'power of vector and scalar via ^', expr = quote(out <- arg1 ^ arg2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'power of vector and scalar via pow', expr = quote(out <- pow(arg1, arg2)), inputDim = c(1,0), outputDim = 1),
-  list(name = 'power of vector and constant via ^', expr = quote(out <- arg1 ^ 2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'power of vector and constant via pow', expr = quote(out <- pow(arg1, 2)), inputDim = c(1,0), outputDim = 1),
-  list(name = 'power of vector and scalar via ^ with positive first arg', expr = quote(out <- exp(arg1) ^ arg2), inputDim = c(1,0), outputDim = 1),
-  list(name = 'power of vector and scalar via pow with positive first arg', expr = quote(out <- pow(exp(arg1), arg2)), inputDim = c(1,0), outputDim = 1)
-  ## list(name = 'modulo of vector and scalar', expr = quote(out <- arg1 %% arg2), inputDim = c(1,0), outputDim = 1) ## FAILS with Eigen casting
-  )
-
-testsMoreMath = list(
-  list(name = 'inverse cloglog of scalar', expr = quote(out <- icloglog(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'cloglog/inverse cloglog of scalar', expr = quote(out <- cloglog(icloglog(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'inverse logit of scalar', expr = quote(out <- ilogit(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'expit of scalar', expr = quote(out <- expit(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'logit/expit of scalar', expr = quote(out <- logit(expit(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'inverse probit of scalar', expr = quote(out <- iprobit(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'inverse probit of scalar via phi', expr = quote(out <- phi(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'probit/iprobit of scalar', expr = quote(out <- probit(iprobit(arg1))), inputDim = 0, outputDim = 0),
-  ###
-  list(name = 'ceiling of scalar', expr = quote(out <- ceiling(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'floor of scalar', expr = quote(out <- floor(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'round of scalar', expr = quote(out <- round(arg1)), inputDim = 0, outputDim = 0),
-  list(name = 'trunc of scalar', expr = quote(out <- trunc(arg1)), inputDim = 0, outputDim = 0),
-  ###
-  list(name = 'gamma of scalar', expr = quote(out <- gamma(abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'lgamma of scalar', expr = quote(out <- lgamma(abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'loggam of scalar', expr = quote(out <- loggam(abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'log1p of scalar', expr = quote(out <- log1p(abs(arg1))), inputDim = 0, outputDim = 0),
-  list(name = 'factorial of scalar', expr = quote(out <- factorial(ceiling(abs(arg1)))), inputDim = 0, outputDim = 0),
-  list(name = 'lfactorial of scalar', expr = quote(out <- lfactorial(ceiling(abs(arg1)))), inputDim = 0, outputDim = 0),
-  ###
-  list(name = 'inverse cloglog of vector', expr = quote(out <- icloglog(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'cloglog/inverse cloglog of vector', expr = quote(out <- cloglog(icloglog(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'inverse logit of vector', expr = quote(out <- ilogit(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'expit of vector', expr = quote(out <- expit(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'logit/expit of vector', expr = quote(out <- logit(expit(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'inverse probit of vector', expr = quote(out <- iprobit(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'inverse probit of vector via phi', expr = quote(out <- phi(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'probit/iprobit of vector', expr = quote(out <- probit(iprobit(arg1))), inputDim = 1, outputDim = 1),
-  ###
-  list(name = 'ceiling of vector', expr = quote(out <- ceiling(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'floor of vector', expr = quote(out <- floor(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'round of vector', expr = quote(out <- round(arg1)), inputDim = 1, outputDim = 1),
-  list(name = 'trunc of vector', expr = quote(out <- trunc(arg1)), inputDim = 1, outputDim = 1),
-  ###
-  list(name = 'gamma of vector', expr = quote(out <- gamma(abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'lgamma of vector', expr = quote(out <- lgamma(abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'loggam of vector', expr = quote(out <- loggam(abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'log1p of vector', expr = quote(out <- log1p(abs(arg1))), inputDim = 1, outputDim = 1),
-  list(name = 'factorial of vector', expr = quote(out <- factorial(ceiling(abs(arg1)))), inputDim = 1, outputDim = 1),
-  list(name = 'lfactorial of vector', expr = quote(out <- lfactorial(ceiling(abs(arg1)))), inputDim = 1, outputDim = 1)
-  )
-
-testsReduction = list(
-  ### vector
-  list(name = 'min of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  list(name = 'max of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  list(name = 'sum of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  list(name = 'mean of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  list(name = 'sd of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  list(name = 'var of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  list(name = 'prod of vector', expr = quote(out <- min(arg1)), inputDim = 1, outputDim = 0),
-  ## list(name = 'norm of vector', expr = quote(out <- norm(arg1)), inputDim = 1, outputDim = 0),  ## norm doesn't work on vector in R
-  ### matrix
-  list(name = 'min of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
-  list(name = 'max of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
-  list(name = 'sum of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
-  list(name = 'mean of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
-##  list(name = 'sd of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
-  list(name = 'var of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0),
-  list(name = 'prod of matrix', expr = quote(out <- min(arg1)), inputDim = 2, outputDim = 0)
-##  list(name = 'norm of matrix', expr = quote(out <- norm(arg1)), inputDim = 2, outputDim = 0, Rcode = quote(out <- norm(arg1, "F"))) ## NIMBLE's C norm is apparently Frobenius, so R and C nimble functions differ => FAILS
-  )
-
-testsComparison = list(
-  ## scalar
-  list(name = 'greater than, scalar', expr = quote(out <- arg1 > arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'equals, scalar', expr = quote(out <- arg1 == arg2), inputDim = c(0,0), outputDim = 0),
-  list(name = 'not equals, scalar', expr = quote(out <- arg1 != arg2), inputDim = c(0,0), outputDim = 0),
-  ## vector
-  ## list(name = 'greater than, vector', expr = quote(out <- arg1 > arg2), inputDim = c(1,1), outputDim = 1), ## FAILS with Eigen issue
-  ## list(name = 'equals, vector', expr = quote(out <- arg1 == arg2), inputDim = c(1,1), outputDim = 1),  ## FAILS with Eigen issue
-  ## list(name = 'not equals, vector', expr = quote(out <- arg1 != arg2), inputDim = c(1,1), outputDim = 1),  ## FAILS with Eigen issue
-  ## logical
-  list(name = 'and operator, scalar', expr = quote(out <- arg1 & arg2), inputDim = c(0,0), outputDim = 0, logicalArgs = c(TRUE, TRUE)),
-  list(name = 'or operator, scalar', expr = quote(out <- arg1 | arg2), inputDim = c(0,0), outputDim = 0, logicalArgs = c(TRUE, TRUE)),
-  list(name = 'not operator, scalar', expr = quote(out <- !arg1), inputDim = c(0), outputDim = 0, logicalArgs = c(TRUE))
-)
-
-
-testsMatrix = list(
-    list(name = 'forwardsolve matrix-vector', expr = quote(out <- forwardsolve(arg1, arg2)), inputDim = c(2, 1), outputDim = 1),
-    list(name = 'forwardsolve matrix-matrix', expr = quote(out <- forwardsolve(arg1, arg2)), inputDim = c(2, 2), outputDim = 2),
-    list(name = 'backsolve matrix-vector', expr = quote(out <- backsolve(arg1, arg2)), inputDim = c(2, 1), outputDim = 1),
-    list(name = 'backsolve matrix-matrix', expr = quote(out <- backsolve(arg1, arg2)), inputDim = c(2, 2), outputDim = 2),
-
-    list(name = 'forwardsolve matrix-vector with indices', expr = quote(out <- forwardsolve(arg1[1:2,1:2], arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
-    list(name = 'forwardsolve matrix-matrix with indices', expr = quote(out <- forwardsolve(arg1[1:2,1:2], arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
-    list(name = 'backsolve matrix-vector with indices', expr = quote(out <- backsolve(arg1[1:2,1:2], arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
-    list(name = 'backsolve matrix-matrix with indices', expr = quote(out <- backsolve(arg1[1:2,1:2], arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
-
-    list(name = 'forwardsolve matrix-vector amid expr', expr = quote(out <- arg2[1:2] + forwardsolve(arg1[1:2,1:2], arg2[1:2] + arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
-    list(name = 'forwardsolve matrix-matrix amid expr', expr = quote(out <- arg2[1:2,1:2] + forwardsolve(arg1[1:2,1:2], arg2[1:2,1:2] + arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
-    list(name = 'backsolve matrix-vector amid expr', expr = quote(out <- arg2[1:2] + backsolve(arg1[1:2,1:2], arg2[1:2] + arg2[1:2])), inputDim = c(2, 1), outputDim = 1),
-    list(name = 'backsolve matrix-matrix amid expr', expr = quote(out <- arg2[1:2,1:2] + backsolve(arg1[1:2,1:2], arg2[1:2,1:2] + arg2[1:2,1:2])), inputDim = c(2, 2), outputDim = 2),
-
-    list(name = 'chol', expr = quote({ A <- arg1; for(i in 1:dim(A)[1]) A[i,i] <- A[i,i] + 10; out <- chol(A) }), inputDim = c(2), outputDim = 2),
-    list(name = 'matrix-vector multiply', expr = quote(out <- arg1 %*% arg2), inputDim = c(2, 1), outputDim = 2),
-    list(name = 'vector-matrix multiply', expr = quote(out <- t(arg1) %*% arg2), inputDim = c(1, 2), outputDim = 2),
-    list(name = 'matrix-matrix multiply', expr = quote(out <- arg1 %*% arg2), inputDim = c(2, 2), outputDim = 2)
-)
-
+source(system.file(file.path('tests', 'mathTestLists.R'), package = 'nimble'))
 
 set.seed(0)
 ans1 <- sapply(testsVaried, test_math)    ## 12
diff --git a/packages/nimble/inst/tests/test-refactorCompilationSteps.R b/packages/nimble/inst/tests/test-refactorCompilationSteps.R
new file mode 100644
index 000000000..492721ee8
--- /dev/null
+++ b/packages/nimble/inst/tests/test-refactorCompilationSteps.R
@@ -0,0 +1,59 @@
+source(system.file(file.path('tests', 'test_utils.R'), package = 'nimble'))
+
+context("Testing of old vs. new generated C++ during refactoring steps")
+
+oldWarnLevel <- options('warn')
+options(warn = -1)
+
+## Known concern: ordering of asRow()/asCol() and intermediates
+
+compareOldAndNewCompilationRC <- function(input) {
+    run <- input$run
+    name <- input$name
+    require(testthat)
+
+    foo <- nimbleFunction(run = run, name = 'foo')
+    nimbleOptions(useRefactoredSizeProcessing = FALSE)
+    nimble:::resetLabelFunctionCreators() ## sets any generated IDs back to 1
+    testProject <- nimble:::nimbleProjectClass(name = 'for_comparison')
+    ## management of compilation through the control list is crude and leads to error we don't care about
+    ## but the project does contain the result, so we can run this, catch the error to keep running
+    ## and extract the result from the project
+    old <- try(compileNimble(foo, project = testProject, control = list(writeFiles = TRUE, compileCpp = FALSE, loadSO = FALSE)))
+    filename <- testProject$RCfunInfos[['foo']][['cppClass']]$filename
+    newfilename <- paste0(filename,'_original')
+    pathedfilename <- file.path(tempdir(), 'nimble_generatedCode', filename)
+    original_pathedfilename <- file.path(tempdir(), 'nimble_generatedCode', newfilename)
+    for(ext in c('.h', '.cpp')) file.copy(paste0(pathedfilename, ext), paste0(original_pathedfilename, ext), overwrite = TRUE)
+    
+    nimbleOptions(useRefactoredSizeProcessing = TRUE)
+    nimble:::resetLabelFunctionCreators() ## sets any generated IDs back to 1
+    testProject <- nimble:::nimbleProjectClass(name = 'for_comparison')
+    new <- try(compileNimble(foo, project = testProject, control = list(writeFiles = TRUE, compileCpp = FALSE, loadSO = FALSE)))
+    filename <- testProject$RCfunInfos[['foo']][['cppClass']]$filename
+    ## we could regenerate it, but might as well read it from the file
+    refactored_pathedfilename <- file.path(tempdir(), 'nimble_generatedCode', filename)
+
+    for(ext in c('.h','.cpp'))
+        compareFilesUsingDiff(paste0(refactored_pathedfilename, ext), paste0(original_pathedfilename, ext),
+                              main = paste0(ext, ' files do not match for: ', name))    
+}
+
+testCases <- list(
+    list(name = 'Y <- x',
+         run = function(x = double(1)) {
+             Y <- x
+         }))
+
+ans <- lapply(testCases, compareOldAndNewCompilationRC)
+
+compareOldAndNewMathTest <- function(input) {
+    runFun <- gen_runFun(input)
+    input$run <- runFun
+    compareOldAndNewCompilationRC(input)
+}
+
+source(system.file(file.path('tests', 'mathTestLists.R'), package = 'nimble'))
+ans2 <- lapply(testsBasicMath, compareOldAndNewMathTest)
+
+options(warn = as.numeric(oldWarnLevel))
diff --git a/packages/nimble/inst/tests/test_utils.R b/packages/nimble/inst/tests/test_utils.R
index 05f64739b..80e139c03 100644
--- a/packages/nimble/inst/tests/test_utils.R
+++ b/packages/nimble/inst/tests/test_utils.R
@@ -897,3 +897,12 @@ compareFilesByLine <- function(trialResults, correctResults, main = "") {
     }, 1:linesToTest, trialResults, correctResults)
     invisible(NULL)
 }
+
+compareFilesUsingDiff <- function(trialFile, correctFile, main = "") {
+    if(main == "") main <- paste0(trialFile, ' and ', correctFile, ' do not match\n')
+    diffOutput <- system2('diff', c(trialFile, correctFile), stdout = TRUE)
+    test_that(paste0(main, paste0(diffOutput, collapse = '\n')),
+              expect_true(length(diffOutput) == 0)
+              )
+    invisible(NULL)
+}