R/nimbleFunction_keywordProcessing.R

Defines functions nimDerivsInfoClass_update_init_impl makeModelDerivsInfo_impl getImmediateParentNodes makeModelDerivsInfo nimDerivsInfoClass_output_init_impl makeOutputNodes nimDerivsInfoClass_init_impl handleScaleAndRateForExpNimble handleScaleAndRateForGamma makeSingleIndexAccessCodeNames makeMapSetupCodeNames varAndIndices2flatIndex getVarAndIndices varAndIndices2mapParts determineNdimFromOneCase makeMapAccessExpr makeSingleIndexAccessExpr matchKeywords_recurse matchKeywordCode matchKeywordCodeMemberFun determineNdimsFromNfproc matchAndFill.call isSymbolType symTypeFromSymTab getSymObj_recurse addNecessarySetupCode addNewCode addBlockFromCppName addSetupCodeNames isSetupCodeGenerated isCodeArgBlank code2Name_fromArgList processKeywords_recurse processKeywordCodeMemberFun processKeyword addDistList2matchFunctions stripArgs ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ] ]

Documented in makeModelDerivsInfo

######	KEYWORD PROCESSING OBJECTS



###		CLASSES

# keywordInfoClass is a class which contains the processor for each keyword
keywordInfoClass <- setRefClass('keywordInfoClass',
                                fields = list(
                                    keyword = 'ANY',
                                    processor = 'ANY'))


# setupCodeTemplateClass is a class that contains the template for generating 
# new setupCode. Objects of this class are used by the function addNecessarySetupCode
setupCodeTemplateClass <- setRefClass('setupCodeTemplateClass',
                                      fields = list(
                                          makeName = 'ANY',
                                          codeTemplate = 'ANY',
                                          makeCodeSubList = 'ANY',
                                          makeOtherNames = 'ANY'),
                                          methods = list(
                                          initialize = function(...){
                                          	makeOtherNames <<- function(name, argList)	return(character(0))
                                          	callSuper(...)
                                          }
                                          ) )


### KEYWORD INFO OBJECTS
		
d_gamma_keywordInfo <- keywordInfoClass(
    keyword = 'dgamma',
    processor = function(code, nfProc, RCfunProc){
        code <- handleScaleAndRateForGamma(code)
	return(code)
    }) 

pq_gamma_keywordInfo <- keywordInfoClass(
    keyword = 'pq_gamma',
    processor = function(code, nfProc, RCfunProc){
        code <- handleScaleAndRateForGamma(code)
	return(code)
    })

rgamma_keywordInfo <- keywordInfoClass(
    keyword = 'rgamma',
    processor = function(code, nfProc, RCfunProc){
        code <- handleScaleAndRateForGamma(code)
        return(code)
    }
)

d_exp_nimble_keywordInfo <- keywordInfoClass(
	keyword = 'dexp_nimble',
	processor = function(code, nfProc, RCfunProc){
		code <- handleScaleAndRateForExpNimble(code)
	return(code)
	}) 

pq_exp_nimble_keywordInfo <- keywordInfoClass(
	keyword = 'pq_exp_nimble',
	processor = function(code, nfProc, RCfunProc){
		code <- handleScaleAndRateForExpNimble(code)
	return(code)
})

rexp_nimble_keywordInfo <- keywordInfoClass(
	keyword = 'rexp_nimble',
	processor = function(code, nfProc, RCfunProc){
		code <- handleScaleAndRateForExpNimble(code)
		return(code)
	}
)

besselK_keywordInfo <- keywordInfoClass(
    keyword = 'besselK',
    processor = function(code, nfProc, RCfunProc) {
        expon.scaledArg <- code$expon.scaled
        if(is.null(expon.scaledArg))
            expon.scaledArg <- FALSE
        if(is.numeric(expon.scaledArg) || is.logical(expon.scaledArg)) {
            code$expon.scaled <- 1 + as.logical(expon.scaledArg)
        } else code$expon.scaled <- substitute(1 + A, list(A = expon.scaledArg))
        return(code)
    }
)


nimSeq_keywordInfo <- keywordInfoClass(
    keyword = 'nimSeq',
    processor = function(code, nfProc, RCfunProc) {
        useBy <- !isCodeArgBlank(code, 'by')
        useLen <- !isCodeArgBlank(code, 'length.out')
        if(useBy && useLen)
            newRunCode <- substitute(nimSeqByLen(FROM, 0, BY, LEN), list(FROM = code$from, BY = code$by, LEN = code$length.out))
        else {
            if(useLen) {
                newRunCode <- substitute(nimSeqLen(FROM, TO, 0, LEN), list(FROM = code$from, TO = code$to, LEN = code$length.out))
            } else {
                byVal <- if(useBy) code$by else 1 ## default by = 1
                newRunCode <- substitute(nimSeqBy(FROM, TO, BY, 0), list(FROM = code$from, TO = code$to, BY = byVal))
            }
        }
        return(newRunCode)
    }
)
    

values_keywordInfo <- keywordInfoClass(
    keyword = 'values',
    processor = function(code, nfProc, RCfunProc){
      if(!isCodeArgBlank(code, 'accessor'))
      	return(code)
      if(isCodeArgBlank(code, 'model'))
      	stop('model argument missing from values call, with no accessor argument supplied')
      
      accessArgList <- list(model = code$model, nodes = code$nodes, logProb = FALSE, logProbOnly = FALSE)

      useAccessorVectorByIndex <- FALSE
      if(hasBracket(accessArgList$nodes)) { 
            useAccessorVectorByIndex <- TRUE
            if(length(accessArgList$nodes) != 3) stop(paste0('Problem with ', deparse(code),'. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- accessArgList$nodes[[3]]
            accessArgList$nodes <- accessArgList$nodes[[2]]
            accessArgList$sortUnique <- FALSE   
      }

      accessName <- modelVariableAccessorVector_setupCodeTemplate$makeName(accessArgList)
      addNecessarySetupCode(accessName, accessArgList, modelVariableAccessorVector_setupCodeTemplate, nfProc)
      if(!useAccessorVectorByIndex)
          newRunCode <- substitute(values(accessor = ACCESS_NAME), 
                               list(ACCESS_NAME = as.name(accessName)))
        else
            newRunCode <- substitute(values(accessor = ACCESS_NAME, accessorIndex = ACCESSVECINDEX),
                                 list(ACCESS_NAME = as.name(accessName), ACCESSVECINDEX = nodesIndexExpr))
      return(newRunCode)
    })                                    

getParam_keywordInfo <- keywordInfoClass(
    keyword = 'getParam',
    processor = function(code, nfProc, RCfunProc) {
        if(!isCodeArgBlank(code, 'nodeFunction'))
            return(code)
        errorContext <- deparse(code)
        nodeFunVec_ArgList <- list(model = code$model, nodes = code$node, includeData = TRUE, sortUnique = TRUE, errorContext = errorContext)
        if(!isCodeArgBlank(code, 'nodeFunctionIndex')) { ## new case: calculate(myNodeFunctionVector, nodeFunctionIndex = i), if myNodeFunctionVector was hand-created in setup code
            if(!isCodeArgBlank(code, 'nodes'))
                stop('nodes argument cannot be provided to getParam if nodeFunctionIndex is specified')
            return(code) ## no modification needed!
        }

        if(isCodeArgBlank(code, 'model'))
            stop('model argument missing from getParam, with no accessor argument supplied')
        if(isCodeArgBlank(code, 'node'))
            stop('node argument missing from getParam, with no accessor argument supplied')

        useNodeFunctionVectorByIndex <- FALSE
        if(hasBracket(nodeFunVec_ArgList$nodes)) { ## like calculate(model, nodes[i]), which could have started as model$calculate(nodes[i])
            useNodeFunctionVectorByIndex <- TRUE
            if(length(nodeFunVec_ArgList$nodes) != 3) stop(paste0('Problem with ', deparse(code),'. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- nodeFunVec_ArgList$nodes[[3]]
            nodeFunVec_ArgList$nodes <- nodeFunVec_ArgList$nodes[[2]]
            nodeFunVec_ArgList$sortUnique <- FALSE
        }

        nodeFunName <- nodeFunctionVector_SetupTemplate$makeName(nodeFunVec_ArgList)

        if(isCodeArgBlank(code, 'param'))
            stop("'param' argument missing from 'getParam', with no accessor argument supplied")

        ## Avoid situations where user has syntax from run code in `param` or of getting anything from global env (issue 1344).
        paramVars <- all.vars(code$param)
        wh <- which(!paramVars %in% nfProc$setupSymTab$getSymbolNames())
        if(length(wh))
            stop("`param` argument in `getParam` contains variables not found in setup code: ", paste(paramVars[wh], collapse = ", "), ".")
            
        paramInfo_ArgList <- list(model = code$model, node = nodeFunVec_ArgList$nodes, param = code$param, hasIndex = useNodeFunctionVectorByIndex) ## use nodeFunVec_ArgList$nodes instead of code$node because nodeFunVec_ArgList$nodes may have been updated if code$nodes has a run-time index.  In that case the paramID will be vector

        paramInfoName <- paramInfo_SetupTemplate$makeName(paramInfo_ArgList)
        paramIDname <- paramInfo_SetupTemplate$makeOtherNames(paramInfoName, paramInfo_ArgList)

        addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_SetupTemplate, nfProc)
        addNecessarySetupCode(paramInfoName, paramInfo_ArgList, paramInfo_SetupTemplate, nfProc)
        if(!useNodeFunctionVectorByIndex)
            newRunCode <- substitute(getParam(nodeFunction = NODEFUNVEC_NAME, paramID = PARAMID_NAME, paramInfo = PARAMINFO_NAME),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName), PARAMID_NAME = as.name(paramIDname), PARAMINFO_NAME = as.name(paramInfoName)))
        else
            newRunCode <- substitute(getParam(nodeFunction = NODEFUNVEC_NAME, paramID = PARAMID_NAME, paramInfo = PARAMINFO_NAME, nodeFunctionIndex = NODEFUNVECINDEX),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName), PARAMID_NAME = as.name(paramIDname), PARAMINFO_NAME = as.name(paramInfoName), NODEFUNVECINDEX = nodesIndexExpr))
        
        return(newRunCode)
    }
)

getBound_keywordInfo <- keywordInfoClass(
    keyword = 'getBound',
    processor = function(code, nfProc, RCfunProc) {
        if(!isCodeArgBlank(code, 'nodeFunction'))
            return(code)
        errorContext <- deparse(code)
        nodeFunVec_ArgList <- list(model = code$model, nodes = code$node, includeData = TRUE, sortUnique = TRUE, errorContext = errorContext)
        if(!isCodeArgBlank(code, 'nodeFunctionIndex')) { ## new case: calculate(myNodeFunctionVector, nodeFunctionIndex = i), if myNodeFunctionVector was hand-created in setup code
            if(!isCodeArgBlank(code, 'nodes'))
                stop('nodes argument cannot be provided to getParam if nodeFunctionIndex is specified')
            return(code) ## no modification needed!
        }

        if(isCodeArgBlank(code, 'model'))
            stop('model argument missing from getParam, with no accessor argument supplied')
        if(isCodeArgBlank(code, 'node'))
            stop('node argument missing from getParam, with no accessor argument supplied')

        useNodeFunctionVectorByIndex <- FALSE
        if(hasBracket(nodeFunVec_ArgList$nodes)) { ## like calculate(model, nodes[i]), which could have started as model$calculate(nodes[i])
            useNodeFunctionVectorByIndex <- TRUE
            if(length(nodeFunVec_ArgList$nodes) != 3) stop(paste0('Problem with ', deparse(code),'. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- nodeFunVec_ArgList$nodes[[3]]
            nodeFunVec_ArgList$nodes <- nodeFunVec_ArgList$nodes[[2]]
            nodeFunVec_ArgList$sortUnique <- FALSE
        }

        nodeFunName <- nodeFunctionVector_SetupTemplate$makeName(nodeFunVec_ArgList)

        if(isCodeArgBlank(code, 'bound'))
            stop("'bound' argument missing from 'getBound', with no accessor argument supplied")
        boundInfo_ArgList <- list(model = code$model, node = nodeFunVec_ArgList$nodes, bound = code$bound) ## use nodeFunVec_ArgList$nodes instead of code$node because nodeFunVec_ArgList$nodes may have been updated if code$nodes has a run-time index.  In that case the boundID will be vector
        boundInfoName <- boundInfo_SetupTemplate$makeName(boundInfo_ArgList)
        boundIDname <- boundInfo_SetupTemplate$makeOtherNames(boundInfoName, boundInfo_ArgList)

        addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_SetupTemplate, nfProc)
        addNecessarySetupCode(boundInfoName, boundInfo_ArgList, boundInfo_SetupTemplate, nfProc)
        if(!useNodeFunctionVectorByIndex)
            newRunCode <- substitute(getBound(nodeFunction = NODEFUNVEC_NAME, boundID = BOUNDID_NAME, boundInfo = BOUNDINFO_NAME),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName), BOUNDID_NAME = as.name(boundIDname), BOUNDINFO_NAME = as.name(boundInfoName)))
        else
            newRunCode <- substitute(getBound(nodeFunction = NODEFUNVEC_NAME, boundID = BOUNDID_NAME, boundInfo = BOUNDINFO_NAME, nodeFunctionIndex = NODEFUNVECINDEX),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName), BOUNDID_NAME = as.name(boundIDname), BOUNDINFO_NAME = as.name(boundInfoName), NODEFUNVECINDEX = nodesIndexExpr))
                                               
        return(newRunCode)
    }
)


calculate_keywordInfo <- keywordInfoClass(
    keyword = 'calculate',
    processor = function(code, nfProc, RCfunProc){
        if(!isCodeArgBlank(code, 'nodeFxnVector'))
            return(code)
        errorContext <- deparse(code)

        buildDerivs <- FALSE
        withDerivsOutputOnly <- FALSE
        if(isTRUE(getNimbleOption("enableDerivs"))) {
          derivControl <- environment(nfProc$nfGenerator)$buildDerivs[[RCfunProc$name]]
          ## There are two cases.
          ## If buildDerivs has an entry, then derivatives are enabled either
          ## explicitly or because of a nimDerivs(model$calculate), i.e. a direct line into derivs for just a model calculation.
          ##
          ## These need to be disambiguated for use below to choose the right setup template case
          buildDerivs <- !is.null(derivControl) | !is.null(code$wrt)
          withDerivsOutputOnly <- buildDerivs & is.null(code$wrt) ## This is the case of *not* nimDerivs(model$calculate)
        }
        
        nodeFunVec_ArgList <- list(model = code$model, nodes = code$nodes, wrtNodes = code$wrt,
                                   includeData = TRUE, sortUnique = TRUE, errorContext = errorContext)

        ##
        if(!isCodeArgBlank(code, 'nodeFunctionIndex')) { ## new case: calculate(myNodeFunctionVector, nodeFunctionIndex = i), if myNodeFunctionVector was hand-created in setup code
            if(!isCodeArgBlank(code, 'nodes'))
                stop('nodes argument cannot be provided to calculate if nodeFunctionIndex is specified')
            return(code) ## no modification needed!
        }
        if(isCodeArgBlank(code, 'model'))
            stop('model argument missing from calculate, with no accessor argument supplied')
        if(isCodeArgBlank(code, 'nodes')){
            LHSnodes_ArgList <- list(model = code$model)
            LHSnodes_name <- allLHSNodes_SetupTemplate$makeName(LHSnodes_ArgList)
            addNecessarySetupCode(LHSnodes_name, LHSnodes_ArgList, allLHSNodes_SetupTemplate, nfProc, allowToCpp = FALSE)
            nodeFunVec_ArgList$nodes = as.name(LHSnodes_name)
        }
        useNodeFunctionVectorByIndex <- FALSE
        if(hasBracket(nodeFunVec_ArgList$nodes)) { ## like calculate(model, nodes[i]), which could have started as model$calculate(nodes[i])
            if(buildDerivs)
                stop(paste0('Derivatives for an indexed node are not allowed.\n',
                            'Use nimDerivs(model$calculate(nodes), ...) instead of \n',
                            'nimDerivs(model$calculate(nodes[i]).\n'))
            useNodeFunctionVectorByIndex <- TRUE
            if(length(nodeFunVec_ArgList$nodes) != 3)
                stop(paste0('Problem with ',
                            deparse(code),
                            '. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- nodeFunVec_ArgList$nodes[[3]]
            nodeFunVec_ArgList$nodes <- nodeFunVec_ArgList$nodes[[2]]
            nodeFunVec_ArgList$sortUnique <- FALSE
        }

        if(!withDerivsOutputOnly) { ## This is regular mode, including without derivs at all and with buildDerivs but not nimDerivs(model$calculate...)
          nodeFunName <- nodeFunctionVector_SetupTemplate$makeName(nodeFunVec_ArgList)	
          addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_SetupTemplate, nfProc)
        } else {
          nodeFunName <- nodeFunctionVector_WithDerivsOutput_SetupTemplate$makeName(nodeFunVec_ArgList)	
          addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_WithDerivsOutput_SetupTemplate, nfProc)
        }
        if(!useNodeFunctionVectorByIndex){
            newRunCode <- substitute(calculate(nodeFxnVector = NODEFUNVEC_NAME),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName)))
        }
        else{
            newRunCode <- substitute(
                calculate(nodeFxnVector = NODEFUNVEC_NAME,
                          nodeFunctionIndex = NODEFUNVECINDEX),
                list(NODEFUNVEC_NAME = as.name(nodeFunName),
                     NODEFUNVECINDEX = nodesIndexExpr))
        }
        return(newRunCode)
    }
)

calculateDiff_keywordInfo <- keywordInfoClass(
    keyword = 'calculateDiff',
    processor = function(code, nfProc, RCfunProc){
        if(!isCodeArgBlank(code, 'nodeFxnVector'))
            return(code)
        errorContext <- deparse(code)
        nodeFunVec_ArgList <- list(model = code$model, nodes = code$nodes, includeData = TRUE, sortUnique = TRUE, errorContext = errorContext)
        if(!isCodeArgBlank(code, 'nodeFunctionIndex')) { ## new case: calculate(myNodeFunctionVector, nodeFunctionIndex = i), if myNodeFunctionVector was hand-created in setup code
            if(!isCodeArgBlank(code, 'nodes'))
                stop('nodes argument cannot be provided to calculateDiff if nodeFunctionIndex is specified')
            return(code) ## no modification needed!
        }
        
        if(isCodeArgBlank(code, 'model'))
            stop('model argument missing from calculateDiff, with no accessor argument supplied')
        if(isCodeArgBlank(code, 'nodes')){
            LHSnodes_ArgList <- list(model = code$model)
            LHSnodes_name <- allLHSNodes_SetupTemplate$makeName(LHSnodes_ArgList)
            addNecessarySetupCode(LHSnodes_name, LHSnodes_ArgList, allLHSNodes_SetupTemplate, nfProc, allowToCpp = FALSE)
            nodeFunVec_ArgList$nodes = as.name(LHSnodes_name)
        }
        useNodeFunctionVectorByIndex <- FALSE
        if(hasBracket(nodeFunVec_ArgList$nodes)) { ## like calculate(model, nodes[i]), which could have started as model$calculate(nodes[i])
            useNodeFunctionVectorByIndex <- TRUE
            if(length(nodeFunVec_ArgList$nodes) != 3) stop(paste0('Problem with ', deparse(code),'. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- nodeFunVec_ArgList$nodes[[3]]
            nodeFunVec_ArgList$nodes <- nodeFunVec_ArgList$nodes[[2]]
            nodeFunVec_ArgList$sortUnique <- FALSE
        }
        
        nodeFunName <- nodeFunctionVector_SetupTemplate$makeName(nodeFunVec_ArgList)	
        addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_SetupTemplate, nfProc)
        if(!useNodeFunctionVectorByIndex)
            newRunCode <- substitute(calculateDiff(nodeFxnVector = NODEFUNVEC_NAME),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName)))
        else
            newRunCode <- substitute(calculateDiff(nodeFxnVector = NODEFUNVEC_NAME, nodeFunctionIndex = NODEFUNVECINDEX),
                                 list(NODEFUNVEC_NAME = as.name(nodeFunName), NODEFUNVECINDEX = nodesIndexExpr))
        return(newRunCode)	
    }
)


simulate_keywordInfo <- keywordInfoClass(
    keyword = 'simulate',
    processor = function(code, nfProc, RCfunProc){
        if(!isCodeArgBlank(code, 'nodeFxnVector')){
            return(substitute(simulate(nodeFxnVector = NODEFXNVECTOR), list(NODEFXNVECTOR = code$nodeFxnVector) ) )
        }
        if(!isCodeArgBlank(code, 'INDEXEDNODEINFO_'))
            return(code)
        errorContext <- deparse(code)
        nodeFunVec_ArgList <- list(model = code$model, nodes = code$nodes, includeData = code$includeData, sortUnique = TRUE, errorContext = errorContext)
        if(!isCodeArgBlank(code, 'nodeFunctionIndex')) { ## new case: calculate(myNodeFunctionVector, nodeFunctionIndex = i), if myNodeFunctionVector was hand-created in setup code
            if(!isCodeArgBlank(code, 'nodes'))
                stop('nodes argument cannot be provided to simulate if nodeFunctionIndex is specified')
            return(code) ## no modification needed!
        }
        if(isCodeArgBlank(code, 'model'))
            stop('model argument missing from simulate, with no accessor argument supplied')
        if(isCodeArgBlank(code, 'nodes')){
            LHSnodes_ArgList <- list(model = code$model)
            LHSnodes_name <- allLHSNodes_SetupTemplate$makeName(LHSnodes_ArgList)
            addNecessarySetupCode(LHSnodes_name, LHSnodes_ArgList, allLHSNodes_SetupTemplate, nfProc, allowToCpp = FALSE)
            nodeFunVec_ArgList$nodes = as.name(LHSnodes_name)
        }
        useNodeFunctionVectorByIndex <- FALSE
        if(hasBracket(nodeFunVec_ArgList$nodes)) { ## like calculate(model, nodes[i]), which could have started as model$calculate(nodes[i])
            useNodeFunctionVectorByIndex <- TRUE
            if(length(nodeFunVec_ArgList$nodes) != 3) stop(paste0('Problem with ', deparse(code),'. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- nodeFunVec_ArgList$nodes[[3]]
            nodeFunVec_ArgList$nodes <- nodeFunVec_ArgList$nodes[[2]]
            nodeFunVec_ArgList$sortUnique <- FALSE # If includeData = FALSE, this can trigger error from nodeFunctionVector if nodes does contain data
        }
  
        nodeFunName <- nodeFunctionVector_SetupTemplate$makeName(nodeFunVec_ArgList)	
        addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_SetupTemplate, nfProc)
        if(!useNodeFunctionVectorByIndex)
            newRunCode <- substitute(simulate(nodeFxnVector = NODEFUNVEC_NAME),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName)))
        else
            newRunCode <- substitute(simulate(nodeFxnVector = NODEFUNVEC_NAME, nodeFunctionIndex = NODEFUNVECINDEX),
                                 list(NODEFUNVEC_NAME = as.name(nodeFunName), NODEFUNVECINDEX = nodesIndexExpr))
        
        return(newRunCode)	
    }
)

getLogProb_keywordInfo <- keywordInfoClass(
    keyword = 'getLogProb',
    processor = function(code, nfProc, RCfunProc){
        if(!isCodeArgBlank(code, 'nodeFxnVector'))
            return(code)
        errorContext <- deparse(code)
        nodeFunVec_ArgList <- list(model = code$model, nodes = code$nodes, includeData = TRUE, sortUnique = TRUE, errorContext = errorContext)
        if(!isCodeArgBlank(code, 'nodeFunctionIndex')) { ## new case: calculate(myNodeFunctionVector, nodeFunctionIndex = i), if myNodeFunctionVector was hand-created in setup code
            if(!isCodeArgBlank(code, 'nodes'))
                stop('nodes argument cannot be provided to getLogProb if nodeFunctionIndex is specified')
            return(code) ## no modification needed!
        }
        if(isCodeArgBlank(code, 'model'))
            stop('model argument missing from getLogProb, with no accessor argument supplied')
        if(isCodeArgBlank(code, 'nodes')){
            LHSnodes_ArgList <- list(model = code$model)
            LHSnodes_name <- allLHSNodes_SetupTemplate$makeName(LHSnodes_ArgList)
            addNecessarySetupCode(LHSnodes_name, LHSnodes_ArgList, allLHSNodes_SetupTemplate, nfProc, allowToCpp = FALSE)
            nodeFunVec_ArgList$nodes = as.name(LHSnodes_name)
        }
        useNodeFunctionVectorByIndex <- FALSE
        if(hasBracket(nodeFunVec_ArgList$nodes)) { ## like calculate(model, nodes[i]), which could have started as model$calculate(nodes[i])
            useNodeFunctionVectorByIndex <- TRUE
            if(length(nodeFunVec_ArgList$nodes) != 3) stop(paste0('Problem with ', deparse(code),'. If you need to index on the nodes argument there should be only one index.'))
            nodesIndexExpr <- nodeFunVec_ArgList$nodes[[3]]
            nodeFunVec_ArgList$nodes <- nodeFunVec_ArgList$nodes[[2]]
            nodeFunVec_ArgList$sortUnique <- FALSE
        }

        nodeFunName <- nodeFunctionVector_SetupTemplate$makeName(nodeFunVec_ArgList)	
        addNecessarySetupCode(nodeFunName, nodeFunVec_ArgList, nodeFunctionVector_SetupTemplate, nfProc)
        if(!useNodeFunctionVectorByIndex)
            newRunCode <- substitute(getLogProb(nodeFxnVector = NODEFUNVEC_NAME),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName)))
        else
            newRunCode <- substitute(getLogProb(nodeFxnVector = NODEFUNVEC_NAME, nodeFunctionIndex = NODEFUNVECINDEX),
                                     list(NODEFUNVEC_NAME = as.name(nodeFunName), NODEFUNVECINDEX = nodesIndexExpr))
        return(newRunCode)	
    }
)    

nimCopy_keywordInfo <- keywordInfoClass(
	keyword = 'nimCopy',
    processor = function(code, nfProc, RCfunProc){
        if(is.null(nfProc)) stop("Can\'t call copy (nimCopy) from a nimbleFunction without setup code")
		possibleObjects <- c('symbolModel', 'symbolModelValues', 'symbolModelVariableAccessorVector', 'symbolModelValuesAccessorVector')
		modelValuesTypes <- c('symbolModelValues', 'symbolModelValuesAccessorVector')
		accessTypes <- c('symbolModelVariableAccessorVector', 'symbolModelValuesAccessorVector')
		from_ArgList <- list(name = code$from, class = symTypeFromSymTab(code$from, nfProc$setupSymTab, options = possibleObjects))
		to_ArgList <- list(name = code$to, class = symTypeFromSymTab(code$to, nfProc$setupSymTab, options = possibleObjects))
		if(is.null(from_ArgList$class)) 
			stop("Error in nimCopy: '", code$from, "' is not a recognized model or modelValues object.") 		
		if(is.null(to_ArgList$class)) 
			stop("Error in nimCopy: '", code$to, "' is not a recognized model or modelValues object.") 		
		if(from_ArgList$class %in% modelValuesTypes){
			if(isCodeArgBlank(code, 'row'))		stop('row argument missing in copy call')
			from_ArgList$row = code$row
		}
		if(to_ArgList$class %in% modelValuesTypes){
			if(isCodeArgBlank(code, 'rowTo')){
				if(isCodeArgBlank(code, 'row'))		stop('row argument missing in copy call')
				else								to_ArgList$row = code$row
			}
			else		to_ArgList$row = code$rowTo
		}
		if(isCodeArgBlank(code, 'nodes')){
			if(from_ArgList$class == 'symbolModel'){
				node_ArgList <- list(model = from_ArgList$name)
				allNodes_name <- allModelNodes_SetupTemplate$makeName( node_ArgList )
				addNecessarySetupCode(allNodes_name, node_ArgList, allModelNodes_SetupTemplate, nfProc, allowToCpp = FALSE)
			}
			else if(from_ArgList$class == 'symbolModelValues'){
				from_ArgList$row = code$row
				mvVar_ArgList <- list(modelValues = from_ArgList$name)
				allNodes_name <- allModelValuesVars_SetupTemplate$makeName(mvVar_ArgList)
				addNecessarySetupCode(allNodes_name, mvVar_ArgList, allModelValuesVars_SetupTemplate, nfProc, allowToCpp = FALSE)
			}
			from_ArgList$nodes <- as.name(allNodes_name)
		}
		else	from_ArgList$nodes <- code$nodes
		
		if(isCodeArgBlank(code, 'nodesTo'))		to_ArgList$nodes <- from_ArgList$nodes
		else									to_ArgList$nodes <- code$nodesTo
				
		if(from_ArgList$class == 'symbolModel'){
                    isMVfrom <- 0 
                    accessFrom_ArgList <- list(model = code$from, nodes = from_ArgList$nodes, logProb = code$logProb, logProbOnly = code$logProbOnly)
                    accessFrom_name <- modelVariableAccessorVector_setupCodeTemplate$makeName(accessFrom_ArgList)
                    addNecessarySetupCode(accessFrom_name, accessFrom_ArgList, modelVariableAccessorVector_setupCodeTemplate, nfProc)
		}
		else if(from_ArgList$class == 'symbolModelValues'){
                    isMVfrom <- 1
                    accessFrom_ArgList <- list(modelValues = code$from, nodes = from_ArgList$nodes, logProb = code$logProb, logProbOnly = code$logProbOnly, row = from_ArgList$row)
                    accessFrom_name <- modelValuesAccessorVector_setupCodeTemplate$makeName(accessFrom_ArgList)
                    addNecessarySetupCode(accessFrom_name, accessFrom_ArgList, modelValuesAccessorVector_setupCodeTemplate, nfProc)
		}
		else if(from_ArgList$class %in% accessTypes) {
                    isMVfrom <- as.integer(from_ArgList$class == 'symbolModelValuesAccessorVector') 
                    accessFrom_name <- as.character(code$from)
                }
        
		if(to_ArgList$class == 'symbolModel'){
                    isMVto <- 0
                    accessTo_ArgList <- list(model = code$to, nodes = to_ArgList$nodes, logProb = code$logProb, logProbOnly = code$logProbOnly)
                    accessTo_name <- modelVariableAccessorVector_setupCodeTemplate$makeName(accessTo_ArgList)
                    addNecessarySetupCode(accessTo_name, accessTo_ArgList, modelVariableAccessorVector_setupCodeTemplate, nfProc)
		}
		else if(to_ArgList$class == 'symbolModelValues'){
                    isMVto <- 1
                    accessTo_ArgList <- list(modelValues = code$to, nodes = to_ArgList$nodes, logProb = code$logProb, logProbOnly = code$logProbOnly, row = to_ArgList$row)
                    accessTo_name <- modelValuesAccessorVector_setupCodeTemplate$makeName(accessTo_ArgList)
                    addNecessarySetupCode(accessTo_name, accessTo_ArgList, modelValuesAccessorVector_setupCodeTemplate, nfProc)
		}
		else if(to_ArgList$class %in% accessTypes) {
                    isMVto <- as.integer(to_ArgList$class == 'symbolModelValuesAccessorVector') 
                    accessTo_name <- as.character(code$to) 
                }
        if(getNimbleOption('useNewNimCopy')) {
            copierVector_ArgList <- list(accessFrom_name = accessFrom_name, accessTo_name = accessTo_name, isMVto = isMVto, isMVfrom = isMVfrom)
            copierVector_name <- copierVector_setupCodeTemplate$makeName(copierVector_ArgList)
            addNecessarySetupCode(copierVector_name, copierVector_ArgList, copierVector_setupCodeTemplate, nfProc) 
        }
        
        if(!getNimbleOption('useNewNimCopy')) {
            ##What happens below is a bit convoluted and really for backwards compatibility 	
            runCode <- substitute(nimCopy(from = FROM_ACCESS, rowFrom = NA, to = TO_ACCESS, rowTo = NA), 
                                  list(FROM_ACCESS = as.name(accessFrom_name), TO_ACCESS = as.name(accessTo_name)))
            if(from_ArgList$class %in% modelValuesTypes)
                runCode$rowFrom = from_ArgList$row
            if(to_ArgList$class %in% modelValuesTypes)
                runCode$rowTo = to_ArgList$row
        } else {
            rowFromArg <- if(from_ArgList$class %in% modelValuesTypes) from_ArgList$row else NA
            rowToArg <- if(to_ArgList$class %in% modelValuesTypes) {
                if(identical(rowFromArg, NA)) {rowFromArg <- 0; unusedArg <- NA} else unusedArg <- 0
                to_ArgList$row
            } else {
                unusedArg <- NA
                NA
            }
            runCode <- substitute(nimCopy(copierVector = COPIER_VECTOR, rowFrom = ROWFROM, rowTo = ROWTO, unused = UNUSED), 
                                  list(COPIER_VECTOR = as.name(copierVector_name),
                                       ROWFROM = rowFromArg, ROWTO = rowToArg, UNUSED  = unusedArg))
        }
        runCode <- runCode[as.character(runCode) != 'NA']
        return(runCode)
    })


doubleBracket_keywordInfo <- keywordInfoClass(
	keyword = '[[', 
    processor = function(code, nfProc, RCfunProc){
        callerCode <- code[[2]]
        if(is.null(nfProc)) stop("No allowed use of [[ in a nimbleFunction without setup code.")
        possibleObjects <- c('symbolModel', 'symbolNimPtrList', 'symbolNimbleFunctionList', 'symbolNimbleList')
        class = symTypeFromSymTab(code[[2]], nfProc$setupSymTab, options = possibleObjects)
        if(is.null(class)){  ##assume that an element of a run-time provided nimbleList is being accessed
          nl_charName <- as.character(callerCode)
          nl_fieldName <-as.character(code[[3]])
          newRunCode <- substitute(nfVar(NIMBLELIST, VARNAME), list(NIMBLELIST = as.name(nl_charName), VARNAME = nl_fieldName))
          return(newRunCode)
        }
        if(class == 'symbolNimPtrList' || class == 'symbolNimbleFunctionList')
            return(code)
        if(class == 'symbolNimbleList'){
          #	Code is of the form 
          #  myNimbleList[['myVar']]
          nl_charName <- as.character(callerCode)
          nl_fieldName <-as.character(code[[3]])
          newRunCode <- substitute(nfVar(NIMBLELIST, VARNAME), list(NIMBLELIST = as.name(nl_charName), VARNAME = nl_fieldName))
          return(newRunCode)
        }
        if(class == 'symbolModel'){
            singleAccess_ArgList <- list(code = code, model = code[[2]], nodeExpr = code[[3]])
            nodeArg <- code[[3]]
            if(is.character(nodeArg)){
                if(length(nodeArg) > 1)
                    stop(paste0("Problem in ",
                                deparse(code), ". ",
                                deparse(code[[3]]),
                                " is too long.  It can only have one element."),
                         call. = FALSE)
                varAndIndices <- nimbleInternalFunctions$getVarAndIndices(nodeArg)
                
                allNDims <- lapply(nfProc$instances, function(x) {
                    model <- eval(singleAccess_ArgList$model,
                                  envir = x)
                    if(length(nodeArg) != 1)
                        stop(paste0("Length of ",
                                    nodeArg,
                                    " requested from ",
                                    deparse(singleAccess_ArgList$model),
                                    " using '[[' is ",
                                    length(nodeArg),
                                    ". It must be 1.")
                           , call. = FALSE)
                    determineNdimFromOneCase(model, varAndIndices)
                } )

                if(length(unique(allNDims)) > 1)
                    stop(paste0('Error for ', deparse(code),
                                '. Inconsistent numbers of dimensions for different instances.'),
                         call. = FALSE)
                nDim <- allNDims[[1]]
                useMap <- nDim > 0
                
                ## ##
                ## ## If input is of the form model[['a']]
                ## ## and a is non-scalar,
                ## ## we treat it like model$a, which handles either
                if(useMap & length(varAndIndices$indices) == 0) {
                    return(
                        keywordList[['$']]$processor(code, nfProc, RCfunProc)
                    )
                }
                ## ## Following line adds up 0 for each scalar index
                ## ## and 1 for each non-scalar index to determine if the
                ## ## accessed node is scalar
                ## nDim <- sum(1 - unlist(lapply(varAndIndices$indices, is.numeric) ) )
                ## useMap <- nDim > 0
            }
            else{
                allNDims <- determineNdimsFromNfproc(singleAccess_ArgList$model, nodeArg, nfProc)
                if(length(unique(allNDims)) > 1) stop(paste0('Error for ', deparse(code), '. Inconsistent numbers of dimensions for different instances.'), call. = FALSE)
                nDim <- allNDims[[1]]
                useMap <- nDim > 0
            }
            if(useMap){
                accessName <- map_SetupTemplate$makeName(singleAccess_ArgList)
                addNecessarySetupCode(accessName, singleAccess_ArgList, map_SetupTemplate, nfProc)	
                ans <- makeMapAccessExpr(accessName, as.name(accessName), nDim)
            }
            else{
                accessName <- singleModelIndexAccess_SetupTemplate$makeName(singleAccess_ArgList)
                addNecessarySetupCode(accessName, singleAccess_ArgList, singleModelIndexAccess_SetupTemplate, nfProc)
                                        #ans <- substitute(ACCESSNAME[MFLATINDEX], list(ACCESSNAME = as.name(accessName), MFLATINDEX = as.name(paste0(accessName, '_flatIndex'))))
                ans <- makeSingleIndexAccessExpr(accessName, as.name(accessName))
            }
            return(ans)			
        }
        stop("Incorrect use of double brackets in: '", deparse(code), "'.")
    })

modelMemberFun_keywordInfo <- keywordInfoClass(
    keyword = 'multiple',
    processor = function(code, nfProc, RCfunProc) {
        ## if we get here it must be model$member(args)
        ## We will turn it into member(model, args)
        newRunCode <- do.call("call",
                              c(list(as.character(code[[1]][[3]]),
                                     code[[1]][[2]]),
                                as.list(code[-1])),
                              quote = TRUE)
        return(newRunCode)
    })

dollarSign_keywordInfo <- keywordInfoClass(
    keyword = '$',
    processor = function(code, nfProc, RCfunProc){
        callerCode <- code[[2]]
        
        if(is.null(nfProc)) { 

            nl_fieldName <-as.character(code[[3]])
            newRunCode <- substitute(nfVar(NIMBLELIST, VARNAME), list(NIMBLELIST = callerCode, VARNAME = nl_fieldName))
            return(newRunCode)
        }
        
        doubleBracketCase <- FALSE
        if(length(callerCode) > 1) {
            if(deparse(callerCode[[1]] == '[[')) {
                doubleBracketCase <- TRUE
                symObj <- getSymObj_recurse(callerCode[[2]], nfProc$setupSymTab)
            }
        }
        if(!doubleBracketCase)
            symObj <- getSymObj_recurse(callerCode, nfProc$setupSymTab)

        class <- class(symObj)[1] ## symObj is allowed to be NULL       

                                        #	This extracts myNimbleFunction from the expression myNimbleFunction$foo()
        if(length(callerCode) > 1){
            if(callerCode[[1]] == '$'){ ## nested NL or NF case
                callerCode <- processKeyword(callerCode, nfProc, RCfunProc)
            }
        }
                                        #       This extracts myNimbleFunctionList from the expression myNimbleFunctionList[[i]]
                                        #       May be a better way to do this
        

        if(is.null(class) || class == 'NULL'){  ##assume that an element of a run-time provided nimbleList is being accessed
            nl_fieldName <-as.character(code[[3]])
            newRunCode <- substitute(nfVar(NIMBLELIST, VARNAME), list(NIMBLELIST = callerCode, VARNAME = nl_fieldName))
            return(newRunCode)				
        }
        if(class == 'symbolNimbleFunctionList'){
             nf_fieldName <-as.character(code[[3]])
            newRunCode <- substitute(nfMethod(NIMBLEFXN, METHODNAME), list(NIMBLEFXN = callerCode, METHODNAME = nf_fieldName))
            return(newRunCode)
        }
        if(class == 'symbolModel'){
            singleAccess_ArgList <- list(code = code, model = callerCode, var = as.character(code[[3]]) )
            accessName <- singleVarAccess_SetupTemplate$makeName(singleAccess_ArgList)
            addNecessarySetupCode(accessName, singleAccess_ArgList, singleVarAccess_SetupTemplate, nfProc)
            return(as.name(accessName))
        }
        if(class == 'symbolNimbleFunction'){
            
                                        #	Code is of the form myNimbleFunction$myMethod
                                        #   or myNimbleFunction$myVar
            
            
                                        #	Note that we have cut off '()' in the case of myMethod, so we must inspect the
                                        #   nested symbol for myMethod to determine whether it is a method or variable
            
            nf_fieldName <-as.character(code[[3]])
            objectSymbol = symObj$nfProc$setupSymTab$getSymbolObject(nf_fieldName)
            if(class(objectSymbol)[[1]] == 'symbolMemberFunction'){
                newRunCode <- substitute(nfMethod(NIMBLEFXN, METHODNAME), list(NIMBLEFXN = callerCode, METHODNAME = nf_fieldName))				
                return(newRunCode)
            }
            else {
                                        # We *assume* that if its not a member function, it should be treated with 
                                        # nfVar
                newRunCode <- substitute(nfVar(NIMBLEFXN, VARNAME), list(NIMBLEFXN = callerCode, VARNAME = nf_fieldName))
                return(newRunCode)
            }
        }
        if(class == 'symbolNimbleList'){
                                        #	Code is of the form 
                                        #  myNimbleList$myVar
            nl_fieldName <-as.character(code[[3]])
            newRunCode <- substitute(nfVar(NIMBLELIST, VARNAME), list(NIMBLELIST = callerCode, VARNAME = nl_fieldName))
            return(newRunCode)
        }
        if(class == 'symbolNimbleFunctionList'){
            
                                        #	Code is of the form myNimbleFunctionList[[i]]$foo	(foo should be a method)
                                        #	At this point, we cannot access variables of a nimble function list, ie
                                        #	myNimbleFunctionList[[i]]$myVariable is not allowed
                                        #	If we add this functionality, we will need to look up what foo as we do
                                        #	for the nimbleFunction case above
            
            nf_name <-code[[2]]
            nf_fieldName <- as.character(code[[3]])
            newRunCode <- substitute(nfMethod(NIMBLEFXN, METHODNAME), list(NIMBLEFXN = nf_name, METHODNAME = nf_fieldName))
            return(newRunCode)				
        }
    }
)
    
singleBracket_keywordInfo <- keywordInfoClass(
	keyword = '[',
    processor = function(code, nfProc, RCfunProc){
        if(is.null(nfProc)) return (code)
        class <- symTypeFromSymTab(code[[2]], nfProc$setupSymTab)
        if(class == 'symbolModelValues'){
            if(length(code) < 4) {
                stop(paste0("incorrect syntax for accessing element of modelValues: ", deparse(code)))
            }
            singleMVAccess_ArgList <- list(code = code, modelValues = code[[2]], var = code[[3]], row = code[[4]])
            accessName <- singleModelValuesAccessor_SetupTemplate$makeName(singleMVAccess_ArgList)
            addNecessarySetupCode(accessName, singleMVAccess_ArgList, singleModelValuesAccessor_SetupTemplate, nfProc)
            if(length(code) == 4)
                indexExpr = code[[4]]
            else
                indexExpr = substitute(1)
            
            return(substitute(ACCESS[INDEX], list(ACCESS = as.name(accessName), INDEX = indexExpr) ) )
        }
        return(code)
    }
)    

length_char_keywordInfo <- keywordInfoClass(
    keyword = 'length',
    processor = function(code, nfProc, RCfunProc) {
        if(is.null(nfProc)) return(code)
        if(length(code) < 2) stop('length() used without an argument')
        class <- symTypeFromSymTab(code[[2]], nfProc$setupSymTab)
        if(class == "symbolString") {
            length_char_ArgList <- list(code = code, string = code[[2]])
            accessName <- length_char_SetupTemplate$makeName(length_char_ArgList)
            addNecessarySetupCode(accessName, length_char_ArgList, length_char_SetupTemplate, nfProc)
            return(substitute(LENGTHNAME, list(LENGTHNAME = as.name(accessName))))
        }
        return(code)
    })

nimIntegrate_keywordInfo <- keywordInfoClass(
    keyword = 'nimIntegrate',
    processor = function(code, nfProc, RCfunProc) {
        if(code$abs.tol == quote(rel.tol)) 
            code$abs.tol = code$rel.tol
        iTols <- which(names(code) %in% c('rel.tol','abs.tol'))
        for(i in iTols) {
            code_i <- code[[i]]
            if(length(code_i) > 1) {
                if(".Machine" %in% all.vars(code_i)) # Must be an expression using .Machine$double.eps or related
                    code[[i]] <- eval(code_i)
            }
        }
        code
    }
)

nimDerivs_keywordInfo <- keywordInfoClass(
  keyword = 'nimDerivs',
  processor = function(code, nfProc, RCfunProc) {
    wrtArgs <- code[['wrt']]
    ## First check to see if nimFxn argument is a method.
    fxnCall <- code[[2]][[1]]
    order <- code[['order']]

    calculateCase <- FALSE
    if(deparse(fxnCall) == 'calculate'){
        calculateCase <- TRUE
    } 
    else if(length(fxnCall) == 3 &&
            (deparse(fxnCall[[1]]) == '$' &&
             deparse(fxnCall[[3]]) == 'calculate')){
        ## re-arrange nimDerivs(model$calculate(...), ...) to nimDerivs(calculate(model, ...), ...)
        code[[2]] <- modelMemberFun_keywordInfo$processor(code[[2]], nfProc)
        code[[2]] <- matchKeywordCode(code[[2]], nfProc)
        calculateCase <- TRUE
    }
    if(calculateCase) {
        innerCode <- code[[2]]
        if(length(code$wrt)==1)
          if(!is.name(code$wrt))
            if(is.na(code$wrt[1]))
              stop("Derivatives of a call to 'calculate()' must have 'wrt' argument specified.")
        innerCode$wrt <- code$wrt
        newCode <- calculate_keywordInfo$processor(innerCode, nfProc, RCfunProc)
        newCode[[1]] <- as.name('nimDerivs_calculate')
        newCode$orderVector <- code$order
        newCode$update <- if(is.null(code[['update']])) TRUE else code[['update']]
        newCode$reset <- if(is.null(code[['reset']])) FALSE else code[['reset']]
        return(newCode)
    }
    ## Not a calculate case:
    if(!is.null(nfProc$origMethods[[deparse(fxnCall)]])) {
      derivMethod <- nfProc$origMethods[[deparse(fxnCall)]]
      derivMethodArgs <- derivMethod$getArgInfo()

      ## Only make the wrt substitution if the names are baked in as character
      wrtArg <- code$wrt
      doPreprocess <- FALSE
      if(is.numeric(wrtArg) | is.logical(wrtArg)) {
        if(any(is.na(wrtArg[1]))) { ## wrt = NULL (default), set to NA, which will become -1 by convertWrtArgToIndices and then c(-1, -1) in the setup code
          doPreprocess <- TRUE
        }
      } else if(is.name(wrtArg)) { ## wrt = some_index_var_maybe_from_setup
        ## do nothing
      } else if(deparse(wrtArg[[1]]) == 'nimC') { ## wrt = c('x', ...)
        if(is.character(wrtArg[[2]]))
          doPreprocess <- TRUE
      } else if(is.character(wrtArg)) { ## wrt = 'x'
        doPreprocess <- TRUE
      } ## In any other setting, it must be an index vector and we leave it alone here.
      ##
      if(doPreprocess) {
        wrt_argList <- list(fxnCall = fxnCall,
                            wrt = code$wrt,
                            #                   vector = wrtArgIndices
                            derivMethodArgs = derivMethodArgs)
        accessName <- wrtVector_setupCodeTemplate$makeName(wrt_argList)
        addNecessarySetupCode(accessName, wrt_argList, wrtVector_setupCodeTemplate, nfProc)
        code[['wrt']] <- substitute(VECNAME,
                                    list(VECNAME = as.name(accessName)))
      }
      ## Check if model and updateNodes are provided.
      ## If so, create a nodeFxnVector for them.
      modelProvided <-  is.name(code[['model']])
      updateNodesProvided <- is.name(code[['updateNodes']])
      if(!updateNodesProvided)
        if(is.character(code[['updateNodes']]))
          updateNodesProvided <- !is.na(code[['updateNodes']][1])
      constantNodesProvided <- is.name(code[['constantNodes']])
      if(!constantNodesProvided)
        if(is.character(code[['constantNodes']]))

          constantNodesProvided <- !is.na(code[['constantNodes']][1])
      if(xor(modelProvided, updateNodesProvided || constantNodesProvided))
        stop('nimDerivs arguments model and at least one of updateNodes or constantNodes must be provided together or not at all.')
      if(modelProvided) { ## At this point that means both were provided
        DMUNargList <- list(model = code$model,
                            updateNodes = if(updateNodesProvided) code[["updateNodes"]] else NULL,
                            constantNodes = if(constantNodesProvided) code[["constantNodes"]] else NULL)
        accessName <- nodeFunctionVector_DerivsModelUpdateNodes_SetupTemplate$makeName(DMUNargList)
        addNecessarySetupCode(accessName, DMUNargList, nodeFunctionVector_DerivsModelUpdateNodes_SetupTemplate, nfProc)
        ## code[['model']] <- NULL
        ## if(updateNodesProvided) code[['updateNodes']] <- NULL
        ## if(constantNodesProvided) code[['constantNodes']] <- NULL
        code[['updateNodesName']] <- accessName
      }
      code[['model']] <- NULL
      code[['updateNodes']] <- NULL
      code[['constantNodes']] <- NULL
    }
    return(code)
  }
)

derivInfo_keywordInfo <- keywordInfoClass(
  keyword = 'derivInfo',
  processor = function(code, nfProc, RCfunProc) {
    fxnCall <- code[[2]][[1]]
    calculateCase <- FALSE
    if(deparse(fxnCall) == 'calculate'){
      calculateCase <- TRUE
    } 
    else if(length(fxnCall) == 3 &&
              (deparse(fxnCall[[1]]) == '$' &&
                 deparse(fxnCall[[3]]) == 'calculate')){
      ## re-arrange nimDerivs(model$calculate(...), ...) to nimDerivs(calculate(model, ...), ...)
      code[[2]] <- modelMemberFun_keywordInfo$processor(code[[2]], nfProc)
      code[[2]] <- matchKeywordCode(code[[2]], nfProc)
      calculateCase <- TRUE
    }
    if(!calculateCase) {
      stop('The first argument to derivInfo must be a call to calculate(model, ...) or model$calculate(...)')
    }
    innerCode <- code[[2]]
    if(is.null(code$wrt))
      stop("derivInfo must provide 'wrt' argument.")
    innerCode$wrt <- code$wrt
    newCode <- calculate_keywordInfo$processor(innerCode, nfProc, RCfunProc)
    return(newCode)
  }
)

#	KeywordList
keywordList <- new.env()
keywordList[['nimSeq']] <- nimSeq_keywordInfo
keywordList[['getParam']] <- getParam_keywordInfo
keywordList[['getBound']] <- getBound_keywordInfo
keywordList[['values']] <- values_keywordInfo
keywordList[['calculate']] <- calculate_keywordInfo
keywordList[['calculateDiff']] <- calculateDiff_keywordInfo
keywordList[['simulate']] <- simulate_keywordInfo
keywordList[['getLogProb']] <- getLogProb_keywordInfo
keywordList[['nimCopy']] <- nimCopy_keywordInfo
keywordList[['[[']] <- doubleBracket_keywordInfo
keywordList[['$']] <- dollarSign_keywordInfo
keywordList[['[']] <- singleBracket_keywordInfo
keywordList[['besselK']] <- besselK_keywordInfo
keywordList[['dgamma']] <- d_gamma_keywordInfo
keywordList[['pgamma']] <- pq_gamma_keywordInfo
keywordList[['qgamma']] <- pq_gamma_keywordInfo
keywordList[['rgamma']] <- rgamma_keywordInfo
keywordList[['dinvgamma']] <- d_gamma_keywordInfo
keywordList[['pinvgamma']] <- pq_gamma_keywordInfo
keywordList[['qinvgamma']] <- pq_gamma_keywordInfo
keywordList[['rinvgamma']] <- rgamma_keywordInfo
keywordList[['dsqrtinvgamma']] <- d_gamma_keywordInfo
keywordList[['psqrtinvgamma']] <- pq_gamma_keywordInfo
keywordList[['qsqrtinvgamma']] <- pq_gamma_keywordInfo
keywordList[['rsqrtinvgamma']] <- rgamma_keywordInfo
keywordList[['ddexp']] <- d_gamma_keywordInfo
keywordList[['pdexp']] <- pq_gamma_keywordInfo
keywordList[['qdexp']] <- pq_gamma_keywordInfo
keywordList[['rdexp']] <- rgamma_keywordInfo
# can be handled the same as gamma, so include although we have dexp_nimble too
keywordList[['dexp']] <- d_gamma_keywordInfo
keywordList[['pexp']] <- pq_gamma_keywordInfo
keywordList[['qexp']] <- pq_gamma_keywordInfo
keywordList[['rexp']] <- rgamma_keywordInfo

keywordList[['dexp_nimble']] <- d_exp_nimble_keywordInfo
keywordList[['pexp_nimble']] <- pq_exp_nimble_keywordInfo
keywordList[['qexp_nimble']] <- pq_exp_nimble_keywordInfo
keywordList[['rexp_nimble']] <- rexp_nimble_keywordInfo

keywordList[['length']] <- length_char_keywordInfo ## active only if argument has type character

keywordList[['nimIntegrate']] <- nimIntegrate_keywordInfo

keywordList[['nimDerivs']] <- nimDerivs_keywordInfo
keywordList[['derivInfo']] <- derivInfo_keywordInfo

keywordListModelMemberFuns <- new.env()
keywordListModelMemberFuns[['calculate']] <- modelMemberFun_keywordInfo
keywordListModelMemberFuns[['simulate']] <- modelMemberFun_keywordInfo
keywordListModelMemberFuns[['calculateDiff']] <- modelMemberFun_keywordInfo
keywordListModelMemberFuns[['getLogProb']] <- modelMemberFun_keywordInfo
keywordListModelMemberFuns[['getParam']] <- modelMemberFun_keywordInfo
keywordListModelMemberFuns[['getBound']] <- modelMemberFun_keywordInfo


matchFunctions <- new.env()
matchFunctions[['setSize']] <- function(var, ..., copy = TRUE, fillZeros = TRUE){} 
matchFunctions[['nimC']] <- nimC
matchFunctions[['nimRep']] <- function(x, times = 1, length.out, each = 1) {}
matchFunctions[['nimSeq']] <- nimSeq
matchFunctions[['nimNumeric']] <- nimNumeric 
matchFunctions[['nimInteger']] <- nimInteger 
matchFunctions[['nimLogical']] <- nimLogical 
matchFunctions[['nimMatrix']] <- nimMatrix 
matchFunctions[['nimArray']] <- nimArray 
matchFunctions[['values']] <- function(model, nodes, accessor){}
matchFunctions[['getParam']] <- getParam
matchFunctions[['getBound']] <- getBound
matchFunctions[['calculate']] <- calculate
matchFunctions[['calculateDiff']] <- calculateDiff
matchFunctions[['simulate']] <- simulate
matchFunctions[['getLogProb']] <- getLogProb
matchFunctions[['nimCopy']] <- function(from, to, nodes, nodesTo, row, rowTo, logProb = FALSE, logProbOnly = FALSE){}
matchFunctions[['double']] <- function(nDim, dim, default, ...){}
matchFunctions[['int']] <- function(nDim, dim, default, ...){}
matchFunctions[['nimOptim']] <- nimOptim
matchFunctions[['nimIntegrate']] <- nimIntegrate
matchFunctions[['nimOptimDefaultControl']] <- nimOptimDefaultControl
matchFunctions[['nimEigen']] <- function(squareMat, symmetric = FALSE, only.values = FALSE){}
matchFunctions[['nimSvd']] <- function(mat, vectors = 'full'){}
matchFunctions[['nimDerivs']] <- function(call = NA,
                                          order = nimC(0,1,2),
                                          wrt = NA,
                                          model = NA, updateNodes = NA, constantNodes = NA,
                                          do_update = TRUE, reset = FALSE) {}
matchFunctions[['derivInfo']] <- derivInfo
matchFunctions[['besselK']] <- function(x, nu, expon.scaled = FALSE){}
matchFunctions[['dgamma']] <- function(x, shape, rate = 1, scale, log = FALSE){}
matchFunctions[['rgamma']] <- function(n, shape, rate = 1, scale){}
matchFunctions[['qgamma']] <- function(p, shape, rate = 1, scale, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['pgamma']] <- function(q, shape, rate = 1, scale, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['dinvgamma']] <- matchFunctions[['dsqrtinvgamma']] <- function(x, shape, scale = 1, rate, log = FALSE){}
matchFunctions[['rinvgamma']] <- matchFunctions[['rsqrtinvgamma']] <- function(n, shape, scale = 1, rate){}
matchFunctions[['qinvgamma']] <- matchFunctions[['qsqrtinvgamma']] <- function(p, shape, scale = 1, rate, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['pinvgamma']] <- matchFunctions[['psqrtinvgamma']] <- function(q, shape, scale = 1, rate, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['dexp']] <- function(x, rate = 1, log = FALSE){}
matchFunctions[['rexp']] <- function(n, rate = 1){}
matchFunctions[['qexp']] <- function(p, rate = 1, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['pexp']] <- function(q, rate = 1, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['dexp_nimble']] <- function(x, rate, scale = 1, log = FALSE){}
matchFunctions[['rexp_nimble']] <- function(n, rate, scale = 1){}
matchFunctions[['qexp_nimble']] <- function(p, rate, scale = 1, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['pexp_nimble']] <- function(q, rate, scale = 1, lower.tail = TRUE, log.p = FALSE){}

matchFunctions[['ddexp']] <- function(x, location, scale = 1, rate, log = FALSE){}
matchFunctions[['rdexp']] <- function(n, location, scale = 1, rate){}
matchFunctions[['qdexp']] <- function(p, location, scale = 1, rate, lower.tail = TRUE, log.p = FALSE){}
matchFunctions[['pdexp']] <- function(q, location, scale = 1, rate, lower.tail = TRUE, log.p = FALSE){}

matchModelMemberFunctions <- new.env()
matchModelMemberFunctions[['calculate']] <- function(nodes) {}
matchModelMemberFunctions[['calculateDiff']] <- function(nodes) {}
matchModelMemberFunctions[['getLogProb']] <- function(nodes) {}
matchModelMemberFunctions[['simulate']] <- function(nodes, includeData = FALSE) {}
matchModelMemberFunctions[['getParam']] <- function(node, param) {}
matchModelMemberFunctions[['getBound']] <- function(node, bound) {}

# remove ncp from signatures
stripArgs <- function(fname, argNames) {
    if(exists(fname)) {
        args <- formals(eval(as.name(fname)))
        args <- args[-which(names(args) %in% argNames)]
        template <- function() {}
        formals(template) <- args
        return(template)
    } else return(NULL)
}

for(distfun in paste0(c('d','p','q','r'), 'beta'))
    matchFunctions[[distfun]] <- stripArgs(distfun, 'ncp')
for(distfun in paste0(c('d','p','q','r'), 'chisq'))
    matchFunctions[[distfun]] <- stripArgs(distfun, 'ncp')
for(distfun in paste0(c('d','p','q','r'), 't'))
    matchFunctions[[distfun]] <- stripArgs(distfun, 'ncp')
for(distfun in paste0(c('d','p','q','r'), 'nbinom'))
    matchFunctions[[distfun]] <- stripArgs(distfun, 'mu')


# the following are standard in terms of both matchFunctions and keywordList
matchDistList <- list('binom', 'cat', 'dirch', 'interval', 'lnorm', 'logis', 'multi', 'mnorm_chol', 'mvt_chol', 'norm', 'pois', 't_nonstandard', 'unif', 'weibull', 'wish_chol', 'invwish_chol', 'car_normal', 'car_proper')

addDistList2matchFunctions <- function(distList, matchFunEnv){
	for(thisDist in distList){
		pFun <- paste0('p', thisDist)
		qFun <- paste0('q', thisDist)
		rFun <- paste0('r', thisDist)
		dFun <- paste0('d', thisDist)
		
                eval(substitute(matchFunctions[[dFun]] <- DFUN, list(DFUN = as.name(dFun))))
                eval(substitute(matchFunctions[[rFun]] <- RFUN, list(RFUN = as.name(rFun))))
                if(exists(qFun))
                    eval(substitute(matchFunctions[[qFun]] <- QFUN, list(QFUN = as.name(qFun))))
                if(exists(pFun))
                    eval(substitute(matchFunctions[[pFun]] <- PFUN, list(PFUN = as.name(pFun))))
	}
}
          

addDistList2matchFunctions(matchDistList, matchFunctions)

#	processKeyword function to be called by nfProc
processKeyword <- function(code, nfProc, RCfunProc){
  thisKeywordInfo <- keywordList[[ as.character(code[[1]]) ]]
  if(!is.null(thisKeywordInfo))
    return(thisKeywordInfo$processor(code, nfProc, RCfunProc))
  return(code)
}

processKeywordCodeMemberFun <- function(code, nfProc, RCfunProc) { ## handle cases like a$b(c) as one unit
    ## this includes nf$method()
    ## nfList[[i]]$method
    ## model$calculate(nodes)
    dollarSignPart <- code[[1]]
    objectPart <- dollarSignPart[[2]]
    
    isModel <- FALSE
    if(length(objectPart) != 1) isModel <- FALSE ## a case like a[[i]]$b(), which can only be a nimbleFunction list
    else {
        symObj <- nfProc$setupSymTab$getSymbolObject(as.character(objectPart))
        if(is.null(symObj)) stop(paste0("In processKeywordCodeMemberFun: not sure what to do with ", deparse(code)))
        if(inherits(symObj, 'symbolModel'))
            isModel <- TRUE
    }
    if(isModel) {
        thisKeywordInfo <- keywordListModelMemberFuns[[ as.character(dollarSignPart[[3]]) ]]
        if(is.null(thisKeywordInfo)) stop(paste0("In processKeywordCodeMemberFun, don't know what do with: ", deparse(code)))
        rearrangedCode <- thisKeywordInfo$processor(code, nfProc)
        rearrangedCode <- matchKeywordCode(rearrangedCode, nfProc)
        return(processKeywords_recurse(rearrangedCode, nfProc, RCfunProc))
    } else {
        ## same as processKeywords_recurse
        ## first line here creates something like nfMethod(model, method)(args)
        ## which is handled as a chainedCall in later processing
        code[[1]] <- processKeywords_recurse(code[[1]], nfProc, RCfunProc)
        cl <- length(code)
        if(cl >= 2) {
            for(i in 2:cl) {
                code[[i]] <- processKeywords_recurse(code[[i]], nfProc, RCfunProc)
            }
        }
        return(code)
    }
}

processKeywords_recurse <- function(code, nfProc = NULL, RCfunProc) {
    cl = length(code)
    if(cl == 1) {
        if(is.call(code)) {
            if(length(code[[1]]) > 1)
                if(deparse(code[[1]][[1]] == '$'))
                    code <- processKeywordCodeMemberFun(code, nfProc, RCfunProc)
                else
                    code[[1]] <- processKeywords_recurse(code[[1]], nfProc, RCfunProc)
        }
        return(code)
    }
    
    if(length(code[[1]]) == 1) {
        code <- processKeyword(code, nfProc, RCfunProc)
    }
    
    cl = length(code)
    
    if(is.call(code)) {
        if(length(code[[1]]) > 1) {
            if(deparse(code[[1]][[1]] == '$')) {
                code <- processKeywordCodeMemberFun(code, nfProc, RCfunProc) ## case like model$calculate(nodes)
                return(code) ## don't recurse on arguments of anything in this category
            }
            code[[1]] <- processKeywords_recurse(code[[1]], nfProc, RCfunProc)
        }
        if(cl >= 2) {
            for(i in 2:cl) {
                code[[i]] <- processKeywords_recurse(code[[i]], nfProc, RCfunProc)
            }
        }
    }
    return(code)
}

#####	SETUPCODE TEMPLATES

wrtVector_setupCodeTemplate <- setupCodeTemplateClass(
    makeName = function(argList){Rname2CppName(paste0('wrtVec_',
                                                      deparse(argList$fxn),
                                                      paste(deparse(argList$wrt), sep='_'),
                                                      '_'))},
  codeTemplate = quote(
  {
      WRTVEC <- nimble:::convertWrtArgToIndices(WRT, #code$wrt
                                                 DERIVMETHODARGS,
                                                 FXNNAME)
      if(length(WRTVEC) == 1)
          WRTVEC <- c(WRTVEC, -1)
  }
  ),
  makeCodeSubList = function(resultName, argList){
      list(WRTVEC = as.name(resultName),
           WRT = argList$wrt
          ,DERIVMETHODARGS = argList$derivMethodArgs
          ,FXNNAME = deparse(argList$fxnCall)
           )##argList$vector)
  }
)

length_char_SetupTemplate <- setupCodeTemplateClass(
    makeName = function(argList) {Rname2CppName(paste0(paste("length", deparse(argList$string), sep='_'), '_KNOWN_'))},
    codeTemplate = quote(LENGTHNAME <- CODE),
    makeCodeSubList = function(resultName, argList) {
        list(LENGTHNAME = as.name(resultName),
             CODE = argList$code)
    })

modelVariableAccessorVector_setupCodeTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are model, nodes and logProb
    makeName = function(argList) {Rname2CppName(paste(deparse(argList$model), deparse(argList$nodes), 'access_logProb', deparse(argList$logProb),'LPO', deparse(argList$logProbOnly), sep = '_'))},
    codeTemplate = quote( ACCESSNAME <- nimble:::modelVariableAccessorVector(MODEL, NODES, logProb = LOGPROB, logProbOnly = LOGPROBONLY) ),
    makeCodeSubList = function(resultName, argList) {
        list(ACCESSNAME = as.name(resultName),
             MODEL = argList$model,
             NODES = argList$nodes,
             LOGPROB = argList$logProb,
             LOGPROBONLY = argList$logProbOnly)
    })

copierVector_setupCodeTemplate <- setupCodeTemplateClass(
    makeName = function(argList) {Rname2CppName(paste0(argList$accessFrom_name, '_', argList$accessTo_name))},
    codeTemplate = quote( COPIERNAME <- nimble:::copierVector(ACCESS_FROM, ACCESS_TO, ISMVFROM, ISMVTO) ),
    makeCodeSubList = function(resultName, argList) {
        list(COPIERNAME = as.name(resultName),
             ACCESS_FROM = as.name(argList$accessFrom_name),
             ACCESS_TO   = as.name(argList$accessTo_name),
             ISMVFROM    = as.integer(argList$isMVfrom),
             ISMVTO      = as.integer(argList$isMVto)) 
    })
    

modelValuesAccessorVector_setupCodeTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are modelValues, nodes and logProb

    makeName = function(argList) {Rname2CppName(paste(deparse(argList[['modelValues']]), deparse(argList$nodes), 'access_logProb', deparse(argList$logProb), 'LPO', deparse(argList$logProbOnly), deparse(argList$row), sep = '_'))},
    codeTemplate = quote( ACCESSNAME <- nimble:::modelValuesAccessorVector(MODEL, NODES, logProb = LOGPROB, logProbOnly = LOGPROBONLY) ),
	makeCodeSubList = function(resultName, argList) {
        list(ACCESSNAME = as.name(resultName),
             MODEL = argList[['modelValues']],
             NODES = argList$nodes,
             LOGPROB = argList$logProb,
             LOGPROBONLY = argList$logProbOnly)
    })


nodeFunctionVector_WithDerivsOutput_SetupTemplate <- setupCodeTemplateClass(
  makeName = function(argList){
    Rname2CppName(paste(deparse(argList$model),
                        deparse(argList$nodes),
                        'nodeFxnVector_WithDerivsOutput', ## The "_derivs_" tag is used later to determine the right class - klugey
                        deparse(argList$includeData),
                        if(argList$sortUnique) "SU" else "notSU",
                        '_derivs_',
                        sep = '_')
                  )
  },
  codeTemplate = quote(
    NODEFXNVECNAME <- nimble:::nodeFunctionVector_WithDerivsOutputNodes(
      model = MODEL,
      calcNodes = CALCNODES,
      excludeData = EXCLUDEDATA,
      sortUnique = SORTUNIQUE)
  ), 
  makeCodeSubList = function(resultName, argList){
    list(NODEFXNVECNAME = as.name(resultName),
         MODEL = argList$model,
         CALCNODES = argList$nodes,
         EXCLUDEDATA = !argList$includeData,
         SORTUNIQUE = argList$sortUnique
         )
  })


nodeFunctionVector_DerivsModelUpdateNodes_SetupTemplate <- setupCodeTemplateClass(
  makeName = function(argList){
    Rname2CppName(paste(deparse(argList$model),
                        deparse(argList[['updateNodes']]),
                        deparse(argList[['constantNodes']]),
                        'nodeFxnVector_DerivsModelUpdateNodes_derivs_', ## The "_derivs_" tag is used later to determine the right class - klugey
                        sep = '_')
                  )
  },
  codeTemplate = quote(
    NODEFXNVECNAME <- nimble:::nodeFunctionVector_DerivsModelUpdateNodes(
      model = MODEL,
      updateNodes = UPDATENODES,
      constantNodes = CONSTANTNODES)
  ), 
  makeCodeSubList = function(resultName, argList){
    list(NODEFXNVECNAME = as.name(resultName),
         MODEL = argList$model,
         UPDATENODES = argList[["updateNodes"]],
         CONSTANTNODES = argList[["constantNodes"]]
         )
  })

nodeFunctionVector_SetupTemplate <- setupCodeTemplateClass(
                                        #Note to programmer: required fields of argList are model, nodes and includeData
    
    makeName = function(argList){
        Rname2CppName(paste(deparse(argList$model),
                            deparse(argList$nodes),
                            'nodeFxnVector_includeData',
                            deparse(argList$includeData),
                            if(argList$sortUnique) "SU" else "notSU",
                            if(is.null(argList$wrtNodes)) '' else '_derivs_', sep = '_')
                      )
    },
    codeTemplate = quote(
        NODEFXNVECNAME <- nimble:::nodeFunctionVector(model = MODEL,
                                                      nodeNames = NODES,
                                                      wrtNodes = WRTNODES,
                                                      excludeData = EXCLUDEDATA,
                                                      sortUnique = SORTUNIQUE,
                                                      errorContext = ERRORCONTEXT)
    ), 
    makeCodeSubList = function(resultName, argList){
        list(NODEFXNVECNAME = as.name(resultName),
             MODEL = argList$model,
             NODES = argList$nodes,
             WRTNODES = argList$wrtNodes,
             EXCLUDEDATA = !argList$includeData,
             SORTUNIQUE = argList$sortUnique,
             ERRORCONTEXT = argList$errorContext
             )
    })


paramInfo_SetupTemplate <- setupCodeTemplateClass(
    #Note to programmer: required fields of argList are model, node and param
    makeName = function(argList){Rname2CppName(paste(deparse(argList$model), deparse(argList$node), deparse(argList$param), 'paramInfo', sep='_'))},
    makeOtherNames = function(name,argList) {Rname2CppName(paste0(name,'_ID'))},
    codeTemplate = quote({
        PARAMINFONAME <- nimble:::makeParamInfo(model = MODEL, nodes = NODE, param = PARAM, vector = HASINDEX )
        PARAMIDNAME <- PARAMINFONAME$paramID
        PARAMINFONAME$paramID <- NULL
       }),
    makeCodeSubList = function(resultName, argList){
        list(PARAMINFONAME = as.name(resultName),
             PARAMIDNAME = as.name(paste0(resultName,'_ID')),
             MODEL = argList$model,
             NODE = argList$node,
             PARAM = argList$param,
             HASINDEX = argList$hasIndex)
    })

boundInfo_SetupTemplate <- setupCodeTemplateClass(
    #Note to programmer: required fields of argList are model, node and param
    makeName = function(argList){Rname2CppName(paste(deparse(argList$model), deparse(argList$node), deparse(argList$bound), 'boundInfo', sep='_'))},
    makeOtherNames = function(name,argList) {Rname2CppName(paste0(name,'_ID'))},
    codeTemplate = quote({
        BOUNDINFONAME <- nimble:::makeBoundInfo(MODEL, NODE, BOUND)
        BOUNDIDNAME <- BOUNDINFONAME$boundID
        BOUNDINFONAME$boundID <- NULL
       }),
    makeCodeSubList = function(resultName, argList){
        list(BOUNDINFONAME = as.name(resultName),
             BOUNDIDNAME = as.name(paste0(resultName,'_ID')),
             MODEL = argList$model,
             NODE = argList$node,
             BOUND = argList$bound)
    })

allLHSNodes_SetupTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are model

	makeName = function(argList){
		Rname2CppName(paste('allLHSnodes', deparse(argList$model), sep = '_'))
	},
	codeTemplate = quote(NODENAMES <- MODEL$getMaps('nodeNamesLHSall')),
	makeCodeSubList = function(resultName, argList){
		list(NODENAMES = as.name(resultName),
			MODEL = argList$model)
	})
	
allModelNodes_SetupTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are model

	makeName = function(argList){
		Rname2CppName(paste('allModelNodes', deparse(argList$model), sep = '_'))
	},
	codeTemplate = quote(NODENAMES <- MODEL$getVarNames()),
	makeCodeSubList = function(resultName, argList){
		list(NODENAMES = as.name(resultName),
			MODEL = argList$model)
	})	
	
allModelValuesVars_SetupTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are modelValues

	makeName = function(argList){
		Rname2CppName(paste('allMVVars', deparse(argList$modelValues), sep = '_'))
	},
	codeTemplate = quote(NODENAMES <- MODELVALUES$getVarNames(includeLogProb = FALSE)),	
		
	makeCodeSubList = function(resultName, argList){
		list(NODENAMES = as.name(resultName),
			MODELVALUES = argList$modelValues)
	})	
	
code2Name_fromArgList <- function(argList)
	Rname2CppName(deparse(argList$code))	
	
singleVarAccess_SetupTemplate <- setupCodeTemplateClass(
	#Note to progammer: required fields of argList are 'code' (raw code to be processed), model and var

	makeName = code2Name_fromArgList,

	codeTemplate = quote(SINGLEACCESSOR <- nimble:::singleVarAccess(MODEL, VAR)),

	makeCodeSubList = function(resultName, argList){
		list(SINGLEACCESSOR = as.name(resultName),
			MODEL = argList$model,
			VAR = argList$var)	
	})	
	
singleModelIndexAccess_SetupTemplate <- setupCodeTemplateClass(
	#Note to progammer: required fields of argList are code, varAndIndices, node (character) and model(expression)
	makeOtherNames = function(name, argList){ paste0(name, '_flatIndex')},
	makeName = code2Name_fromArgList,
	
	codeTemplate = quote({
		VARANDINDICES <- nimble:::nimbleInternalFunctions$getVarAndIndices(NODEVARNAME)
		NEWVARNAME <- as.character(VARANDINDICES$varName)
		MFLATINDEX <- nimble:::varAndIndices2flatIndex(VARANDINDICES, MODELVAREXPR$getVarInfo(NEWVARNAME))
		VARACCESSOR <- nimble:::singleVarAccess(MODELVAREXPR, NEWVARNAME, useSingleIndex = TRUE)
	}),
	makeCodeSubList = function(resultName, argList){
		list(VARACCESSOR = as.name(resultName),
			VARANDINDICES = as.name(paste0(resultName, '_varAndIndices') ),
			NEWVARNAME = as.name(paste0(resultName, '_newVarName')),
			NODEVARNAME = argList$nodeExpr,
			MFLATINDEX = as.name(paste0(resultName, '_flatIndex')),
			MODELVAREXPR = argList$model)
	})
	
map_SetupTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are code, model
	makeName  = code2Name_fromArgList,
	makeOtherNames = function(name, argList){
		output <- character()
		output[1] = paste0(name, '_strides')
		output[2] = paste0(name, '_offset')
		output[3] = paste0(name, '_sizes')
		return(output)
	},
	codeTemplate = quote({
		VARANDINDICES <- nimble:::nimbleInternalFunctions$getVarAndIndices(NODEVARNAME)
		NEWVARNAME <- as.character(VARANDINDICES$varName)
                map_SetupTemplate_vi <- MODEL$getVarInfo(NEWVARNAME)
		map_SetupTemplate_mapParts <- nimble:::varAndIndices2mapParts(VARANDINDICES, map_SetupTemplate_vi$maxs, map_SetupTemplate_vi$nDim)
		MSTRIDES <- map_SetupTemplate_mapParts$strides
		MOFFSET <- map_SetupTemplate_mapParts$offset
		MSIZES <- map_SetupTemplate_mapParts$sizes
		VARACCESSOR <- nimble:::singleVarAccess(model, NEWVARNAME)
	}),
	makeCodeSubList = function(resultName, argList){
		list(VARACCESSOR = as.name(resultName),
		NODEVARNAME =	argList$nodeExpr,
		NEWVARNAME = as.name(paste0(resultName, '_newVarName')),
		VARANDINDICES = as.name(paste0(resultName, '_varAndIndices')),
		MODEL = argList$model,
		MSTRIDES = as.name(paste0(resultName, '_strides')),
		MOFFSET = as.name(paste0(resultName, '_offset')),

                     MSIZES = as.name(paste0(resultName, '_sizes')))
	})
	
singleModelValuesAccessor_SetupTemplate <- setupCodeTemplateClass(
	#Note to programmer: required fields of argList are modelValues, var, row, code
	makeName = code2Name_fromArgList,
	codeTemplate = quote({
		MVACCESS <- nimble:::singleModelValuesAccess(MODELVALUES, VAR)
	}),
	makeCodeSubList = function(resultName, argList){
		list(MVACCESS = as.name(resultName),
		MODELVALUES = argList$modelValues,
		VAR = argList$var)
	}
)





#### KEYWORD PROCESSING UTILITIES


isCodeArgBlank <- function(code, arg){
	#return(nchar(code[[arg]]) == 0)
	return(is.null(code[[arg]]))
}

# Utility functions to make things a little neater
isSetupCodeGenerated <- function(name, nfProc)
    name %in% nfProc$newSetupOutputNames
addSetupCodeNames <- function(name, otherNames, nfProc)
    nfProc$newSetupOutputNames <- c(name, otherNames, nfProc$newSetupOutputNames)
addBlockFromCppName <- function(name, nfProc)
    nfProc$blockFromCppNames <- c(name, nfProc$blockFromCppNames)
addNewCode <- function(name, subList, template, nfProc)
    nfProc$newSetupCode[[name]] <- eval(substitute(substitute(TEMPLATE, subList), list(TEMPLATE = template$codeTemplate) ) )

addNecessarySetupCode <- function(name, argList, template, nfProc, allowToCpp = TRUE){
    if(is.null(nfProc)) stop("Trying to add setup code for a nimbleFunction with no setup code.")
    if(!isSetupCodeGenerated(name, nfProc)){
        addSetupCodeNames(name, template$makeOtherNames(name, argList), nfProc)
        subList <- template$makeCodeSubList(name, argList)
        addNewCode(name, subList, template, nfProc)
        if(!allowToCpp) addBlockFromCppName(name, nfProc) ## ignores makeOtherNames for now
    }
}

getSymObj_recurse <- function(code, symTab, recurse = FALSE) { #code will be like a$b$c with expectation all are NF or NL
    if(length(code) > 1) {
        if(deparse(code[[1]]) != '$') return(NULL) ## This is valid if we have makeNewNimbleListObject(...)$member or foo()$member
        firstArg <- code[[2]]
        memberName <- deparse(code[[3]])
        symTab <- getSymObj_recurse(firstArg, symTab, recurse = TRUE) ## when recursing, return the symTab
    } else {
        memberName <- deparse(code)
    }
    symObject <- if(is.null(symTab)) NULL else symTab$getSymbolObject(memberName)
    if(recurse) {
        if(inherits(symObject, 'symbolNimbleFunction')) return(symObject$nfProc$setupSymTab)
        if(inherits(symObject, 'symbolNimbleList')) return(symObject$nlProc$symTab) ## can only be known if it was created in setup code
        if(is.null(symObject)) return(NULL) ## will assume to be nimbleList, only kind of nested data structure that could be of unknown type right now
    } else {
        return(symObject)
    }
    stop(paste0('Problem (ii) working through ', deparse(code)), .call = FALSE)
}

symTypeFromSymTab <- function(codeName, symTab, options = character(0) ){
    if(is.language(codeName))
        codeName <- as.character(codeName)
    if(length(codeName) > 1)
        return('NULL')
    class <- class(symTab$symbols[[codeName]])[1]
    if(length(options) == 0)
        return(class)
    if(!(class %in% options))
        return(NULL)  ## nimbleList that was not constructed in setup code 
    return(class)
}

isSymbolType <- function(symTab, varName, symType)
	inherits(symTab$symbols[[varName]], symType)

matchAndFill.call <- function(def, call){
    ##   Note re: matchAndFill.call(function(a = 1, ..., b = 2){}, quote(foo(b = 1, 2, 3))): behavior on this is due to match.call
    theseFormals <- formals(def)
    formalNames <- names(theseFormals) # formalArgs are the arguments that are defined, i.e. does NOT include anything that is from the args "..."
    theseFormals <- theseFormals[nchar(theseFormals) > 0]
    matchedCall <- match.call(def, call) # problem with match.call for our needs is it omits formals that were not provided
    missingArgs <- which(!(names(theseFormals) %in% names(matchedCall)))
    for(ind in missingArgs){ ## this puts back in anything omitted, but order may become wrong
        name <- names(theseFormals)[ind]
        matchedCall[[name]] <- (if(is.null(theseFormals[[name]])) list(NULL) else theseFormals[[name]])
    }

    newCall <- matchedCall[1]

    if(is.null(names(matchedCall))) names(matchedCall) <- c("CALL_", rep("", length(matchedCall) - 1)) ## strangely assigning all "" values results in NULL
    indexAdditionalArgs <- which(!(names(matchedCall)[-1] %in% formalNames))
    
    for(thisArgName in formalNames){					# This is to get the order of the arguments correctly and to insert unmatched arguemnts to ... location if appropriate
        if(thisArgName == '...') {
            for(thisIndex in indexAdditionalArgs) {
                thisName <- names(matchedCall)[thisIndex+1]
                if(thisName=="")
                    newCall[[length(newCall) + 1]] <- matchedCall[[thisIndex + 1]]
                else {
                    newCall[[thisName]] <- matchedCall[[thisName]]
                }
            }
        } else {        
            thisArg <- matchedCall[[thisArgName]]
            if(!is.null(thisArg))
                newCall[[thisArgName]] <- thisArg
        }
    }
    
    return(newCall)
}


determineNdimsFromNfproc <- function(modelExpr, varOrNodeExpr, nfProc) {
    allNDims <- lapply(nfProc$instances, function(x) {
        model <- eval(modelExpr, envir = x)
        if(length(varOrNodeExpr) > 1)
            stop("One must request a node from a model using syntax like `model[[node]]` and not syntax such as `model[[nodes[i]]]`. For the latter case use `values()` instead.")
        if(!exists(as.character(varOrNodeExpr), x, inherits = FALSE) ) {
            stop(paste0('Problem accessing node or variable ', deparse(varOrNodeExpr), '.'), call. = FALSE)
        }
        lab <- eval(varOrNodeExpr, envir = x)
        if(length(lab) != 1)
            stop(paste0("Length of ",
                        deparse(varOrNodeExpr),
                        " requested from ",
                        deparse(modelExpr),
                        " using '[[' is ",
                        length(lab),
                        ". It must be 1." )
               , call. = FALSE)
        varAndIndices <- nimbleInternalFunctions$getVarAndIndices(lab)
        determineNdimFromOneCase(model, varAndIndices)
    } )
    return(allNDims)
}

## from a$b(), goal is to get symbolObject for a
## from a$b$c(), goal is to get symbolObject for b
matchKeywordCodeMemberFun <- function(code, nfProc) {  ## handles cases like a$b(c) as one unit so the member function template for b can be looked up
    dollarSignPart <- code[[1]] ## we already checked that code[[1]][[1]] is '$' before calling this function
    leftSide <- dollarSignPart[[2]]  ## could have further recursing needed. at the moment, nimbleFunction list can only appear at the end of nesting (member data are not part of inheritance)
    rightSide <- dollarSignPart[[3]]
    memFunName <- deparse(rightSide)
    nestedLeftSide <- FALSE

    ## Step 1: Find the symbol object for whatever is on left-hand side of $
    ## This may involve recursion
    ## Only case where the symbol object can be missing (NULL) is nimbleListDef$new() if nimbleListDef if from global env
    if(!is.null(nfProc)) { ## If nfProc is null, we are in an RCfunction and the only valid case is nimbleListDef$new()
        if(length(leftSide) != 1) {
            if(deparse(leftSide[[1]])=='$') {
                symObj <- getSymObj_recurse(leftSide, nfProc$setupSymTab)
                nestedLeftSide <- TRUE
            } else {
                nfNestedCall <- leftSide[[1]]
                if(length(nfNestedCall) != 1) stop(paste0("Cannot handle this expression: ", deparse(code)))
                if(deparse(nfNestedCall) != '[[') stop(paste0("Cannot handle this expression: ", deparse(code)))
                leftLeftSide <- leftSide[[2]]
                if(length(leftLeftSide) != 1) {
                    if(deparse(leftLeftSide[[1]]) == '$') { ##a$b[[i]]$foo()
                        symObj <- getSymObj_recurse(leftLeftSide, nfProc$setupSymTab)
                        nestedLeftSide <- TRUE
                    } else stop('Problem processing something like a$b[[i]]$foo()')
                } else {
                    nfListName <- deparse(leftSide[[2]])
                    if(nfProc$setupSymTab$symbolExists(nfListName)) { ## look in symbol table
                        symObj <- nfProc$setupSymTab$getSymbolObject(nfListName)
                    }
                }
                
            }
        } else {
            symTab <- nfProc$setupSymTab
            symObj <- symTab$getSymbolObject(deparse(leftSide))
        }
    } else {
        symTab <- symObj <- NULL
    }
    
    if(memFunName=='new') { ## this is unique because in non-nested mode, this can be looking for a nlDef in global environment (or possibly elsewhere, but not dealt with)
        ## symObj can be null here
        if(is.null(symObj)) {
            if(nestedLeftSide) stop('Cannot find nested nimbleList definition')
            nlGenName <- deparse(leftSide)
            if(exists(nlGenName, where = globalenv())) {
                possibleNLgen <- get(nlGenName, envir = globalenv())

                if(is.nlGenerator(possibleNLgen)) {
                    thisFunctionMatch <- makeNimbleListTemplateWithBlankFirstArg(nl.getListDef(possibleNLgen))
                } else {
                    stop(paste0('problem with ', deparse(code)))
                }
            } else {
                stop(paste0('problem with ', deparse(code)))
            }
        } else {
            thisFunctionMatch <- makeNimbleListTemplateWithBlankFirstArg(nl.getListDef(symObj$nlProc$nlGenerator))
        }
        for(i in seq_along(code[-1])) code[[i+1]] <- matchKeywords_recurse(code[[i+1]], nfProc)
        code[[ length(code) + 1]] <- leftSide      ## add '.LEFTSIDE = leftSide' arg to code
        names(code)[length(code)] <- '.LEFTSIDE'
        code[[1]] <- as.name('makeNewNimbleListObject') ## modify first arg of code to be desired name of call
        return(matchAndFill.call(thisFunctionMatch, code ) ) ## should create makeNewNimbleListObject( nimbleListCallMaybeNested, var1, var2, etc.)
    }

    if(is.null(symObj)) stop('Problem looking up object')
    
    if(symObj$type == 'nimbleFunction') {
        thisRCfunProc <- symObj$nfProc$RCfunProcs[[memFunName]] 
        if(is.null(thisRCfunProc)) stop(paste0("Cannot handle this expression (member function may not exist): ", deparse(code)), call. = FALSE)
        thisFunctionMatch <- thisRCfunProc$RCfun$template
        return(matchAndFill.call(thisFunctionMatch, code ) )
    } 
    if(inherits(symObj, 'symbolModel')) {
        if(nestedLeftSide) stop('Access to a model cannot be nested.')
        thisFunctionMatch <- matchModelMemberFunctions[[ memFunName ]]
        if(is.null(thisFunctionMatch)) stop(paste0("Cannot handle this expression (looks like a model with an invalid member function call?): ", deparse(code)))
        return(matchAndFill.call(thisFunctionMatch, code) )
    }
    if(inherits(symObj, 'symbolNimbleFunctionList')) {
        thisBaseClass <- symObj$baseClass
        thisFunctionMatch <- environment(symObj$baseClass)$methodList[[memFunName]]$template
        return(matchAndFill.call(thisFunctionMatch, code ) )
    }
    stop(paste0("Cannot handle this expression: ", deparse(code))) 
}


matchKeywordCode <- function(code, nfProc){
    callName <- as.character(code[[1]])
    thisFunctionMatch <- matchFunctions[[ callName ]]
    ## see if this is a member function of an nf object
    if(!is.null(nfProc)) {
        modCallName <- callName
        if(nfProc$setupSymTab$symbolExists(modCallName)) {
            symObj <- nfProc$setupSymTab$getSymbolObject(modCallName)
            if(inherits(symObj, "symbolMemberFunction")) {
                thisRCfunProc <- nfProc$RCfunProcs[[modCallName]]
                if(is.null(thisRCfunProc)) stop(paste0("Cannot handle this expression (looks like a member function but something is wrong): ", deparse(code)), call. = FALSE)
                thisFunctionMatch <- thisRCfunProc$RCfun$template
                return(matchAndFill.call(thisFunctionMatch, code ) )
            }
        }
    }
    
    ## see if this is a call to an RCfunction
    if(is.null(thisFunctionMatch)) {
        if(exists(callName)) {
            callObj <- get(callName)
            if(is.rcf(callObj)) {
                thisFunctionMatch <- callObj
            }
        }
    }
    
    if(!is.null(thisFunctionMatch))
        return(matchAndFill.call(thisFunctionMatch, code ) )
    return(code)
}

matchKeywords_recurse <- function(code, nfProc = NULL) {
    cl = length(code)
    if(cl == 1){ ## There are no arguments
        if(is.call(code)){  
            if(length(code[[1]]) > 1) {
                if(deparse(code[[1]][[1]]) == '$') code <- matchKeywordCodeMemberFun(code, nfProc)
                else
                    code[[1]] <- matchKeywords_recurse(code[[1]], nfProc) ## recurse on the "a$b" part of a$b() (or the "a(b)" part of a(b)()), etc
            } else
                code <- matchKeywordCode(code, nfProc)
        }
        return(code)
    }
    if(length(code[[1]]) == 1) ## a simple call like a(b,c), not a$b(c)
        code <- matchKeywordCode(code, nfProc)
    
    if(is.call(code)) {
        if(length(code[[1]]) > 1) {
            if(deparse(code[[1]][[1]]) == '$') code <- matchKeywordCodeMemberFun(code, nfProc) ## handle a$b(c) as one unit
            else code[[1]] <- matchKeywords_recurse(code[[1]], nfProc) ## handle "a(b)" part of a(b)(c), which is probably *never* triggered
        }
        if(cl >= 2) { ## recurse through arguments
            for(i in 2:cl) {
                code[[i]] <- matchKeywords_recurse(code[[i]], nfProc)
            }
        }
    }
    return(code)
}


makeSingleIndexAccessExpr <- function(newName, newNameExpr) {
    codeNames <- makeSingleIndexAccessCodeNames(newName)
    subList <- lapply(codeNames, as.name)
    ans <- substitute( name[MflatIndex], c(list(name = newNameExpr), subList))
    ans
}

## want map(name, nDim, offset, sizelist, stridelist)
## this is a unique case, where sizelist and stridelist are just lists
## stuck in there
makeMapAccessExpr <- function(newName, newNameExpr, nDim) { ## newNameExpr not used any more!
    codeNames <- makeMapSetupCodeNames(newName)
    subList <- lapply(codeNames, as.name)
    if(nDim == 0) { ## not sure this can happen
        sizeExprList <- strideExprList <- list()
    }
    if(nDim == 1) {
        sizeExprList <- list(substitute(Msizes, subList))
        strideExprList <- list(substitute(Mstrides, subList))
    }
    if(nDim >= 2) {
        sizeExprList <- rep(list( substitute(Msizes[1], subList)), nDim)
        for(i in 1:nDim) sizeExprList[[i]][[3]] <- i
        strideExprList <- rep(list( substitute(Mstrides[1], subList)), nDim)
        for(i in 1:nDim) strideExprList[[i]][[3]] <- i
    }
    ans <- substitute(map( name, nDim, Moffset, sizes, strides),
                      c(subList, list(nDim = nDim, name = newName, sizes = sizeExprList, strides = strideExprList)))
    ans
}

determineNdimFromOneCase <- function(model, varAndIndices) {
    varInfo <- try(model$getVarInfo(as.character(varAndIndices$varName)))
    if(inherits(varInfo, 'try-error')) stop(paste0('In determineNdimFromOneCase: error in extracting varInfo for ', varAndIndices$varName), call. = FALSE)
    varNdim <- varInfo$nDim
    if(length(varAndIndices$indices) == 0) return(varNdim)
    if(length(varAndIndices$indices) != varNdim) {
        stop(paste0('Error, wrong number of dimensions in a node label for ', varAndIndices$varName, '.  Expected ',varNdim,' indices but got ', length(varAndIndices$indices),'.'))
    }
    dropNdim <- sum(unlist(lapply(varAndIndices$indices, is.numeric)))
    return(varNdim - dropNdim)
}



## steps here are similar to makeMapExprFromBrackets, but that uses exprClasses

varAndIndices2mapParts <- function(varAndIndices, sizes, nDim) {
    indices <- varAndIndices$indices
    ## put together offsetExpr, sizeExprs, strideExprs
    ## need sizes to get strides
    if(length(sizes) == 0) sizes <- 1
    if(nDim > 0 & length(indices)==0) {
        blockBool <- rep(TRUE, nDim)
        firstIndexRexprs <- rep(list(1), nDim)
        targetSizes <- sizes
    } else {
        nDim <- length(indices) ## may be redundant/moot
        firstIndexRexprs <- vector('list', nDim)
        targetSizes <- integer(nDim)
        blockBool <- rep(FALSE, nDim)
        for(i in seq_along(indices)) {
            if(is.blank(indices[[i]])) {
                blockBool[i] <- TRUE
                firstIndexRexprs[[i]] <- 1
                targetSizes[i] <- sizes[i]
            }
            else if(is.numeric(indices[[i]])) {
                firstIndexRexprs[[i]] <- indices[[i]]
            } else {
                ## better be :
                if(indices[[i]][[1]] != ":") stop("error, expecting : here")
                blockBool[i] <- TRUE
                firstIndexRexprs[[i]] <- indices[[i]][[2]]
                targetSizes[i] <- indices[[i]][[3]] - indices[[i]][[2]] + 1
            }
        }
    }
    strides <- c(1, cumprod(sizes[-length(sizes)]))
    sourceStrideRexprs <- as.list(strides)
    targetOffsetRexpr <- makeOffsetRexpr(firstIndexRexprs, sourceStrideRexprs)
    targetStrides <- strides[blockBool]
    targetSizes <- targetSizes[blockBool]
    list(offset = eval(targetOffsetRexpr),
         sizes = targetSizes,
         strides = targetStrides)
}


getVarAndIndices <- function(code) {
    if(is.character(code)) code <- parse(text = code, keep.source = FALSE)[[1]]
    if(length(code) > 1) {
        if(code[[1]] == '[') {
            varName <- code[[2]]
            indices <- as.list(code[-c(1,2)])
        } else {
            stop(paste('Error:', deparse(code), 'is a malformed node label.'))
        }
    } else {
        varName <- code
        indices <- list()
    }
    list(varName = varName, indices = indices)
}

## This takes the indices field returned by getVarAndIndices and turns it into a matrix
## e.g. from getVarAndIndices('x[1:3, 2:4]'), we have varName = 'x' and indices = list(quote(1:3), quote(2:4))
## indexExprs2matrix takes the indices and makes [1 3; 2 4]

varAndIndices2flatIndex <- function(varAndIndices, varInfo) {
    if(length(varInfo$maxs) == 0) return(1) ## A -1 is done automatically, later, so here we should stay in R's 1-based indexing
    sizes <- varInfo$maxs
    strides <- c(1, cumprod(sizes[-length(sizes)]))
    flatIndex <- 1 + sum((unlist(varAndIndices$indices)-1) * strides)
    flatIndex
}


makeMapSetupCodeNames <- function(baseName) {
    list(Mstrides = paste0(baseName, '_strides'),
         Msizes = paste0(baseName, '_sizes'),
         Moffset = paste0(baseName, '_offset'))
}


makeSingleIndexAccessCodeNames <- function(baseName) {
    list(MflatIndex = paste0(baseName, '_flatIndex'))
}

handleScaleAndRateForGamma <- function(code){
  ## also handles core R dexp, as well as ddexp, and dinvgamma
  scaleArg <- code$scale
  rateArg <- code$rate
  if(is.null(scaleArg) && is.null(rateArg))	stop('neither scale nor rate defined in dgamma, invgamma, dexp, or ddexp')
  codeName <- deparse(code[[1]])
  dist <- substring(codeName, 2, nchar(codeName))
  if(dist == 'invgamma' || dist == 'sqrtinvgamma') {  ## For [drpq]invgamma
    if(is.null(rateArg)) {
      rateArg <- substitute(1.0/(A), list(A = code$scale)) 
      code$scale <- rateArg
      names(code)[which(names(code) == 'scale')] <- 'rate'  # to preserve correct order
    }
    code$scale <- NULL
  } else if(dist == 'dexp') { ## This is for [drpq]dexp
    ## The logic here is trickier.  scale has a default value and is canonical (what is needed for C).
    ## If rate was provided, then by the time we are here, matchFunctions has been used (matchAndFill.call),
    ## so there will also be a scale, since it has a default
    setScaleArg <- FALSE
    if(is.null(scaleArg)) setScaleArg <- TRUE ## scale is not expected to be null ever (see previous comment), but this is defensive.
    if(!is.null(scaleArg) & !is.null(rateArg)) setScaleArg <- TRUE ## Both are there, so set scale from the provided rate
    if(setScaleArg) {
      code$scale <- NULL
      scaleArg <- substitute(1.0/(A), list(A = code$rate)) 
      code$rate <- scaleArg
      names(code)[which(names(code) == 'rate')] <- 'scale'  # to preserve correct order
    }
    code$rate <- NULL
  } else { # dgamma
    if(is.null(scaleArg)) {
      scaleArg <- substitute(1.0/(A), list(A = code$rate)) 
      code$rate <- scaleArg
      names(code)[which(names(code) == 'rate')] <- 'scale'  # to preserve correct order
    }
    code$rate <- NULL
  }
  return(code)
}

handleScaleAndRateForExpNimble <- function(code){
    scaleArg <- code$scale
    rateArg <- code$rate
    if(is.null(scaleArg) && is.null(rateArg))	stop('neither scale nor rate defined in dexp_nimble')
    if(is.null(rateArg)) {
        rateArg <- substitute(1.0/(A), list(A = code$scale)) 
        code$scale <- rateArg
        names(code)[which(names(code) == 'scale')] <- 'rate'  # to preserve correct order
    }
    code$scale <- NULL
    return(code)
}

## wrtNodes are nodes with respect to which derivatives will be taken
##     (i.e. denominator of dy/dx).  wrtNodes may or may not be part of
##     calcNodes
## extraInputNodes are nodes needed in calculation of derivatives but
##     are not wrtNodes.  They include two categories: (i) parents of
##     any calcNodes that are not themselves calcNodes or wrtNodes;
##     (ii) any stochastic calcNodes, because the node values are needed.
##     However, stochastic calcNodes that are in constantNodes are not
##     included in extraInputNodes.
## constantNodes are nodes that are assumed to be constant for all
##     derivative calls throughout the life of the nimbleFunction.
##     They will be baked in to the CppAD tape and cannot then be changed.
##     They will typically be model data.
## modelOutputNodes are nodes whose values are calculated as part of the
##     tape "forward-zero" stage (value calculation) and need to be stored
##     in the model.  It appears more efficient to copy them in to the model
##     than to use the regular model$calculate() itself.
##     modelOutputNodes will include all deterministic nodes in calcNodes
##     as well as the logProb_ node for any stochastic nodes in calcNodes.
## calcNodes is the same as calcNodes for model$calculate(calcNodes).  It is
##     the ordered sequence of nodes to be calculated
nimDerivsInfoClass_init_impl <- function(.self
                                       , wrtNodes
                                       , calcNodes
                                       , model) {
    .self$model <- model

    ## wrt nodes
    wrtNodes <- model$expandNodeNames(wrtNodes, returnScalarComponents = TRUE)
    wrtNodesAccessor <- modelVariableAccessorVector(model,
                                               wrtNodes,
                                               logProb = FALSE)
    .self$wrtMapInfo <- makeMapInfoFromAccessorVectorFaster(wrtNodesAccessor)

    # See comment in makeModelDerivsInfo for explanation of why both of the next
    # two lines are necessary.
    calcNodes <- model$expandNodeNames(calcNodes)
    calcNodes <- model$expandNodeNames(calcNodes, returnScalarComponents = TRUE)
    derivsInfo <- makeModelDerivsInfo_impl(model,
                                      wrtNodes,
                                      calcNodes,
                                      dataAsConstantNodes = TRUE)
    constantNodes <- derivsInfo$constantNodes
    updateNodes <- derivsInfo$updateNodes
    
    extraInputNodesAccessor <- modelVariableAccessorVector(model,
                                                           updateNodes,
                                                           logProb = FALSE)
    .self$extraInputMapInfo <-
        makeMapInfoFromAccessorVectorFaster(extraInputNodesAccessor)

    constantNodesAccessor <- modelVariableAccessorVector(model,
                                                         constantNodes,
                                                         logProb = FALSE)
    .self$constantMapInfo <- makeMapInfoFromAccessorVectorFaster(constantNodesAccessor)

    
    ## output nodes: deterministic nodes in calcNodes plus logProb nodes
    ##   but not the actual data nodes.

    modelOutputNodes <- makeOutputNodes(model, calcNodes)
    
    modelOutputNodesAccessor <- modelVariableAccessorVector(model,
                                                            modelOutputNodes,
                                                            logProb = FALSE)
    .self$modelOutputMapInfo <-
        makeMapInfoFromAccessorVectorFaster(modelOutputNodesAccessor)
    NULL
}

makeOutputNodes <- function(model,
                            calcNodes) {
  ## Need to do `isDeterm` on nodes, not node components (issue 1431).
  ## Presuming (and based on a simple test) that by now,
  ## input `calcNodes` would include all encompassing nodes
  ## and that this next step would not introduce any additional components.  
  calcNodes <- model$expandNodeNames(calcNodes)
  logProbCalcNodeNames <- model$modelDef$nodeName2LogProbName(calcNodes)
  isDetermCalcNodes <- model$isDeterm(calcNodes, nodesAlreadyExpanded = TRUE)
  modelOutputNodes <- c(model$expandNodeNames(calcNodes[isDetermCalcNodes], returnScalarComponents = TRUE),
                        logProbCalcNodeNames)
  modelOutputNodes
}

nimDerivsInfoClass_output_init_impl <- function(.self,
                                                calcNodes,
                                                model) {
  .self$model <- model
  wrtNodesAccessor <- modelVariableAccessorVector(model,
                                                  character(),
                                                  logProb = FALSE)
  .self$wrtMapInfo <- makeMapInfoFromAccessorVectorFaster(wrtNodesAccessor)
  
  constantNodesAccessor <- modelVariableAccessorVector(model,
                                                       character(),
                                                       logProb = FALSE)
  .self$constantMapInfo <- makeMapInfoFromAccessorVectorFaster(constantNodesAccessor)

  extraInputNodesAccessor <- modelVariableAccessorVector(model,
                                                         character(),
                                                         logProb = FALSE)
  .self$extraInputMapInfo <-
    makeMapInfoFromAccessorVectorFaster(extraInputNodesAccessor)

  modelOutputNodes <- makeOutputNodes(model, calcNodes)
  modelOutputNodesAccessor <- modelVariableAccessorVector(model,
                                                          modelOutputNodes,
                                                          logProb = FALSE)
  .self$modelOutputMapInfo <-
    makeMapInfoFromAccessorVectorFaster(modelOutputNodesAccessor)

  NULL
}

#' Information on model structure used for derivatives
#'
#' Inspect structure of a nimble model to determine nodes needed as "update"
#' and/or "constant" entries in usage of nimDerivs. This will typically be used
#' in the setup code of a nimbleFunction.
#'
#' @param model a nimble model object, such as returned from \code{nimbleModel}.
#' @param wrtNodes a character vector of node names in the model with respect to
#'   which derivatives will be taken through a call to \code{nimDerivs} (same as
#'   \code{derivs}).
#' @param calcNodes a character vector of node names in the model that will be
#'   used in \code{model$calculate(calcNodes)} while derivatives are being
#'   recorded.
#' @param dataAsConstantNodes logical indicating whether data nodes in the model
#'   should automatically be treated as "constant" entries (TRUE) or "update"
#'   entries (FALSE). Defaults to TRUE.
#'
#' @details In the compilable parts of a \code{nimbleFunction} (i.e. \code{run}
#'   or other method code, not \code{setup} code), a call like
#'   \code{nimDerivs(foo(x), ...)} records derivatives of \code{foo(x)}. If
#'   \code{foo} contains any calls to \code{model$calculate(calcNodes)}, it may
#'   be necessary to provide auxiliary information about the model in further
#'   arguments to \code{nimDerivs}, specifically the \code{model},
#'   \code{updateNodes} and \code{constantNodes} arguments.
#'   `makeModelDerivsInfo` is a utility to set up that information for typical
#'   use cases. It returns a list with elements \code{updateNodes} and
#'   \code{constantNodes} to be passed as arguments of the same name to
#'   \code{nimDerivs} (along with passing the \code{model} as the \code{model}
#'   argument).
#'
#' The reason auxiliary information is needed is that recording of derivatives
#' uses a different model than for regular calculations. Together,
#' \code{updateNodes} and \code{constantNodes} should contain all nodes whose
#' values are needed for the model calculations being recorded and that are not
#' part of \code{wrtNodes}. These may include parents of nodes that are in
#' \code{calcNodes} but are not themselves in \code{calcNodes}, as well as the
#' values of stochastic nodes in \code{calcNodes}, which are needed to calculate
#' the corresponding log probabilities. \code{updateNodes} will have their
#' values updated from the regular model every time that recorded derivative
#' calculations are used. \code{constantNodes} will not be updated every time,
#' which means their values will be permanently fixed either the first time the
#' call to `nimDerivs` is invoked or on any subsequent call that has
#' \code{reset=TRUE}. Use of \code{constantNodes} can be slightly more
#' efficient, but one must be careful to be aware that values will not be
#' updated unless \code{reset=TRUE}. See the automatic differentiation section
#' of the User Manual for more information.
#'
#' In the above explanation, care must be taken to understand what should be
#' included in \code{wrtNodes}. In a typical use case, some arguments to
#' \code{foo} are put into the model using \code{values(model, nodes) <-
#' some_foo_arguments}. Next there is typically a call to
#' \code{model$calculate(calcNodes)}. Here the \code{nodes} are considered
#' "with-respect-to" nodes because derivative tracking will follow the arguments
#' of \code{foo}, including when they are put into a model and hence used in
#' \code{model$calculate}. Therefore these \code{nodes} should be the
#' \code{wrtNodes} for \code{makeModelDerivsInfo}.
#'
#' @return A list with elements \code{updateNodes} and \code{constantNodes}.
#'   These shouls be provided as the same-named arguments to \code{nimDerivs}
#'   (same as \code{derivs}).
#'
#' When using double-taping of derivatives (i.e. \code{foo} contains another
#' call to \code{nimDerivs}), both calls to \code{nimDerivs} should include the
#' \code{model}, \code{updateNodes}, and \code{constantNodes} arguments.
#'
#' @export
makeModelDerivsInfo <- function(model,
                           wrtNodes,
                           calcNodes,
                           dataAsConstantNodes = TRUE) {
  wrtNodes <- model$expandNodeNames(wrtNodes, returnScalarComponents = TRUE)
  # This ensures that elements of a non-scalar node become the entire non-scalar node.
  calcNodes <- model$expandNodeNames(calcNodes)
  # And then this splits into scalar components.
  # E.g. if calcNodes is 'mu[1]' but 'mu[1:3]' is a vector node,
  # the above call gets `mu[1:3]` and then the below call splits it.
  calcNodes <- model$expandNodeNames(calcNodes, returnScalarComponents = TRUE) 
  makeModelDerivsInfo_impl(model,
                      wrtNodes,
                      calcNodes,
                      dataAsConstantNodes)
}

getImmediateParentNodes <- function(nodes, model) { 
  ## adapted from BUGS_modelDef creation of edgesFrom2To
  maps <- model$modelDef$maps
  maxNodeID <- length(maps$vertexID_2_nodeID) ## should be same as length(maps$nodeNames)
  
  edgesLevels <- if(maxNodeID > 0) 1:maxNodeID else numeric(0)
  fedgesTo <- factor(maps$edgesTo, levels = edgesLevels) ## setting levels ensures blanks inserted into the splits correctly
  edgesTo2From <- split(maps$edgesFrom, fedgesTo)
  nodeIDs <- model$expandNodeNames(nodes, returnType = "ids")
  fromIDs <- sort(unique(unlist(edgesTo2From[nodeIDs])))
  fromNodes <- maps$graphID_2_nodeName[fromIDs]
  fromNodes
}

makeModelDerivsInfo_impl <- function(model,
                                     wrtNodes,
                                     calcNodes,
                                     dataAsConstantNodes = TRUE) {
  ## Gymnastics to convert to actual nodes for computational efficiency are needed because `setdiff` needs to
  ## operate on node components because `wrt` is in terms of components not nodes.
  nonWrtCalcNodes <- setdiff(calcNodes, wrtNodes)  # Node components here, since `calcNodes`, `wrtNodes` are as components upon input.
  nonWrtCalcActualNodes <- model$expandNodeNames(nonWrtCalcNodes)  # Nodes here. Can be costly for large multivar nodes (say 3 sec. for a 1m-element node).
  nonWrtStochCalcNodes <- nonWrtCalcActualNodes[ model$isStoch(nonWrtCalcActualNodes,
                               nodesAlreadyExpanded = TRUE) ] # Run `isStoch` on nodes not components (issue #1431). 

  parentNodes <- getImmediateParentNodes(calcNodes, model)
  parentNodes <- model$expandNodeNames(parentNodes, returnScalarComponents = TRUE)
  neededParentNodes <- setdiff(parentNodes, c(wrtNodes, nonWrtCalcNodes)) 
    
  extraInputNodes <- model$expandNodeNames(c(neededParentNodes,
                                             nonWrtStochCalcNodes)) # Need as nodes so `isData` below run on nodes for efficiency.
  constantNodes <- character()
  if(dataAsConstantNodes) {
    boolData <- model$isData(extraInputNodes)
    constantNodes <- model$expandNodeNames(extraInputNodes[boolData], returnScalarComponents = TRUE, sort = TRUE) 
    extraInputNodes <- model$expandNodeNames(extraInputNodes[!boolData], returnScalarComponents = TRUE, sort = TRUE)
    ## `wrtNodes` components could have crept in when initializing `extraInputNodes` based on nodes.
    extraInputNodes <- setdiff(extraInputNodes, wrtNodes)
    constantNodes <- setdiff(constantNodes, wrtNodes)  
  } else {
      extraInputNodes <- model$expandNodeNames(extraInputNodes, returnScalarComponents = TRUE, sort = TRUE)
      ## `wrtNodes` components could have crept in when initializing `extraInputNodes` based on nodes.
      extraInputNodes <- setdiff(extraInputNodes, wrtNodes)
  }
  list(updateNodes = extraInputNodes,
       constantNodes = constantNodes)
}

nimDerivsInfoClass_update_init_impl <- function(.self,
                                                updateNodes = NULL,
                                                constantNodes = NULL,
                                                model) {  
  if(is.null(updateNodes)) updateNodes <- character()
  if(is.null(constantNodes)) constantNodes <- character()

  .self$model <- model
  wrtNodesAccessor <- modelVariableAccessorVector(model,
                                                  character(),
                                                  logProb = FALSE)
  .self$wrtMapInfo <- makeMapInfoFromAccessorVectorFaster(wrtNodesAccessor)
  
  constantNodesAccessor <- modelVariableAccessorVector(model,
                                                       constantNodes,
                                                       logProb = FALSE)
  .self$constantMapInfo <- makeMapInfoFromAccessorVectorFaster(constantNodesAccessor)

  extraInputNodesAccessor <- modelVariableAccessorVector(model,
                                                         updateNodes,
                                                         logProb = FALSE)
  .self$extraInputMapInfo <-
    makeMapInfoFromAccessorVectorFaster(extraInputNodesAccessor)
  
  modelOutputNodesAccessor <- modelVariableAccessorVector(model,
                                                          character(),
                                                          logProb = FALSE)
  .self$modelOutputMapInfo <-
    makeMapInfoFromAccessorVectorFaster(modelOutputNodesAccessor)

  NULL
}

nimDerivsInfoClass <- setRefClass(
    'nimDerivsInfoClass',
    fields = list(
        wrtMapInfo = 'ANY'
      , extraInputMapInfo = 'ANY'
      , modelOutputMapInfo = 'ANY'
      , constantMapInfo = 'ANY'
      , model = 'ANY'
    ),
    methods = list(
        initialize = function(wrtNodes = NULL,
                              calcNodes = NULL,
                              updateNodes = NULL,
                              constantNodes = NULL,
                              thisModel = NULL,
                              case = "nimDerivsCalculate",
                              ...) {
          if(!isTRUE(thisModel$modelDef$buildDerivs)) {
            stop(paste0("To use automatic differentiation involving model calculations,\n",
                        "include 'buildDerivs = TRUE' in the call to 'nimbleModel' or set\n",
                        "'nimbleOptions(buildModelDerivs = TRUE)' before calling 'nimbleModel'."), call. = FALSE)
          }
          switch(case,
                 nimDerivsCalculate = nimDerivsInfoClass_init_impl(.self = .self,
                                                                   wrtNodes = wrtNodes,
                                                                   calcNodes = calcNodes,
                                                                   model = thisModel),
                 
                 updateOnly = nimDerivsInfoClass_update_init_impl(.self = .self,
                                                                  updateNodes = updateNodes,
                                                                  constantNodes = constantNodes,
                                                                  model = thisModel),
                 
                 outputOnly = nimDerivsInfoClass_output_init_impl(.self = .self,
                                                                  calcNodes = calcNodes,
                                                                  model = thisModel)
                 )
        }
    )
)

Try the nimble package in your browser

Any scripts or data that you put into this service are public.

nimble documentation built on June 22, 2024, 9:49 a.m.