R/nimbleFunction_Rexecution.R
In nimble: MCMC, Particle Filtering, and Programmable Hierarchical Modeling

Documented in ADbreak any_na any_nan asCol asRow calculate calculateDiff checkInterrupt declare getBound getLogProb getNodeFunctionIndexedInfo getParam makeBoundInfo makeParamInfo nfMethod nfVar nimArray nimC nimCat nimCopy nimInteger nimIntegrate nimLogical nimMatrix nimNumeric nimOptim nimOptimDefaultControl nimPrint nimRep nimSeq nimStop nimSwitch optimDefaultControl pow_int rankSample run.time simulate values

###		These functions are used for calculate/sim/getLP for the nodeFunctionVectors
###		Can either enter model, nodes or model_nodes

#' NIMBLE language function to break tracking of derivatives
#'
#' This function is used in a method of a nimbleFunction that has derivatives enabled.  It returns its value but breaks tracking of derivatives.
#'
#' @param x scalar value
#'
#' @details
#' This funcion only works with scalars.
#'
#' @export
ADbreak <- function(x) x

#' Power function for integer-valued exponent
#'
#' pow function with exponent required to be integer
#'
#' @param a Base
#' @param b Exponent
#'
#' @return a^b
#'
#' @details This is required in nimble models and nimbleFunctions if derivatives will be tracked
#' but tracked only with respect to a, such that b might be any (positive, 0, or negative) integer.
#' This contrasts with pow(a, b) (equivalent to a^b), which requires b > 0 if derivatives will be
#' tracked, even if they will only be requested with respect to a.
#' 
#' @export
pow_int <- function(a, b) a^round(b) # b should be integer.  rounding it makes finite-element derivatives 0, as needed.

#' NIMBLE language functions for R-like vector construction
#'
#' The functions \code{c}, \code{rep}, \code{seq}, \code{which}, \code{diag}, \code{length}, \code{seq_along}, \code{is.na}, \code{is.nan}, \code{any}, and \code{all} can be used in nimbleFunctions and compiled using \code{compileNimble}.
#'
#' @name nimble-R-functions
#'
#' @param ... values to be concatenated.
#' @param x vector of values to be replicated (\code{rep}), or logical array or vector (\code{which}), or object whose length is wanted (\code{length}), or input value (\code{diag}), or vector of values to be tested/checked (\code{is.na}, \code{is.nan}, \code{any}, \code{all}).
#' @param from starting value of sequence.
#' @param to end value of sequence.
#' @param by increment of the sequence.
#' @param length.out desired length of the sequence.
#'
#' @aliases nimC nimRep nimSeq c rep seq which diag length seq_along is.na is.nan any all
#'
#' @details
#' For \code{c}, \code{rep}, \code{seq}, these functions are NIMBLE's version of similar R functions, e.g., \code{nimRep} for \code{rep}.   In a \code{nimbleFunction}, either the R name (e.g., \code{rep}) or the NIMBLE name (e.g., \code{nimRep}) can be used.  If the R name is used, it will be converted to the NIMBLE name. For \code{which}, \code{length}, \code{diag}, \code{seq_along}, \code{is.na}, \code{is.nan}, \code{any}, \code{all} simply use the standard name without \code{"nim"}. These functions largely mimic (see exceptions below) the behavior of their R counterparts, but they can be compiled in a \code{nimbleFunction} using \code{compileNimble}.
#' 
#' \code{nimC} is NIMBLE's version of \code{c} and behaves identically.
#'
#' \code{nimRep} is NIMBLE's version of \code{rep}.  It should behave identically to \code{rep}.  There are no NIMBLE versions of \code{rep.int} or \code{rep_len}.
#'
#' \code{nimSeq} is NIMBLE's version of \code{seq}.  It behaves like \code{seq} with support for \code{from}, \code{to}, \code{by} and \code{length.out} arguments.  The \code{along.with} argument is not supported.  There are no NIMBLE versions of \code{seq.int}, \code{seq_along} or \code{seq_len}, with the exception that \code{seq_along} can take a nimbleFunctionList as an argument to provide the index range of a for-loop (\href{https://r-nimble.org/html_manual/cha-welcome-nimble.html}{User Manual} Ch. 13). 
#'
#' \code{which} behaves like the R version but without support for \code{arr.ind} or \code{useNames} arguments.
#'
#' \code{diag} behaves like the R version but without support for the \code{nrow} and \code{ncol} arguments.
#'
#' \code{length} behaves like the R version.
#' 
#' \code{seq_along} behaves like the R version.
#'
#' \code{is.na} behaves like the R version but does not correctly handle \code{NA} values from R that are type 'logical', so convert these using \code{as.numeric()} before passing from R to NIMBLE.
#' 
#' \code{is.nan} behaves like the R version, but treats \code{NA} of type 'double' as being \code{NaN} and \code{NA} of type 'logical' as not being \code{NaN}. 
#' 
#' \code{any} behaves like the R version but takes only one argument and treats NAs as \code{FALSE}.
#'
#' \code{all} behaves like the R version but takes only one argument and treats NAs as \code{FALSE}.
#'
NULL

#' @rdname nimble-R-functions
#' @export
nimC <- function(...) {
    c(...)
}

#' @rdname nimble-R-functions
#' @export
nimRep <- function(x, ...) {
    rep(x, ...)
}

#' @rdname nimble-R-functions
#' @export
nimSeq <- function(from, to, by, length.out) { ## this creates default arguments filled in at keyword matching that are then useful in cppOutput step to determine if C++ should use nimSeqBy or nimSeqLen
    haveBy <- !missing(by)
    haveLen <- !missing(length.out)
    if(haveBy)
        if(haveLen)
            seq(from, to, by, length.out)
        else
            seq(from, to, by)
    else
        if(haveLen)
            seq(from, to, length.out = length.out)
        else
            seq(from, to)
}

#' Explicitly declare objects created in setup code to be preserved and compiled as member data
#'
#' Normally a nimbleFunction determines what objects from setup code need to be preserved for run code or other member functions.  \code{setupOutputs} allows explicit declaration for cases when an object created in setup code is not used in member functions.
#' 
#' @name setupOutputs
#'
#' @param ... An arbitrary set of names
#'
#' @details
#' Normally any object created in \code{setup} whose name appears in \code{run} or another member function is included in the saved results of setup code.  When the nimbleFunction is compiled, such objects will become member data of the resulting C++ class.  If it is desired to force an object to become member data even if it does not appear in a member function, declare it using \code{setupOutputs}.  E.g., \code{setupOutputs(a, b)} declares that \code{a} and \code{b} should be preserved.
#'
#' The \code{setupOutputs} line will be removed from the setup code.  It is really a marker during nimbleFunction creation of what should be preserved.
#' 
NULL

#' Halt execution of a nimbleFunction function method.  Part of the NIMBLE language
#'
#'
#' @param msg Character object to be output as an error message
#'
#' @aliases stop
#'
#' @author Perry de Valpine
#' @export
#' @details
#' The NIMBLE \code{stop} is similar to the native R \code{stop}, but it takes only one argument, the error message to be output.  During uncompiled NIMBLE execution, \code{nimStop} simply calls R's stop funtion. During compiled execution it calls the error function from the R headers.  \code{stop} is an alias for \code{nimStop} in the NIMBLE language  
nimStop <- function(msg) stop(msg, call. = FALSE)
# we use call.=FALSE because otherwise the error msg indicates the
# error itself occurs in nimStop() and not in the calling frame

#' Time execution of NIMBLE code
#'
#' @param code code to be timed
#'
#' @author NIMBLE Development Team
#' @details
#' Function for use in nimbleFunction run code; when nimbleFunctions are run in R, this simply wraps \code{system.time}.
#' @export
run.time <- function(code) {
    as.numeric(system.time(code)[3])
}

#' Check for interrupt (e.g. Ctrl-C) during nimbleFunction execution. Part of the NIMBLE language.
#'
#' @author Perry de Valpine
#' @export
#' @details
#' During execution of nimbleFunctions that take a long time, it is nice to occassionally check if the user has entered an interrupt and bail out of execution if so.  This function does that.  During uncompiled nimbleFunction execution, it does nothing.  During compiled execution, it calls R_checkUserInterrupt() of the R headers. 
checkInterrupt <- function() {}

#' Turn a numeric vector into a single-row or single-column matrix
#'
#' Turns a numeric vector into a matrix that has 1 row or 1 column.  Part of NIMBLE language.
#'
#' @aliases asCol
#'
#' @param x       Numeric to be turned into a single row or column matrix
#'
#' @author Perry de Valpine
#' @export
#' @details
#' In the NIMBLE language, some automatic determination of how to turn vectors into single-row or single-column matrices is done.
#' For example, in \code{A \%*\% x}, where A is a matrix and x a vector, x will be turned into a single-column matrix unless
#' it is known at compile time that A is a single column, in which case x will be turned into a single-row matrix.
#' However, if it is desired that x be turned into a single row but A cannot be determined at compile time to be a single column,
#' then one can use \code{A \%*\% asRow(x)} to force this conversion.
asRow <- function(x) {
    matrix(x, nrow = 1)
}

## Aliased in asRow
#' @rdname asRow
#' @export
asCol <- function(x) {
    matrix(x, ncol = 1)
}

#' Make an object of information about a model-parameter pairing for getParam.  Used internally
#'
#' Creates a simple getParam_info object, which has a list with a paramID and a type
#'
#' @param model A model such as returned by \code{\link{nimbleModel}}.
#'
#' @param nodes A character string naming one one or more stochastic nodes, such as "mu", "c('mu', 'beta[2]')", or "eta[1:3, 2]".  getParam only works for one node at a time, but if it is indexed (nodes[i]), then makeParamInfo sets up the information for the entire vector nodes.  The processing pathway is used by the NIMBLE compiler.
#'
#' @param param A character string naming a parameter of the distribution followed by node, such as "mean", "rate", "lambda", or whatever parameter names are relevant for the distribution of the node.
#'
#' @param vector A logical indicating whether nodes should definitely be treated as a vector in compiled code, even if it has length = 1.  For type consistency, the compiler needs this option.  If nodes has length > 1, this argument is ignored.
#' 
#' @export
#' @details This is used internally by \code{\link{getParam}}.  It is not intended for direct use by a user or even a nimbleFunction programmer.
makeParamInfo <- function(model, nodes, param, vector = FALSE) {
    ## updating to allow nodes to be a vector. getParam only works for a scalar but in a case like nodes[i] the param info is set up for the entire vector.

    ## this allows for(i in seq_along(nodes)) a <- a + model$getParam(nodes[i], 'mean') through compilation even if some instances have nodes empty and so won't be called.
  if(length(nodes) == 0) return(structure(c(list(paramID = integer()), type = NA, nDim = NA), class = 'getParam_info'))
  
  if(length(param) != 1) stop(paste0(paste0('Problem with param(s) ', paste0(param, collapse = ','), ' while setting up getParam for node ', nodes,
                                            '\nOnly one parameter is allowed.')))
  
  nodeIDs <- model$expandNodeNames(nodes, returnType = 'ids')
  nodeDeclIDs <- model$modelDef$maps$graphID_2_declID[nodeIDs]
  declID2nodeIDs <- split(nodeIDs, nodeDeclIDs)
  numDeclIDs <- length(declID2nodeIDs)
  paramIDs <- integer(numDeclIDs)
  types <- character(numDeclIDs)
  nDims <- integer(numDeclIDs)
  for(i in seq_along(declID2nodeIDs)) {
    nodeIDsFromOneDecl <- declID2nodeIDs[[i]]
    firstNodeName <- model$modelDef$maps$graphID_2_nodeName[nodeIDsFromOneDecl[1]]
    dist <- model$getDistribution(firstNodeName)
    distInfo <- getDistributionInfo(dist)
    paramIDs[i] <- distInfo$paramIDs[param]
    if(is.na(paramIDs[i]))
      stop(paste0("getParam: parameter '", param, "' not found in distribution ",
                  dist, "."))
    types[i] <- distInfo$types[[param]]$type
    nDims[i] <- distInfo$types[[param]]$nDim
  }
  if(length(unique(types)) != 1 || length(unique(nDims)) != 1)
    stop('cannot have an indexed vector of nodes used in getParam if they have different types or dimensions for the same parameter.')
    
  ## on C++ side, we always work with double
  if(types[1] %in% c('integer', 'logical')) types[1] <- 'double'

  if(length(paramIDs) == 1) { ## We could shortcut on this case earlier
    if(vector)
      paramIDvec <- c(-1L, paramIDs[1]) # paramIDs is length 1 anyway
    else
      paramIDvec <- paramIDs[1]
  } else {
    if(length(unique(paramIDs)) == 1) {
      # They are all the same.
      # We encode this in a sparse way
      paramIDvec <- c(-1L, paramIDs[1])
    } else {
      ## Otherwise, create a full vector of paramIDs
      paramIDvec <- rep(1L, length(nodeIDs))
      sourceIDs <- split(seq_along(nodeIDs), nodeDeclIDs)
      for(i in seq_along(declID2nodeIDs)) { #unsplit() would be another approach to this step.
        paramIDvec[sourceIDs[[i]] ] <- paramIDs[i]
      }
    }
  }
    
  ans <- c(list(paramID = paramIDvec), type = types[1], nDim = nDims[1])
  class(ans) <- 'getParam_info'
  ans
}

defaultParamInfo <- function() {
    ans <- list(paramID = integer(), type = 'double', nDim = 0)
    class(ans) <- 'getParam_info'
    ans
}


#' Get value of a parameter of a stochastic node in a model
#'
#' Get the value of a parameter for any single stochastic node in a model.
#'
#' @param model A NIMBLE model object
#'
#' @param node  The name of a stochastic node in the model
#'
#' @param param The name of a parameter for the node
#'
#' @param nodeFunctionIndex For internal NIMBLE use only
#'
#' @param warn For internal NIMBLE use only
#'
#' @export
#' @details
#' Standard usage is as a method of a model, in the form \code{model$getParam(node, param)}, but the usage as a simple function with the model as the first argument as above is also allowed.
#'
#' For example, suppose node 'x[1:5]' follows a multivariate
#' normal distribution (dmnorm) in a model declared by BUGS code.
#' model$getParam('x[1:5]', 'mean') would return the current value of
#' the mean parameter (which may be determined from other nodes).  The
#' parameter requested does not have to be part of the
#' parameterization used to declare the node.  Rather, it can be any
#' parameter known to the distribution.  For example, one can request
#' the scale or rate parameter of a gamma distribution, regardless of
#' which one was used to declare the node.
getParam <- function(model, node, param, nodeFunctionIndex, warn = TRUE) {
    if(missing(param)) { ## already converted by keyword conversion
        stop('Missing either the node or the parameter. Please specify both the node and the parameter of interest. This error may also be triggered if this case of getParam (after keyword replacement) has not been updated for R execution with newNodeFunction system')
        ## nodeFunctionIndex would only be used here, when we make this part work
        nodeFunction <- model
        paramInfo <- node
    } else {
        ## not already converted; this is regular execution
        if(length(node) != 1) stop("getParam only works for one node at a time, but ", length(node), " were provided.")
        tmp <- model$expandNodeNames(node)
        if(length(tmp) > 1) stop("getParam only works for one node at a time, but ", node, " includes multiple nodes.")
        if(!length(tmp)) stop("node ", node, " does not exist.")
        ## makeParamInfo, called by nodeFunctionVector, will check on length of param
        ## nodeFunctionIndex should never be used.
        nfv <- nodeFunctionVector(model, node)
        indexingInfo <- nfv$indexingInfo
        declID <- indexingInfo$declIDs[1] ## should only be one
        nodeFunction <- model$nodeFunctions[[ declID ]] 
        paramInfo <- makeParamInfo(model, node, param)
        if(is.na(paramInfo$type)) stop(paste('getParam called with empty or invalid node:', as.character(node)))
    }
    paramID <- paramInfo$paramID
    if(paramID[1]==-1)
        paramID <- paramID[2]
    nDim <- paramInfo$nDim
    type <- paramInfo$type
    unrolledIndicesMatrixRow <- model$modelDef$declInfo[[declID]]$unrolledIndicesMatrix[ indexingInfo$unrolledIndicesMatrixRows[1], ]
    if(nDim > 2) {
        if(warn)
            messageIfVerbose("  [Warning] `getParam` is not available for parameters of dimension greater than two, but other calculations with models are unaffected.")
        return(NULL)
    }
    funName <- paste0('getParam_',nDim,'D_',type)

    useCompiledNonNestedInterface <- inherits(model, 'CmodelBaseClass') & !getNimbleOption('buildInterfacesForCompiledNestedNimbleFunctions')

    if(useCompiledNonNestedInterface) {
        ans <- model$nodeFunctions[[ declID ]][[1]]$callMemberFunction(model$nodeFunctions[[ declID ]][[2]], funName, paramID, unrolledIndicesMatrixRow)
    } else 
        ans <- eval(substitute(nodeFunction$FUNNAME(paramID, unrolledIndicesMatrixRow), list(FUNNAME = as.name(funName))))
    return(ans)
}

#' Make an object of information about a model-bound pairing for getBound.  Used internally
#'
#' Creates a simple getBound_info object, which has a list with a boundID and a type.
#' Unlike \code{makeParamInfo} this is more bare-bones, but keeping it for parallelism
#' with \code{getParam}.
#'
#' @param model A model such as returned by \code{\link{nimbleModel}}.
#'
#' @param nodes A character string naming a stochastic nodes, such as \code{'mu'}.
#'
#' @param bound A character string naming a bound of the distribution, either \code{'lower'} or \code{'upper'}.
#'
#' @export
#' @details This is used internally by \code{\link{getBound}}.  It is not intended for direct use by a user or even a nimbleFunction programmer.
makeBoundInfo <- function(model, nodes, bound) {
    # note: would be good to work through this and makeParamInfo to ensure that both give good error messages and handle situation where user tries to pass in multiple nodes or multiple params/bounds
    # note: also, it should be fine to combine symbolTable and sizeProcessing stuff for param and bound to avoid repetitive code
    if(length(nodes) > 1) stop("'getBound' only set up to work with a single node'") # I think this is consistent with current getParam behavior
    distInfo <- getDistributionList(model$getDistribution(nodes))
    boundIDvec <- c('lower'=1,'upper'=2)[bound]
    typeVec <- unlist(lapply(distInfo, function(x) x$types[['value']]$type)) 
    ## on C++ side, we always work with double
    if(typeVec[1] %in% c('integer', 'logical')) typeVec[1] <- 'double'
    
    # at moment have bound be scalar regardless of dimension of parameter, as assumed to be same value for all elements of a multivariate distribution; we are not fully handling truncation/bounds for multivariate nodes anyway:
    nDimVec <- 0  # unlist(lapply(distInfo, function(x) x$types[['value']]$nDim))
    ans <- c(list(boundID = boundIDvec, type = typeVec, nDim = nDimVec))
    class(ans) <- 'getBound_info'
    ans
}

#' Get value of bound of a stochastic node in a model
#'
#' Get the value of the lower or upper bound for a single stochastic node in a model.
#'
#' @param model A NIMBLE model object
#'
#' @param node  The name of a stochastic node in the model
#'
#' @param bound Either \code{'lower'} or \code{'upper'} indicating the desired bound for the node
#'
#' @param nodeFunctionIndex For internal NIMBLE use only
#'
#' @export
#' @details
#' Standard usage is as a method of a model, in the form \code{model$getBound(node, bound)}, but the usage as a simple function with the model as the first argument as above is also allowed.
#'
#' For nodes that do not involve truncation of the distribution
#' this will return the lower or upper bound of the distribution, which
#' may be a constant or for a limited number of distributions a parameter
#' or functional of a parameter (at the moment in NIMBLE, the only case
#' where a bound is a parameter is for the uniform distribution. For nodes
#' that are truncated, this will return the desired bound, which may be
#' a functional of other quantities in the model or may be a constant. 
getBound <- function(model, node, bound, nodeFunctionIndex) {
    if(!bound %in% c('lower', 'upper'))
        stop("getBound: 'bound' must be either 'lower' or 'upper'")
    nfv <- nodeFunctionVector(model, node)
    indexingInfo <- nfv$indexingInfo
    declID <- indexingInfo$declIDs[1] ## should only be one
    nodeFunction <- model$nodeFunctions[[ declID ]] 
    boundInfo <- makeBoundInfo(model, node, bound)
    boundID <- boundInfo$boundID
    nDim <- boundInfo$nDim
    type <- boundInfo$type
    unrolledIndicesMatrixRow <- model$modelDef$declInfo[[declID]]$unrolledIndicesMatrix[ indexingInfo$unrolledIndicesMatrixRows[1], ]
    funName <- paste0('getBound_',nDim,'D_',type)

    useCompiledNonNestedInterface <- inherits(model, 'CmodelBaseClass') & !getNimbleOption('buildInterfacesForCompiledNestedNimbleFunctions')

    if(useCompiledNonNestedInterface) {
        ans <- model$nodeFunctions[[ declID ]][[1]]$callMemberFunction(model$nodeFunctions[[ declID ]][[2]], funName, boundID, unrolledIndicesMatrixRow)
    } else 
        ans <- eval(substitute(nodeFunction$FUNNAME(boundID, unrolledIndicesMatrixRow), list(FUNNAME = as.name(funName))))
    return(ans)
}


#' @export
nimSwitch <- function(paramID, IDoptions, ...) {
    dotsList <- eval(substitute(alist(...)))
    iUse <- which(IDoptions == paramID)
    eval(dotsList[[iUse]], envir = parent.frame())
    invisible(NULL)
}

rCalcNodes <- function(model, nfv){ ##nodeFunctionVector
    l_Prob = 0
    model <- nfv$model
    useCompiledNonNestedInterface <- inherits(model, 'CmodelBaseClass') & !getNimbleOption('buildInterfacesForCompiledNestedNimbleFunctions')
    indexingInfo <- nfv$indexingInfo
    declIDs <- indexingInfo$declIDs
    numNodes <- length(declIDs)
    if(numNodes < 1) return(l_Prob)
    unrolledIndicesMatrixRows <- indexingInfo$unrolledIndicesMatrixRows
    for(i in 1:numNodes) {
        declID <- declIDs[i]
        unrolledIndicesMatrixRow <- model$modelDef$declInfo[[declID]]$unrolledIndicesMatrix[ unrolledIndicesMatrixRows[i], ]
        if(useCompiledNonNestedInterface) {
            l_Prob = l_Prob + model$nodeFunctions[[ declID ]][[1]]$callMemberFunction(model$nodeFunctions[[ declID ]][[2]], 'calculate', unrolledIndicesMatrixRow)
        } else
            l_Prob = l_Prob + model$nodeFunctions[[ declID ]]$calculate(unrolledIndicesMatrixRow) ## must use nodeFunctions to have declID ordering
    }
    return(l_Prob)
}

getNodeFunctionIndexedInfo <- function(indexedNodeInfo, iCol) indexedNodeInfo[iCol]

rCalcDiffNodes <- function(model, nfv){
    l_Prob <- 0
    model <- nfv$model
    useCompiledNonNestedInterface <- inherits(model, 'CmodelBaseClass') & !getNimbleOption('buildInterfacesForCompiledNestedNimbleFunctions')
    indexingInfo <- nfv$indexingInfo
    declIDs <- indexingInfo$declIDs
    numNodes <- length(declIDs)
    if(numNodes < 1) return(l_Prob)
    unrolledIndicesMatrixRows <- indexingInfo$unrolledIndicesMatrixRows
    for(i in 1:numNodes) {
        declID <- declIDs[i]
        unrolledIndicesMatrixRow <- model$modelDef$declInfo[[declID]]$unrolledIndicesMatrix[ unrolledIndicesMatrixRows[i], ]
        if(useCompiledNonNestedInterface) {
            l_Prob = l_Prob + model$nodeFunctions[[ declID ]][[1]]$callMemberFunction(model$nodeFunctions[[ declID ]][[2]], 'calculateDiff', unrolledIndicesMatrixRow)
        } else
            l_Prob = l_Prob + model$nodeFunctions[[ declID ]]$calculateDiff(unrolledIndicesMatrixRow) ## must use nodeFunctions to have declID ordering
    }
    return(l_Prob)
}

#' calculate, calculateDiff, simulate, or get the current log probabilities (densities) a set of nodes in a NIMBLE model
#'
#' calculate, calculateDiff, simulate, or get the current log probabilities (densities) of one or more nodes of a NIMBLE model and (for calculate and getLogProb) return the sum of their log probabilities (or densities).  Part of R and NIMBLE.
#' @name nodeFunctions
#' 
#' @param model        A NIMBLE model, either the compiled or uncompiled version
#' @param nodes        A character vector of node names, with index blocks allowed, such as 'x', 'y[2]', or 'z[1:3, 2:4]'
#' @param nodeFxnVector An optional vector of nodeFunctions on which to operate, in lieu of \code{model} and \code{nodes}
#' @param nodeFunctionIndex For internal NIMBLE use only
#' @param includeData  A logical argument specifying whether \code{data} nodes should be simulated into (only relevant for \code{\link{simulate}})
#' @author NIMBLE development team
#' @details
#' Standard usage is as a method of a model, in the form \code{model$calculate(nodes)}, but the usage as a simple function with the model as the first argument as above is also allowed.
#' 
#' These functions expands the nodes and then process them in the model in the order provided.  Expanding nodes means turning 'y[1:2]' into c('y[1]','y[2]') if y is a vector of scalar nodes.
#' Calculation is defined for a stochastic node as executing the log probability (density) calculation and for a deterministic node as calculating whatever function was provided on the right-hand side of the model declaration.
#'
#' Difference calculation (calculateDiff) executes the operation(s) on the model as calculate, but it returns the sum of the difference between the new log probabilities and the previous ones.
#' 
#' Simulation is defined for a stochastic node as drawing a random value from its distribution, and for deterministic node as equivalent to calculate.
#'
#' getLogProb collects and returns the sum of the log probabilities of nodes, using the log probability values currently stored in the model (as generated from the most recent call to calculate on each node)
#'
#' These functions can be used from R or in NIMBLE run-time functions that will be compiled.  When executed in R (including when an uncompiled nimbleFunction is executed), they can be slow because the nodes are expanded each time.  When compiled in NIMBLE, the nodes are expanded only once during compilation, so execution will be much faster.
#'
#' It is common to want the nodes to be provided in topologically sorted order, so that they will be calculated or simulated following the order of the model graph.  Functions such as model$getDependencies(nodes, ...) return nodes in topologically sorted order.  They can be directly sorted by model$topologicallySortNodes(nodes), but if so it is a good idea to expand names first by model$topologicallySortNodes(model$expandNodeNames(nodes))
#'
#' @return calculate and getLogProb return the sum of the log probabilities (densities) of the calculated nodes, with a contribution of 0 from any deterministic nodes
#'
#' @return calculateDiff returns the sum of the difference between the new and old log probabilities (densities) of the calculated nodes, with a contribution of 0 from any deterministic nodes.
#'
#' simulate returns NULL.
#' 
#' @rdname nodeFunctions
#' @export
calculate <- function(model, nodes, nodeFxnVector, nodeFunctionIndex)	
{
    if(!missing(nodeFxnVector)){
        return(rCalcNodes(model, nodeFxnVector))
    }
    if(inherits(model, 'modelBaseClass') ){
        if(missing(nodes) ) 
            nodes <- model$getMaps('nodeNamesLHSall')
        nfv <- nodeFunctionVector(model, nodes)
        return(rCalcNodes(model, nfv))
    }	
}

#' @rdname nodeFunctions
#' @export
calculateDiff <- function(model, nodes, nodeFxnVector, nodeFunctionIndex)		
{
    if(!missing(nodeFxnVector)){
        return(rCalcDiffNodes(model, nodeFxnVector))
    }
    if(inherits(model, 'modelBaseClass') ){
        if(missing(nodes) ) 
            nodes <- model$getMaps('nodeNamesLHSall')
        nfv <- nodeFunctionVector(model, nodes)
        return(rCalcDiffNodes(model, nfv))
    }	
}

rGetLogProbsNodes <- function(model, nfv){
    l_Prob = 0
    model <- nfv$model
    useCompiledNonNestedInterface <- inherits(model, 'CmodelBaseClass') & !getNimbleOption('buildInterfacesForCompiledNestedNimbleFunctions')
    indexingInfo <- nfv$indexingInfo
    declIDs <- indexingInfo$declIDs
    numNodes <- length(declIDs)
    if(numNodes < 1) return(l_Prob)
    unrolledIndicesMatrixRows <- indexingInfo$unrolledIndicesMatrixRows
    for(i in 1:numNodes) {
        declID <- declIDs[i]
        unrolledIndicesMatrixRow <- model$modelDef$declInfo[[declID]]$unrolledIndicesMatrix[ unrolledIndicesMatrixRows[i], ]
        if(useCompiledNonNestedInterface) {
            l_Prob = l_Prob + model$nodeFunctions[[ declID ]][[1]]$callMemberFunction(model$nodeFunctions[[ declID ]][[2]], 'getLogProb', unrolledIndicesMatrixRow)
        } else
            l_Prob = l_Prob + model$nodeFunctions[[ declID ]]$getLogProb(unrolledIndicesMatrixRow) ## must use nodeFunctions to have declID ordering
    }
    return(l_Prob)
}

#' @rdname nodeFunctions
#' @export
getLogProb <- function(model, nodes, nodeFxnVector, nodeFunctionIndex)		
{
    if(!missing(nodeFxnVector)){
        return(rGetLogProbsNodes(model, nodeFxnVector))
    }
    if( inherits(model, "modelBaseClass") ){		
        if(missing(nodes) ) 
            nodes <- model$getMaps('nodeNamesLHSall')
        nfv <- nodeFunctionVector(model, nodes)
    	return(rGetLogProbsNodes(model, nfv))
    }        
}


rSimNodes <- function(model, nfv){
    model <- nfv$model
    useCompiledNonNestedInterface <- inherits(model, 'CmodelBaseClass') & !getNimbleOption('buildInterfacesForCompiledNestedNimbleFunctions')
    indexingInfo <- nfv$indexingInfo
    declIDs <- indexingInfo$declIDs
    numNodes <- length(declIDs)
    if(numNodes < 1) return()
    unrolledIndicesMatrixRows <- indexingInfo$unrolledIndicesMatrixRows
    for(i in 1:numNodes) {
        declID <- declIDs[i]
        unrolledIndicesMatrixRow <- model$modelDef$declInfo[[declID]]$unrolledIndicesMatrix[ unrolledIndicesMatrixRows[i], ]
        if(useCompiledNonNestedInterface) {
            model$nodeFunctions[[ declID ]][[1]]$callMemberFunction(model$nodeFunctions[[ declID ]][[2]], 'simulate', unrolledIndicesMatrixRow)
            
        } else
            model$nodeFunctions[[ declID ]]$simulate(unrolledIndicesMatrixRow) ## must use nodeFunctions to have declID ordering
    }
}

#' @rdname nodeFunctions
#' @export
simulate <- function(model, nodes, includeData = FALSE, nodeFxnVector, nodeFunctionIndex)		
{
    if(!missing(nodeFxnVector)){
        rSimNodes(model, nodeFxnVector)
    }
    else if( inherits(model, "modelBaseClass") ) {
        if(missing(nodes) ) 
            nodes <- model$getMaps('nodeNamesLHSall')
        nfv <- nodeFunctionVector(model, nodes, excludeData = !includeData)
        rSimNodes(nfv$model, nfv)
    }
}


# Fill a vector with flattened values from a set of model nodes
#
# Take a vector of node names for a model and fill a vector by concatenating their values.  Works in R and NIMBLE.
#
# @param vals        the variable in the calling function to receive the values
# @param model       a NIMBLE model object, either compiled or uncompiled
# @param nodes       a vector of node names, allowing index blocks that will be expanded
#
# @author NIMBLE development team
# @export
# @details
# Calling \code{getValues(P, model, nodes)} will modify P in the calling function.  This is being deprecated by
# \code{P <- values(model, nodes)}, but the development syntax of \code{getValues(P, model, nodes)} is still supported.
#
# The result P will be the concatenation of values from the nodes requested.  When requested nodes are from matrices or arrays, the values will be flattened into a vector following column-wise order.
#
# The reverse of \code{getValues(P, model, nodes)} is \code{setValues(P, model, nodes)}, which is being deprecated by \code{values(model, nodes) <- P}.
#
# These functions work in R and in NIMBLE run-time code that can be compiled.
#
# @return NULL, but this function works by the side-effect of modifying P in the calling environment.
getValues <- function(vals, model, nodes, envir = parent.frame()) {
    valsExp = substitute(vals)
    access = modelVariableAccessorVector(model, nodes, logProb = FALSE)
    output = getValuesAccess(access)
    if(is.name(valsExp))
        assign(as.character(valsExp), output, envir = parent.frame(n = 1) ) ## again this assumes valsExp has no indexing -- need to fix
    else {
        varName <- valsExp[[2]]
        orig <- get(varName, envir = envir)
        assignTarget <- valsExp
        assignTarget[[2]] <- quote(orig)
        eval(substitute(AT <- output, list(AT = assignTarget)))
        assign(varName, orig, envir = envir)
    }
}

getValuesAccess <- function(access) {
    fromCode <- makeGetCodeFromAccessorVector(access)
    if(length(fromCode)==0) return(numeric())
    sourceFromObject <- access[[1]]
    unlist(lapply(fromCode, function(i) eval(i)))

}


setValuesAccess <- function(input, access) {
    toCode <- makeSetCodeFromAccessorVector(access)
    mapInfo <- makeMapInfoFromAccessorVector(access) ## inefficient to do both of these, but we need the lengths!
    if(length(toCode)==0) return(NULL)
    sourceToObject <- access[[1]]
    nextIndex <- 0
    ## not easy to check lengths any more
    for(i in 1:length(toCode)) {
        nextLength <- mapInfo[[i]]$length
        oneValue <- input[nextIndex + (1:nextLength)]
        eval(toCode[[i]])
        nextIndex <- nextIndex + nextLength
    }
    invisible(NULL)
}


# Fill a set of nodes in a model from a vector of values
#
# Take a vector of node names for a model and fill them sequentially from the values in a vector.  Works in R and NIMBLE.
#
# @param input        a vector of values to be put in the model's nodes
# @param model       a NIMBLE model object, either compiled or uncompiled
# @param nodes       a vector of node names, allowing index blocks that will be expanded
#
# @author NIMBLE development team
# @export
# @details
# Calling \code{setValues(P, model, nodes)} will place values from P, in order, into the nodes provided for the model.  This is deprecated by
# \code{values(model, nodes) <- P}, but the development syntax of \code{setValues(P, model, nodes)} is still supported.
#
# When provided nodes are from matrices or arrays, the values will filled following column-wise order.
#
# The reverse of \code{setValues(P, model, nodes)} is \code{getValues(P, model, nodes}, which is deprecated by \code{P <- values(model, nodes)}.
#
# These functions work in R and in NIMBLE run-time code that can be compiled.
#
# @return NULL, but this function works by the side-effect of modifying the model.
setValues <- function(input, model, nodes){
	access = modelVariableAccessorVector(model, nodes, logProb = FALSE)
	setValuesAccess(input, access)
}


#' Access or set values for a set of nodes in a model
#'
#' Access or set values for a set of nodes in a NIMBLE model. 
#'
#' @aliases values values<-
#' 
#' @description
#' Get or set values for a set of nodes in a model
#'
#' @param model       a NIMBLE model object, either compiled or uncompiled
#' @param nodes       a vector of node names, allowing index blocks that will be expanded
#' @param value       value to set the node(s) to
#' @param accessorIndex For internal NIMBLE use only
#' 
#' @author NIMBLE development team
#' @export
#' @details
#' Calling \code{values(model, nodes)} returns a vector of the concatenation of values from the nodes requested
#' \code{P <- values(model, nodes)} is a newer syntax for \code{getValues(P, model, values)}, which still works and modifies P in the calling environment.
#'
#' Calling \code{values(model, nodes) <- P} sets the value of the nodes in the model, in sequential order from the vector P.
#'
#' In both uses, when requested nodes are from matrices or arrays, the values will be handled following column-wise order.
#'
#' The older function \code{getValues(P, model, nodes)} is equivalent to \code{P <- values(model, nodes)}, and the older function \code{setValues(P, model, nodes)} is equivalent to \code{values(model, nodes) <- P}
#'
#'
#' These functions work in R and in NIMBLE run-time code that can be compiled.
#'
#' @return A vector of values concatenated from the provided nodes in the model
values <- function(model, nodes, accessorIndex){
	ans <- NA
	getValues(ans, model, nodes, parent.frame())
	ans
}

## # @rdname values

#' @rdname values 
#' @export
`values<-` <- function(model, nodes, accessorIndex, value){
	setValues(value, model, nodes)
	return(model)
}

#' Copying function for NIMBLE
#' 
#' Copies values from a NIMBLE model or modelValues object to another NIMBLE model or modelValues. Work in R and NIMBLE.  The NIMBLE keyword \code{copy} is identical to \code{nimCopy}
#' 
#' @param from		Either a NIMBLE model or modelValues object
#' @param to		Either a NIMBLE model or modelValues object
#' @param nodes		Vector of one or more node names of object \code{from} that will be copied from
#' @param nodesTo	Vector of one or more node names of object \code{to} that will be copied to. If \code{nodesTo} is \code{NULL}, will automatically be set to \code{nodes}
#' @param row		If \code{from} is a modelValues, the row that will be copied from
#' @param rowTo		If \code{to} is a modelValues, the row which will be copied to. If \code{rowTo} is \code{NA}, will automatically be set to \code{row}
#' @param logProb	A logical value indicating whether the log probabilities of the given nodes should also be copied (i.e. if \code{nodes = 'x'}
#' and \code{logProb = TRUE}, then both \code{'x'} and \code{'logProb_x'} will be copied)
#' @param logProbOnly   A logical value indicating whether only the log probabilities of the given nodes should be copied (i.e. if \code{nodes = 'x'}
#' and \code{logProbOnly = TRUE}, then only \code{'logProb_x'} will be copied)
#'
#' @aliases copy
#'
#' @author Clifford Anderson-Bergman
#' @export
#' @details
#'
#' This function copies values from one or more nodes (possibly including log probabilities for nodes) between models and modelValues objects. For modelValues objects, the row must be specified. This function allows one to conveniently copy multiple nodes, avoiding having to write a loop. 
#'
#' @examples
#'	# Building model and modelValues object
#' simpleModelCode <- nimbleCode({
#'	for(i in 1:100)
#'		x[i] ~ dnorm(0,1)
#'})
#' rModel <- nimbleModel(simpleModelCode)
#' rModelValues <- modelValues(rModel)
#'
#' #Setting model nodes
#' rModel$x <- rnorm(100)
#' #Using nimCopy in R.
#' nimCopy(from = rModel, to = rModelValues, nodes = 'x', rowTo = 1)
#' 
#' #Use of nimCopy in a simple nimbleFunction
#' cCopyGen <- nimbleFunction(
#' 	setup = function(model, modelValues, nodeNames){},
#' 	run = function(){
#' 		nimCopy(from = model, to = modelValues, nodes = nodeNames, rowTo = 1)
#' 	}
#' )
#' 
#' rCopy <- cCopyGen(rModel, rModelValues, 'x')
#' \dontrun{
#' cModel <- compileNimble(rModel)
#' cCopy <- compileNimble(rCopy, project = rModel)
#' cModel[['x']] <- rnorm(100)
#' 
#' cCopy$run() ## execute the copy with the compiled function
#' }
nimCopy <- function(from, to, nodes = NULL, nodesTo = NULL, row = NA, rowTo = NA, logProb = FALSE, logProbOnly = FALSE){
    if(is.null(nodes) )
        nodes = from$getVarNames(includeLogProb = logProb) ## allNodeNames(from)
    if( inherits(from, "modelBaseClass") ){
        accessFrom = modelVariableAccessorVector(from, nodes, logProb = logProb, logProbOnly = logProbOnly)
    } else
        if(inherits(from, "modelValuesBaseClass") || inherits(from, "CmodelValues")) {
            accessFrom = modelValuesAccessorVector(from, nodes, logProb = logProb, logProbOnly = logProbOnly)
            if(is.na(row))
                stop("Error: need to supply 'row' for a modelValues copy")
            ##accessFrom$setRow(row) ## NEW ACCESSORS
        }
        else stop('argument "from" in nimCopy is neither a model nor modelValues')

    if( inherits(to, "modelBaseClass") ){
        if(is.null(nodesTo) ) 
            accessTo = modelVariableAccessorVector(to, nodes, logProb = logProb, logProbOnly = logProbOnly)
        else
            accessTo = modelVariableAccessorVector(to, nodesTo, logProb = logProb, logProbOnly = logProbOnly)
    } else
        if(inherits(to, "modelValuesBaseClass") || inherits(to, "CmodelValues")) {
            if(is.null(nodesTo) ) 
                accessTo = modelValuesAccessorVector(to, nodes, logProb = logProb, logProbOnly = logProbOnly)
            else
                accessTo = modelValuesAccessorVector(to, nodesTo, logProb = logProb, logProbOnly = logProbOnly)
            if(is.na(rowTo))
                rowTo = row
            ##accessTo$setRow(rowTo) ## NEW ACCESSORS
        } else stop('argument "to" in nimCopy is neither a model nor modelValues')

    sourceToObject <- accessTo[[1]]
    sourceFromObject <- accessFrom[[1]]
    setCode <- makeSetCodeFromAccessorVector(accessTo)
    getCode <- makeGetCodeFromAccessorVector(accessFrom)
    lengthTo <- length(setCode)
    if(length(getCode) != lengthTo)
    stop('unequal number of entries in nimCopy') 
    if(lengthTo > 0){
        for(i in 1:lengthTo){
            oneValue <- eval(getCode[[i]]) ## may have row hardwired in
            eval(setCode[[i]]) ## oneValue is hard-wired in. may have rowTo hardwired in
        }
    }
}

#' Access or set a member variable of a nimbleFunction
#' 
#' Internal way to access or set a member variable of a nimbleFunction created during \code{setup}.  Normally in NIMBLE code you would use \code{nf$var} instead of \code{nfVar(nf, var)}. 
#'
#' @aliases nfVar nfVar<-
#'
#' @description
#' Access or set a member variable of a specialized nimbleFunction, i.e. a variable passed to or created during the \code{setup} function that is used in run code or preserved by \code{setupOutputs}.  Works in R for any variable and in NIMBLE for numeric variables.
#'
#' @param nf      a specialized nimbleFunction, i.e. a function returned by executing a function returned from \code{nimbleFunction} with \code{setup} arguments
#' @param varName a character string naming a variable in the \code{setup} function.
#' @param value value to set the variable to.
#' @author NIMBLE development team
#' @export
#' @details
#' When \code{nimbleFunction} is called and a \code{setup} function is provided, then \code{nimbleFunction} returns a function.  That function is a generator that should be called with arguments to the \code{setup} function and returns another function with \code{run} and possibly other member functions.  The member functions can use objects created or passed to \code{setup}.  During internal processing, the NIMBLE compiler turns some cases of \code{nf$var} into \code{nfVar(nf, var)}. These provide direct access to setup variables (member data).  \code{nfVar} is not typically called by a NIMBLE user or programmer.
#'
#'
#' For internal access to methods of \code{nf}, see \code{\link{nfMethod}}.
#' 
#' For more information, see \code{?nimbleFunction} and the NIMBLE \href{https://r-nimble.org/html_manual/cha-welcome-nimble.html}{User Manual}.
#'
#' @return whatever varName is in the nimbleFunction nf.
#' @examples
#' nfGen1 <- nimbleFunction(
#'     setup = function(A) {
#'     B <- matrix(rnorm(4), nrow = 2)
#'     setupOutputs(B) ## preserves B even though it is not used in run-code
#'    },
#'    run = function() {
#'       print('This is A', A, '\n')
#' })
#'
#' nfGen2 <- nimbleFunction(
#'   setup = function() {
#'     nf1 <- nfGen1(1000)
#'   },
#'   run = function() {
#'       print('accessing A:', nfVar(nf1, 'A'))
#'       nfVar(nf1, 'B')[2,2] <<- -1000
#'       print('accessing B:', nfVar(nf1, 'B'))
#'    })
#'        
#' nf2 <- nfGen2()
#' nf2$run()
nfVar <- function(nf, varName) {
    refClassObj <- nf_getRefClassObject(nf)
    v <- refClassObj[[varName]]
    return(v)
}

#' @rdname nfVar
#' @export
`nfVar<-` <- function(nf, varName, value) {
    refClassObj <- nf_getRefClassObject(nf)
    refClassObj[[varName]] <- value
    return(nf)
}

#' access (call) a member function of a nimbleFunction
#' 
#' Internal function for accessing a member function (method) of a nimbleFunction.  Normally a user will write \code{nf$method(x)} instead of \code{nfMethod(nf, method)(x)}.
#'
#' @param nf          a specialized nimbleFunction, i.e. one that has already had setup parameters processed
#' @param methodName  a character string giving the name of the member function to call
#'
#' @author NIMBLE development team
#' @export
#' @details
#' nimbleFunctions have a default member function called \code{run}, and may have other member functions provided via the \code{methods} argument to \code{nimbleFunction}.
#' As an internal step, the NIMBLE compiler turns \code{nf$method(x)} into \code{nfMethod(nf, method)(x)}, but a NIMBLE user or programmer would not normally need to use \code{nfMethod} directly.
#'
#' @return a function that can be called.
nfMethod <- function(nf, methodName) {
    # refClassObj <- nf_getRefClassObject(nf)
    # methodCall <- substitute(refClassObj$METHODNAME, list(METHODNAME = methodName))
    # eval(methodCall)
    eval(substitute(nf_getRefClassObject(nf)$METHODNAME, list(METHODNAME = methodName)))
}

#' Generates a weighted sample (with replacement) of ranks
#'
#' Takes a set of non-negative \code{weights} (do not need to sum to 1) and 
#' returns a sample with \code{size} elements of the integers \code{1:length(weights)}, where the probability of being sampled is proportional
#' to the value of \code{weights}. An important note is that the output vector
#' will be sorted in ascending order. Also, right now it works slightly odd syntax (see example below). Later releases of NIMBLE will contain more natural syntax. 
#'
#' If invalid weights provided (i.e. negative weights or weights sum to 1), sets output = rep(1, size) and prints warning. 
#' \code{rankSample} can be used inside nimble functions.
#'
#' @param weights		A vector of numeric weights. Does not need to sum to 1, but must be non-negative
#' @param size			Size of sample
#' @param output		An R object into which the values will be placed. See example below for proper use
#' @param silent Logical indicating whether to suppress logging information
#' @author	Clifford Anderson-Bergman
#' @export
#' @details		
#' \code{rankSample} first samples from the joint distribution \code{size} uniform(0,1) distributions by conditionally sampling from the rank statistics. This leads to 
#' a sorted sample of uniform(0,1)'s. Then, a cdf vector is constructed from weights. Because the sample of uniforms is sorted, \code{rankSample} walks
#' down the cdf in linear time and fills out the sample.
#'  
#' @examples
#' set.seed(1)
#' sampInts = NA	#sampled integers will be placed in sampInts
#' rankSample(weights = c(1, 1, 2), size = 10, sampInts)
#' sampInts
#'# [1] 1 1 2 2 2 2 2 3 3 3
#' rankSample(weights = c(1, 1, 2), size = 10000, sampInts)
#' table(sampInts)
#' #sampInts
#' #   1    2    3 
#' #2434 2492 5074 
#'
#' #Used in a nimbleFunction
#' sampGen <- nimbleFunction(setup = function(){
#' 	x = 1:2
#' },
#' run = function(weights = double(1), k = integer() ){
#' 	rankSample(weights, k, x)
#' 	returnType(integer(1))
#' 	return(x)
#' })
#' rSamp <- sampGen()
#' rSamp$run(1:4, 5)
#' #[1] 3 3 4 4 4
rankSample <- function(weights, size, output, silent = FALSE) {
    ##cat('in R version rankSample\n')
    assign(as.character(substitute(output)), .Call(C_rankSample, as.numeric(weights), as.integer(size), as.integer(output), as.logical(silent)), envir = parent.frame())
}

#' print function for use in nimbleFunctions
#'
#' @param ...  an abitrary set of arguments that will be printed in sequence.
#'
#' @details The keyword \code{print} in nimbleFunction run-time code will be automatically turned into \code{nimPrint}.  This is a function that prints its arguments in order using \code{cat} in R, or using \code{std::cout} in C++ code generated by compiling nimbleFunctions.
#' Non-scalar numeric objects can be included, although their output will be formatted slightly different in uncompiled and compiled nimbleFunctions.
#'
#' For numeric values, the number of digits printed is controlled by the system option \code{nimbleOptions('digits')}.
#'
#' @aliases print
#'
#' @examples
#' ans <- matrix(1:4, nrow = 2) ## R code, not NIMBLE code
#' nimPrint('Answer is ', ans) ## would work in R or NIMBLE
#'
#' @seealso \code{\link{cat}}
#' @export
nimPrint <- function(...) {
    items <- list(...)
    for(i in seq_along(items)) {if(is.array(items[[i]])) print(items[[i]]) else cat(items[[i]])}
    cat('\n')
}

#' cat function for use in nimbleFunctions
#'
#' @param ...  an arbitrary set of arguments that will be printed in sequence.
#'
#' @aliases cat
#'
#' @details
#'
#' \code{cat} in nimbleFunction run-code imitates the R function \code{\link{cat}}.  It prints its arguments in order.  No newline is inserted, so include \code{"\n"} if one is desired.
#'
#' When an uncompiled nimbleFunction is executed, R's \code{cat} is used.  In a compiled nimbleFunction, a C++ output stream is used that will generally format output similarly to R's \code{cat}. Non-scalar numeric objects can be included, although their output will be formatted slightly different in uncompiled and compiled nimbleFunctions.
#'
#' In nimbleFunction run-time code, \code{cat} is identical to \code{print} except the latter appends a newline at the end.
#'
#' \code{nimCat} is the same as \code{cat}, and the latter is converted to the former when a nimbleFunction is defined.
#'
#' For numeric values, the number of digits printed is controlled by the system option \code{nimbleOptions('digits')}.
#'
#' @examples
#' ans <- matrix(1:4, nrow = 2) ## R code, not NIMBLE code
#' nimCat('Answer is ', ans) ## would work in R or NIMBLE
#'
#' @seealso \code{\link{print}}
#' @export
nimCat <- function(...) {
    items <- list(...)
    for(i in seq_along(items)) {if(is.array(items[[i]])) print(items[[i]]) else cat(items[[i]])}
}


#' Creates numeric, integer or logical vectors for use in nimbleFunctions
#'
#' In a \code{nimbleFunction}, \code{numeric}, \code{integer} and \code{logical} are identical to \code{nimNumeric}, \code{nimInteger} and \code{nimLogical}, respectively.
#'
#' @aliases nimInteger nimLogical numeric integer logical
#'
#' @param length the length of the vector (default = 0)
#' @param value value(s) for initializing the vector (default = 0).  This may be a vector, matrix or array but will be used as a vector.
#' @param init logical, whether to initialize elements of the vector (default = TRUE)
#' @param fillZeros logical, whether to initialize any elements not filled by (possibly recycled) \code{value} with 0 (or FALSE for \code{nimLogical}) (default = TRUE)
#' @param recycle logical, whether \code{value} should be recycled to fill the entire \code{length} of the new vector (default = TRUE)
#'
#' @details
#' These functions are similar to R's \code{\link{numeric}}, \code{\link{integer}}, \code{\link{logical}} functions, but they can be used in a nimbleFunction and then compiled using \code{compileNimble}.  Largely for compilation purposes, finer control is provided over initialization behavior.  If \code{init = FALSE}, no initialization will be done, and \code{value}, \code{fillZeros} and \code{recycle} will be ignored.  If \code{init=TRUE} and \code{recycle=TRUE}, then \code{fillZeros} will be ignored, and \code{value} will be repeated (according to R's recycling rule) as much as necessary to fill a vector of length \code{length}.  If \code{init=TRUE} and \code{recycle=FALSE}, then if \code{fillZeros=TRUE}, values of 0 (or FALSE for \code{nimLogical}) will be filled in after \code{value} up to length \code{length}.  Compiled code will be more efficient if unnecessary initialization is not done, but this may or may not be noticeable depending on the situation.
#' 
#' When used in a \code{nimbleFunction} (in \code{run} or other member function), \code{numeric}, \code{integer} and \code{logical} are immediately converted to \code{nimNumeric}, \code{nimInteger} and \code{nimLogical}, respectively.  
#' 
#' @author Daniel Turek, Christopher Paciorek, Perry de Valpine
#' @aliases numeric
#' @seealso \code{\link{nimMatrix}}, \code{\link{nimArray}}
#' @export
nimNumeric <- function(length = 0, value = 0, init = TRUE, fillZeros = TRUE, recycle = TRUE) {
    fillValue <- makeFillValue(value, 'double', init)
    makeReturnVector(fillValue, length, recycle)
}

#' @rdname nimNumeric
#' @export
nimInteger <- function(length = 0, value = 0, init = TRUE, fillZeros = TRUE, recycle = TRUE) {
    fillValue <- makeFillValue(value, 'integer', init)
    makeReturnVector(fillValue, length, recycle)
}

#' @rdname nimNumeric
#' @export
nimLogical <- function(length = 0, value = 0, init = TRUE, fillZeros = TRUE, recycle = TRUE) {
    fillValue <- makeFillValue(value, 'logical', init)
    makeReturnVector(fillValue, length, recycle)
}

#' Creates matrix or array objects for use in nimbleFunctions
#' 
#' In a \code{nimbleFunction}, \code{matrix} and \code{array} are identical to \code{nimMatrix} and \code{nimArray}, respectively
#'
#' @aliases nimArray matrix array
#'
#' @param value value(s) for initialization (default = 0).  This can be a vector, matrix or array, but it will be used as a vector.
#' @param nrow the number of rows in a matrix (default = 1)
#' @param ncol the number of columns in a matrix (default = 1)
#' @param dim vector of dimension sizes in an array (default = \code{c(1, 1)})
#' @param init logical, whether to initialize values (default = \code{TRUE})
#' @param fillZeros logical, whether to initialize any elements not filled by (possibly recycled) \code{value} with 0 (or \code{FALSE} for \code{nimLogical}) (default = \code{TRUE})
#' @param recycle logical, whether \code{value} should be recycled to fill the entire contents of the new object (default = \code{TRUE})
#' @param type character representing the data type, i.e. \code{'double'}, \code{'integer'}, or \code{'logical'} (default = \code{'double'})
#' @param nDim number of dimensions in an array.  This is only necessary for \code{compileNimble} if the length of \code{dim} cannot be determined during compilation.
#'
#' @details
#' These functions are similar to R's \code{\link{matrix}} and \code{\link{array}} functions, but they can be used in a nimbleFunction and compiled using \code{compileNimble}.  Largely for compilation purposes, finer control is provided over initialization behavior, similarly to \code{\link{nimNumeric}}, \code{\link{nimInteger}}, and \code{\link{nimLogical}}. If \code{init = FALSE}, no initialization will be done, and \code{value}, \code{fillZeros} and \code{recycle} will be ignored.  If \code{init=TRUE} and \code{recycle=TRUE}, then \code{fillZeros} will be ignored, and \code{value} will be repeated (according to R's recycling rule) as much as necessary to fill the object.  If \code{init=TRUE} and \code{recycle=FALSE}, then if \code{fillZeros=TRUE}, values of 0 (or FALSE for \code{nimLogical}) will be filled in after \code{value}.  Compiled code will be more efficient if unnecessary initialization is not done, but this may or may not be noticeable depending on the situation.
#'
#' When used in a \code{nimbleFunction} (in \code{run} or other member function), \code{matrix} and \code{array} are immediately converted to \code{nimMatrix} and \code{nimArray}, respectively.
#'
#' The \code{nDim} argument is only necessary for a use like \code{dim <- c(2, 3, 4); A <- nimArray(0, dim = dim, nDim = 3)}.  It is necessary because the NIMBLE compiler must determine during compilation that \code{A} will be a 3-dimensional numeric array.  However, the compiler doesn't know for sure what the length of \code{dim} will be at run time, only that it is a vector.  On the other hand,   \code{A <- nimArray(0, dim = c(2, 3, 4))} is allowed because the compiler can directly determine that a vector of length three is constructed inline for the \code{dim} argument.
#' 
#' @author Daniel Turek and Perry de Valpine
#' @seealso \code{\link{nimNumeric}} \code{\link{nimInteger}} \code{\link{nimLogical}}
#' @export
nimMatrix <- function(value = 0, nrow = NA, ncol = NA, init = TRUE, fillZeros = TRUE, recycle = TRUE, type = 'double') {
    ## the -1's are used because nimble does not allow both missingness and default value
    ## but R's matrix function relies on both possibilities
    fillValue <- makeFillValue(value, type, init)
    mnrow <- missing(nrow) || is.na(nrow)
    mncol <- missing(ncol) || is.na(ncol)
    if(mnrow)
        if(mncol) {
            base::matrix(fillValue)
        } else {
            nrow <- ceiling( length(fillValue) / ncol )
            fillValue <- makeReturnVector(fillValue, nrow * ncol, recycle)
            base::matrix(fillValue, ncol = ncol, nrow = nrow)
        }
    else
        if(mncol) {
            ncol <- ceiling( length(fillValue) / nrow )
            fillValue <- makeReturnVector(fillValue, nrow * ncol, recycle)
            base::matrix(fillValue, nrow = nrow, ncol = ncol)
        } else {
            fillValue <- makeReturnVector(fillValue, ncol*nrow, recycle)
            base::matrix(fillValue, nrow = nrow, ncol = ncol)
        }
}


#' @rdname nimMatrix
#' @export
nimArray <- function(value = 0, dim = c(1, 1), init = TRUE, fillZeros = TRUE, recycle = TRUE, nDim, type = 'double') {
    if(!missing(nDim)) dim <- dim[1:nDim]
    fillValue <- makeFillValue(value, type, init)
    fillValue <- makeReturnVector(fillValue, prod(dim), recycle)
    if(length(dim) == 1) fillValue
    else base::array(fillValue, dim)
}

makeFillValue <- function(value, type, init) {
    fillValue <- if(init) value else 0
    fillValueTyped <- switch(type,
                             double = as.numeric(fillValue),
                             integer = as.integer(fillValue),
                             logical = as.logical(fillValue),
                             stop('unknown type argument'))
    return(fillValueTyped)
}

makeReturnVector <- function(fillValue, length, recycle) {
    if(length(fillValue) == 1) {
        if(recycle)
            rep(fillValue, length)
        else
            c(fillValue, as(rep(0, max(length-1, 0)), class(fillValue)))
    }
    else {
        if(length(fillValue) != length) {
            if(length(fillValue) < length) {
                ##warning(paste0("Not enough values provided for vector of length ",length, ".")) 
                if(recycle)
                    rep(fillValue, length.out = length)
                else
                    c(fillValue, as(rep(0, length-length(fillValue)), class(fillValue)))
            } else {
                ##warning(paste0("Too many values provided for vector of length ",length, ".")) 
                fillValue[1:length]
            }
        } else {
            fillValue
        }
    }
}


#' Explicitly declare a variable in run-time code of a nimbleFunction
#'
#' Explicitly declare a variable in run-time code of a nimbleFunction, for cases when its dimensions cannot be inferred before it is used.  Works in R and NIMBLE.
#'
#' @param name     Name of a variable to declare, without quotes
#' @param def      NIMBLE type declaration, of the form \code{TYPE(nDim, sizes)}, where \code{TYPE} is \code{integer}, \code{double}, or \code{logical}, \code{nDim} is the number of dimensions, and \code{sizes} is an optional vector of sizes concatenated with \code{c}.  If \code{nDim} is omitted, it defaults to 0, indicating a scalar.  If sizes are provided, they should not be changed subsequently in the function, including by assignment.  Omitting \code{nDim} results in a scalar.  For \code{logical}, only scalar is currently supported.
#'
#' @author NIMBLE development team
#' @export
#' @details
#' In a run-time function of a nimbleFunction (either the \code{run} function or a function provided in \code{methods} when calling \code{nimbleFunction}), the dimensionality and numeric type of a variable is inferred when possible from the statement first assigning into it.  E.g. \code{A <- B + C} infers that \code{A} has numeric types, dimensions and sizes taken from \code{B + C}.  However, if the first appearance of \code{A} is e.g. \code{A[i] <- 5}, \code{A} must have been explicitly declared.  In this case, \code{declare(A, double(1))} would make \code{A} a 1-dimensional (i.e. vector) double.
#'
#' When sizes are not set, they can be set by a call to \code{setSize} or by assignment to the whole object.  Sizes are not automatically extended if assignment is made to elements beyond the current sizes.  In compiled nimbleFunctions doing so can cause a segfault and crash the R session.
#'
#'
#' This part of the NIMBLE language is needed for compilation, but it also runs in R.  When run in R, is works by the side effect of creating or modifying \code{name} in the calling environment.
#'
#' @examples
#' declare(A, logical())             ## scalar logical, the only kind allowed
#' declare(B, integer(2, c(10, 10))) ## 10 x 10 integer matrix
#' declare(C, double(3))             ## 3-dimensional double array with no sizes set.
declare <- function(name, def){
    defCode <- substitute(def)
    name <- substitute(name)
    if(exists(as.character(name), parent.frame(), inherits = FALSE)) return(invisible(NULL))
    value <- if(defCode[[1]] == 'logical') FALSE else 0
    if(length(defCode) == 1){  ## no arg, like double()
        assign(as.character(name), value, envir = parent.frame() )
        return()
    }
    nDim = eval(defCode[[2]], envir = parent.frame() )
    if(nDim == 0 ){ ## like double(0)
        assign(as.character(name), value, envir = parent.frame() )
        return()
    }
    dims = rep(1, nDim)
    if(length(defCode) == 3) ## notation like double(2, c(3, 5))
        dims = eval(defCode[[3]], envir = parent.frame() )
    else {
        if(length(defCode) == 2 + nDim) {
            dims <- numeric(nDim)
            for(i in 1:nDim)
                dims[i] <- eval(defCode[[2 + i]], envir = parent.frame())
        }
    }
    if(length(dims) != nDim)
        stop('in declare, dimensions are not declared properly')
    if(prod(dims)==1) {
        if(length(dims) == 1)
            return(assign(as.character(name), value, envir = parent.frame()))
    }
    assign(as.character(name), array(value, dim = dims), envir = parent.frame() )
}

#' Determine if any values in a vector are NA or NaN
#'
#' NIMBLE language functions that can be used in either compiled or uncompiled
#' nimbleFunctions to detect if there are any NA or NaN values in a vector.
#'
#' @param x vector of values
#'
#' @aliases any_nan
#' @author NIMBLE Development Team
#' @export
any_na <- function(x) any(is.na(x))

## Deprecated as of v. 1.1.0
is.na.vec <- function(x)
    stop("`is.na.vec` has been removed from NIMBLE; please use `any_na`")

#' @rdname any_na
#' @export
any_nan <- function(x) any(is.nan(x))

## Deprecated as of v. 1.1.0
is.nan.vec <- function(x)
    stop("`is.nan.vec` has been removed from NIMBLE; please use `any_nan`")

#' @export
nimRound <- round

#' General-purpose Optimization
#'
#' NIMBLE wrapper around R's builtin \code{\link{optim}}, with flexibility for
#' additional methods.
#'
#' @param par Initial values for the parameters to be optimized over.
#' @param fn A function to be minimized (or maximized), with first argument the
#'   vector of parameters over which minimization is to take place. It should
#'   return a scalar result.
#' @param gr A function to return the gradient for the "BFGS", "CG" and
#'   "L-BFGS-B" methods. If not provided, a finite-difference approximation to
#'   derivatives will be used.
#' @param he A function to return the Hessian matrix of second derivatives. Used
#'   (but not required) in "nlminb" or (optionally) user-provided methods.
#' @param ... IGNORED
#' @param method The method to be used. See `Details` section of
#'   \code{\link{optim}}. One of: "Nelder-Mead", "BFGS", "CG", "L-BFGS-B", or
#'   "nlminb". Note that the R methods "SANN", "Brent" are not supported. It is
#'   also possible to provide a new method; see details.
#' @param lower Vector or scalar of lower bounds for parameters.
#' @param upper Vector or scalar of upper bounds for parameters.
#' @param control A list of control parameters. See \code{Details} section of
#'   \code{\link{optim}}. For code in a nimbleFunction to be compiled, this must
#'   be an \code{optimControlNimbleList}, which has fields for most elements in
#'   the control list for R's \code{optim}.
#' @param hessian Logical. Should a Hessian matrix be returned?
#'
#' @details This function for use in nimbleFunctions for compilation by
#'   \code{compileNimble} provides capabilities similar to R's \code{optim} and
#'   \code{nlminb}. For the supported methods provided by \code{optim}, a
#'   compiled nimbleFunction will directly call the C code used by R for these
#'   methods.
#'
#' If \code{optim} appears in a nimbleFunction, it will be converted to
#' \code{nimOptim}.
#'
#' Note that if a gradient function (\code{gr}) is not provided, \code{optim}
#'   provides a finite difference approximation for use by optimization methods
#'   that need gradients. nimble's compiled version of \code{nimOptim} does the
#'   same thing, although results might not be completely identical.
#'
#' For \code{method="nlminb"}, a compiled nimbleFunction will run R's
#' \code{nlminb} directly in R, with \code{fn}, \code{gr} (if provided) and
#' \code{he} (if provided) that call back into compiled code.
#'
#' An experimental feature is the capability to provide one's own optimization
#' method in R and register it for use by \code{nimOptim}. One must write a
#' function that takes arguments \code{par}, \code{fn}, \code{gr}, \code{he},
#' \code{lower}, \code{upper}, \code{control}, and \code{hessian}. The function
#' must return a list with elements \code{par}, \code{value},
#' \code{convergence}, \code{counts}, \code{evaluations}, \code{message}, and
#' \code{hessian} (which may be NULL). If \code{hessian=TRUE} but the function
#' does not return a matrix in the \code{hessian} element of its return list,
#' \code{nimOptim} will fill in that element using finite differences of the
#' gradient.
#'
#' The \code{control} list passed from a nimbleFunction to the
#' optimization function will include a minimum of options, including
#' \code{abstol}, \code{reltol}, \code{maxit}, and \code{trace}. Other options
#' for a specific method may be set within the custom optimization function but
#' cannot be passed from \code{nimOptim}.
#'
#'  The elements \code{parscale} and \code{fnscale} in \code{control} are used in
#' a special way. They are implemented by \code{nimOptim} such that for *any*
#' the method is expected to do minimization and \code{nimOptim} will arrange
#' for it to minimize \code{fn(par)/fnscale} in the parameter space
#' \code{par/parscale}.
#'
#' An optimizer \code{fun} may be registered by
#' \code{nimOptimMethod("method_name", fun)}, and then "\code{method_name}" can
#' be used as the \code{method} argument to \code{nimOptim} to use \code{fun}.
#' An optimizer may be found by \code{nimOptimMethod("method_name")} and may be
#' removed by \code{nimOptimMethod("method_name", NULL)}.
#'
#' Support for \code{method="nlminb"} is provided in this way, and can be
#' studied as an example via \code{nimOptimMethod("nlminb")}.
#'
#' The system for providing one's own optimizer is not considered stable and is
#' subject to change in future versions.
#'
#' @return \code{\link{optimResultNimbleList}}
#' @seealso \code{\link{optim}}
#' @export
#' @examples
#' \dontrun{
#' objectiveFunction <- nimbleFunction(
#'     run = function(par = double(1)) {
#'         return(sum(par) * exp(-sum(par ^ 2) / 2))
#'         returnType(double(0))
#'     }
#' )
#' optimizer <- nimbleFunction(
#'     run = function(method = character(0), fnscale = double(0)) {
#'         control <- optimDefaultControl()
#'         control$fnscale <- fnscale
#'         par <- c(0.1, -0.1)
#'         return(optim(par, objectiveFunction, method = method, control = control))
#'         returnType(optimResultNimbleList())
#'     }
#' )
#' cOptimizer <- compileNimble(optimizer)
#' cOptimizer(method = 'BFGS', fnscale = -1)
#' }
nimOptim <- function(par, fn, gr = "NULL", he = "NULL", ..., method = "Nelder-Mead", lower = -Inf, upper = Inf,
                     control = nimOptimDefaultControl(), hessian = FALSE) {
    ## Tweak parameters.
    if(identical(gr, "NULL")) gr <- NULL
    if(identical(he, "NULL")) he <- NULL
  ## We don't handle creation of finite difference versions of gr and he here
  ## as we do in compiled code. That is a compiled/uncompiled behavior discrepancy.
    defaultControl <- nimOptimDefaultControl()
    Rcontrol <- list()
    ## Only enter non-default values into Rcontrol.
    ## R's optim will fill in the rest.
    ## This should match compiled handling of control list
    ## parameters in NimOptimProblem::solve in nimOptim.cpp.
    for (name in defaultControl$nimbleListDef$types$vars) {
        if (!identical(control[[name]], defaultControl[[name]])) {
            Rcontrol[[name]] <- control[[name]]
        }
    }
    if(method %in% c("Nelder-Mead", "CG", "BFGS", "L-BFGS-B")) {
      result <- optim(par, fn, gr = gr, ..., method = method,
                      lower = lower, upper = upper, control = Rcontrol, hessian = hessian)
    } else {
      # NB: Compiled nimOptim applies parscale and fnscale internally,
      # so they are applied for all optimizers in the same way.
      # Hence we have to do that here.
      fnscale <- control$fnscale
      if(is.null(fnscale) || is.na(fnscale))
        fnscale <- 1
      working_parscale <- control$parscale
      if(is.null(working_parscale) ||
           ((length(working_parscale)==1) && is.na(working_parscale[1])))
        working_parscale <- rep(1, length(par))
      adj_fn <- function(par, ...) fn(par*working_parscale, ...) / fnscale
      if(is.function(gr))
        adj_gr <- function(par, ...) working_parscale * gr(par*working_parscale, ...) / fnscale
      else adj_gr <- NULL
      if(hessian || is.function(he))
        working_parscale_mat <- outer(working_parscale, working_parscale)
      if(is.function(he)) {
        adj_he <- function(par, ...) working_parscale_mat * he(par*working_parscale, ...) / fnscale
      }
      else adj_he <- NULL
      init_par <- par / working_parscale
      result <- custom_optim_inner(method, init_par, adj_fn, adj_gr, adj_he,
                                   lower, upper, control, hessian)
      hessian_returned <- is.numeric(result$hessian) && length(result$hessian)
      if(hessian && !hessian_returned)
        result$hessian <- optimHess(result$par, adj_fn, adj_gr) # ignore ndeps
      hessian_returned <- is.numeric(result$hessian) && length(result$hessian)

      result$value <- result$value * fnscale
      result$par <- result$par * working_parscale
      if(hessian_returned)
        result$hessian <- fnscale * result$hessian / working_parscale_mat
      ## Rcontrol$type <- NULL
      ## optimr_methods <- optimx::ctrldefault(4) # imitated optimx::checkallsolvers
      ## avail_in_optimr <- suppressWarnings(
      ##   optimx::checksolver(method, optimr_methods$allmeth, optimr_methods$allpkg)
      ## )
      ## if(!is.null(avail_in_optimr)) {
      ##   resultRaw <- try(optimx::optimr(par, fn, gr = gr, he = he, method = method, lower = lower,
      ##                                   upper = upper, hessian=hessian))
      ## } else {
      ##   optimizer <- nimOptimMethod(method)
      ##   if(is.null(optimizer)) stop("optimizer ", method, " not found. See help(nimOptimMethod).")
      ##   resultRaw <- try(optimizer(par, fn, gr = gr, he = he, lower = lower,
      ##                                   upper = upper, control = Rcontrol, hessian=hessian))
      ## }
      ## if(inherits(resultRaw, "try-error")) stop("problem with optimizer ", method)

      ## result <- list()
      ## for(name in c("par", "value", "convergence", "counts", "message")) {
      ##   result[[name]] <- resultRaw[[name]]
      ##   attributes(result[[name]]) <- NULL
      ## }
    }
    nimResult <- do.call(optimResultNimbleList$new, result)
    # Tweak result value to exactly match C++ behavior.
    nimResult$counts <- unname(nimResult$counts)
    if (is.null(nimResult$message)) {
        nimResult$message <- ''
    }
    if (is.null(nimResult$hessian)) {
        nimResult$hessian <- matrix(nrow = 0, ncol = 0)
    }
    return(nimResult)
}

# nimble_custom_optim_ function is called only from C++
custom_optim <- function(method, par, lower, upper, control,
                         hessian, use_gr, use_he, extptr,...) {
    
  fnsym <- nimbleUserNamespace$sessionSpecificDll$CALL_NimOptimProblem_fn # getNativeSymbolInfo("CALL_NimOptimProblem_fn")
  fn <- function(p) eval(call('.Call', fnsym, p, extptr))
  if(use_gr) {
    grsym <- nimbleUserNamespace$sessionSpecificDll$CALL_NimOptimProblem_gr # getNativeSymbolInfo("CALL_NimOptimProblem_gr")
    gr <- function(p) eval(call('.Call', grsym, p, extptr)) 
  } else gr <- NULL
  if(use_he) {
    hesym <- nimbleUserNamespace$sessionSpecificDll$CALL_NimOptimProblem_he # getNativeSymbolInfo("CALL_NimOptimProblem_he")
    he <- function(p) eval(call('.Call', hesym, p, extptr))
  } else he <- NULL
  custom_optim_inner(method, par, fn, gr, he, lower, upper, control, hessian)
}

custom_optim_inner <- function(method, par, fn, gr, he, lower, upper, control,
                               hessian) {
                                        # To do:
  # uncompiled R execution
  # the control list here will have been packed up like this from the optimControlNimbleList:
  # list(parscale = control$parscale, fnscale = control$fnscale, maxit = control$maxit,
  #      trace = control$trace, abstol = control$abstol, reltol = control$reltol)
  ## We experimented with giving optimx a priority of call, but we are not doing that for now.
  ## optimr_methods <- optimx::ctrldefault(4) # imitated optimx::checkallsolvers
  ## avail_in_optimr <- suppressWarnings(
  ##   optimx::checksolver(method, optimr_methods$allmeth, optimr_methods$allpkg)
  ## )
  ## if(!is.null(avail_in_optimr)) {
  ##   resultRaw <- try(optimx::optimr(par, fn, gr = gr, he = he, method = method, lower = lower,
  ##                                   upper = upper, hessian=hessian))
  ## } else {
  optimizer <- nimOptimMethod(method)
  if(is.null(optimizer)) stop("optimizer ", method, " not found. See help(nimOptimMethod).")
  resultRaw <- try(optimizer(par, fn, gr = gr, he = he, lower = lower,
                             upper = upper, control = control, hessian=hessian))
  ##  }
  if(inherits(resultRaw, "try-error")) stop("problem with optimizer ", method)

  result <- list()
  for(name in c("par", "value", "convergence", "counts", "message", "hessian")) {
    result[[name]] <- resultRaw[[name]]
    attributes(result[[name]]) <- NULL
  }
  result
}

#' Creates a deafult \code{control} argument for \code{\link{optim}} (just an empty list).
#'
#' @return an empty \code{list}.
#' @seealso \code{\link{nimOptim}}, \code{\link{optim}}
#' @export
optimDefaultControl <- function() {
    return(list())
}

#' Creates a default \code{control} argument for \code{\link{nimOptim}}.
#'
#' @return \code{\link{optimControlNimbleList}}
#' @seealso \code{\link{nimOptim}}, \code{\link{optim}}
#' @export
nimOptimDefaultControl <- function() {
  control <- optimControlNimbleList$new()
  control$trace <- 0
  control$fnscale <- 1
  control$parscale <- NA # Must be filled in to length of par
  control$ndeps <- NA    # Ditto
  control$maxit <- NA  ## The default value depends on method.
  control$abstol = -Inf
  control$reltol <- sqrt(.Machine$double.eps)
  control$alpha <- 1.0
  control$beta <- 0.5
  control$gamma <- 2.0
  control$REPORT <- NA # Method dependent and not used in compiled version
  control$type <- 1
  control$lmm <- 5
  control$factr <- 1e7
  control$pgtol <- 0
  control$tmax <- 10
  control$temp <- 10
  return(control)
}

#' Integration of One-Dimensional Functions
#'
#' NIMBLE wrapper around R's builtin \code{\link{integrate}}. Adaptive quadrature
#' of functions of one variable over a finite or infinite interval.
#'
#' @param f nimbleFunction of one input for which the integral is desired.
#' See below for details on requirements for how \code{f} must be defined.
#' @param lower an optional scalar lower bound for the input of the function.
#' @param upper an optional scalar upper bound for the input of the function.
#' @param param additional parameter(s) to the function
#'               that are fixed with respect to the integration. If \code{f}
#'               takes no additional arguments (beyond the variable of
#'               integration), this must be provided but need not be used in
#'               \code{f}. Can be of length one or more.
#'               parameters.
#' @param subdivisions the maximum number of subintervals.
#' @param rel.tol relative accuracy requested.
#' @param abs.tol absolute accuracy requested.
#' @param stop.on.error logical. If \code{TRUE} (the default) an error stops the
#'        function. Otherwise some errors will give a result with the error
#'        code given in the third element of the result vector.
#'
#' @return A vector with three values, the first the estimate of the integral,
#' the second an estimate of the modulus of the absolute error, and the third
#' a result code corresponding to the \code{message} returned by \code{integrate}.
#' The numerical result code can be interpreted as follows:
#' \itemize{
#' \item \code{0}: "OK"
#' \item \code{1}: "maximum number of subdivisions reached"
#' \item \code{2}: "roundoff error was detected"
#' \item \code{3}: "extremely bad integrand behaviour"
#' \item \code{4}: "roundoff error is detected in the extrapolation table"
#' \item \code{5}: "the integral is probably divergent"
#' \item \code{6}: "the input is invalid"
#' }
#' 
#' @details
#' The function \code{f} should take two arguments, the first of type
#' \code{double(1)}, i.e., vector. \code{f} should be vectorized
#' in that it should also return a \code{double(1)} object, containing
#' the result of applying the function to each element of the first
#' argument. (The result can be calculated using vectorized NIMBLE code or
#' using a loop.) The second argument is required to also be of type
#' \code{double(1)}, containing any additional parameter(s) to the function
#' that are not being integrated over. This argument can be unused in
#' the function if the function does not need additional parameters.
#' Note that this must be of type \code{double(1)} even if \code{param}
#' contains a single element (NIMBLE will manage the lengths behind the
#' scenes).
#'
#' Note that unlike with R's \code{integrate}, additional parameters
#' must be passed as part of a vector, specified via \code{param},
#' and cannot be passed as individual named arguments.
#' 
#' @seealso \code{\link{integrate}}
#' @author Christopher Paciorek, Paul van Dam-Bates, Perry de Valpine
#' @aliases integrate
#' @export
#' @examples
#' integrand <- nimbleFunction(
#'    run = function(x = double(1), theta = double(1)) {
#'        return(x*theta[1])   
#'    returnType(double(1))
#'  }
#')
#'
#' fun <- nimbleFunction(
#'    run = function(theta = double(0), lower = double(0), upper = double(0)) {
#'        param = c(theta, 0)  # cannot be scalar, so pad with zero.
#'        output = integrate(integrand, lower, upper, param)
#'        returnType(double(1))
#'        return(output)
#'  })
#'
#' fun(3.1415927, 0, 1)
#' \dontrun{
#' cfun <- compileNimble(fun)
#' cfun(3.1415927, 0, 1)
#' }
#' 
nimIntegrate <- function(f, lower, upper, param, subdivisions = 100L,
                         rel.tol = .Machine$double.eps^0.25, abs.tol = .Machine$double.eps^0.25, stop.on.error = TRUE) {
  # R's integrate has default abs.tol = rel.tol, but nimble has now way to do that kind of chain of defaults
  # in compiled code.
  output <- rep(0, 3)
    messages <- c("OK", "maximum number of subdivisions reached", 
        "roundoff error was detected", "extremely bad integrand behaviour", 
        "roundoff error is detected in the extrapolation table", 
        "the integral is probably divergent", "the input is invalid")
    result <- integrate( f, lower, upper, param, subdivisions = subdivisions,
              rel.tol = rel.tol, abs.tol = abs.tol,
              stop.on.error = stop.on.error)
    output[1] <- result$value
    output[2] <- result$abs.error
    output[3] <- which(result$message == messages) - 1
    return(output)           
}
Any scripts or data that you put into this service are public.
nimble documentation built on Sept. 11, 2024, 7:10 p.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
nimble
MCMC, Particle Filtering, and Programmable Hierarchical Modeling

R/nimbleFunction_Rexecution.R
In nimble: MCMC, Particle Filtering, and Programmable Hierarchical Modeling

Try the nimble package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

nimble MCMC, Particle Filtering, and Programmable Hierarchical Modeling

R/nimbleFunction_Rexecution.R In nimble: MCMC, Particle Filtering, and Programmable Hierarchical Modeling

Try the nimble package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

nimble
MCMC, Particle Filtering, and Programmable Hierarchical Modeling

R/nimbleFunction_Rexecution.R
In nimble: MCMC, Particle Filtering, and Programmable Hierarchical Modeling
