diff --git a/packages/.gitignore b/packages/.gitignore index 2b0023b39..a7c4e0931 100644 --- a/packages/.gitignore +++ b/packages/.gitignore @@ -14,3 +14,5 @@ Eigen_local* config.* *.Rproj profile/ + +.vscode diff --git a/packages/nimble/NAMESPACE b/packages/nimble/NAMESPACE index 943031b9d..f8b0b2207 100644 --- a/packages/nimble/NAMESPACE +++ b/packages/nimble/NAMESPACE @@ -18,6 +18,7 @@ S3method(as.list, modelValuesBaseClass) S3method(length, nimPointerList) export(calc_dmnormConjugacyContributions) export(calc_dmnormAltParams) +export(calc_dmnorm_inv_ld_AltParams) export(calc_dwishAltParams) export(calc_dcatConjugacyContributions) export(CAR_calcM) @@ -103,6 +104,7 @@ export(dinvgamma) export(dinvwish_chol) export(dlkj_corr_cholesky) export(dmnorm_chol) +export(dmnorm_inv_ld) export(dmulti) export(dmvt_chol) export(dsqrtinvgamma) @@ -193,6 +195,7 @@ export(optimDefaultControl) export(optimResultNimbleList) export(parameterTransform) export(pdexp) +export(PDinverse_logdet) export(pexp_nimble) export(phi) export(pinvgamma) @@ -225,6 +228,7 @@ export(rinvgamma) export(rinvwish_chol) export(rlkj_corr_cholesky) export(rmnorm_chol) +export(rmnorm_inv_ld) export(rmulti) export(rmvt_chol) export(rsqrtinvgamma) diff --git a/packages/nimble/R/BUGS_BUGSdecl.R b/packages/nimble/R/BUGS_BUGSdecl.R index 1328a0b4c..65105c01f 100644 --- a/packages/nimble/R/BUGS_BUGSdecl.R +++ b/packages/nimble/R/BUGS_BUGSdecl.R @@ -77,6 +77,7 @@ nimblePreevaluationFunctionNames <- c('+', 'asCol', 'logdet', 'chol', + 'PDinverse_logdet', 'inverse', 'forwardsolve', 'backsolve', diff --git a/packages/nimble/R/BUGS_modelDef.R b/packages/nimble/R/BUGS_modelDef.R index 78ece52b2..c8285c040 100644 --- a/packages/nimble/R/BUGS_modelDef.R +++ b/packages/nimble/R/BUGS_modelDef.R @@ -831,6 +831,17 @@ modelDefClass$methods(reparameterizeDists = function() { BUGSdecl <- declInfo[[i]] ## grab this current BUGS declation info object if(BUGSdecl$type == 'determ') next ## skip deterministic nodes code <- BUGSdecl$code ## grab the original code + if(BUGSdecl$distributionName == "dmnorm" && buildDerivs && getNimbleOption('useADdmnorm')) { + if(length(BUGSdecl$code) > 2 && "cholesky" %in% names(BUGSdecl$code[[3]])) { + messageIfVerbose(" [Note] Detected use of `cholesky` parameterization of `dmnorm` with a\n", + " derivative-enabled model. AD-optimized `dmnorm` is only available\n", + " for the `prec` or `cov` parameterizations. NIMBLE will use a version\n", + " of `dmnorm` not optimized for AD, which may result in inefficiency.") + } else { + BUGSdecl$distributionName <- "dmnormAD" + BUGSdecl$valueExpr[[1]] <- quote(dmnormAD) + } + } valueExpr <- BUGSdecl$valueExpr ## grab the RHS (distribution) distName <- BUGSdecl$distributionName #as.character(valueExpr[[1]]) if(!(distName %in% getAllDistributionsInfo('namesVector'))) stop('unknown distribution name: ', distName) ## error if the distribution isn't something we recognize @@ -1059,6 +1070,16 @@ liftedCallsGetIndexingFromArgumentNumbers <- list( CAR_calcEVs3 = c(3) ) +liftedCallsGetIndexingOther <- list( + ## This is general in that it finds the number of elements of the matrix, + ## but the input shouldn't be anything other than square. + PDinverse_logdet = function(argList) { + getlen <- function(arg) length(eval(arg)) + list(substitute(1:N, list(N = prod(sapply(argList[[1]][3:length(argList[[1]])], getlen))+1))) + } +) + + modelDefClass$methods(liftExpressionArgs = function() { ## overwrites declInfo (*and adds*), lifts any expressions in distribution arguments to new nodes newDeclInfo <- list() @@ -1125,6 +1146,7 @@ isExprLiftable <- function(paramExpr, type = NULL) { callText <- getCallText(paramExpr) if(callText == 'chol') return(TRUE) ## do lift calls to chol(...) if(callText == 'inverse') return(TRUE) ## do lift calls to inverse(...) + if(callText == 'PDinverse_logdet') return(TRUE) ## do lift calls to PDinverse_logdet(...) if(callText == 'CAR_calcNumIslands') return(TRUE) ## do lift calls to CAR_calcNumIslands(...) if(callText == 'CAR_calcC') return(TRUE) ## do lift calls to CAR_calcC(...) if(callText == 'CAR_calcM' ) return(TRUE) ## do lift calls to CAR_calcM(...) @@ -1148,6 +1170,8 @@ isExprLiftable <- function(paramExpr, type = NULL) { addNecessaryIndexingToNewNode <- function(newNodeNameExpr, paramExpr, indexVarExprs) { if(is.call(paramExpr) && safeDeparse(paramExpr[[1]], warn = TRUE) %in% names(liftedCallsGetIndexingFromArgumentNumbers)) return(addNecessaryIndexingFromArgumentNumbers(newNodeNameExpr, paramExpr, indexVarExprs)) + if(is.call(paramExpr) && safeDeparse(paramExpr[[1]], warn = TRUE) %in% names(liftedCallsGetIndexingOther)) + return(addNecessaryIndexingOther(newNodeNameExpr, paramExpr, indexVarExprs)) usedIndexVarsList <- indexVarExprs[indexVarExprs %in% all.vars(paramExpr)] # this extracts any index variables which appear in 'paramExpr' vectorizedIndexExprsList <- extractAnyVectorizedIndexExprs(paramExpr) # creates a list of any vectorized (:) indexing expressions appearing in 'paramExpr' neededIndexExprsList <- c(usedIndexVarsList, vectorizedIndexExprsList) @@ -1165,6 +1189,15 @@ addNecessaryIndexingFromArgumentNumbers <- function(newNodeNameExpr, paramExpr, newNodeNameExprIndexed[3:(2+length(neededIndexExprsList))] <- neededIndexExprsList return(newNodeNameExprIndexed) } +addNecessaryIndexingOther <- function(newNodeNameExpr, paramExpr, indexVarExprs) { + paramExprCallName <- as.character(paramExpr[[1]]) + neededIndexExprsList <- liftedCallsGetIndexingOther[[paramExprCallName]](as.list(paramExpr[-1])) + newNodeNameExprIndexed <- substitute(NAME[], list(NAME = newNodeNameExpr)) + newNodeNameExprIndexed[3:(2+length(neededIndexExprsList))] <- neededIndexExprsList + return(newNodeNameExprIndexed) +} + + extractAnyVectorizedIndexExprs <- function(expr) { if(!(':' %in% all.names(expr))) return(list()) if(!is.call(expr)) return(list()) diff --git a/packages/nimble/R/MCMC_conjugacy.R b/packages/nimble/R/MCMC_conjugacy.R index a4a22753d..991e8a58b 100644 --- a/packages/nimble/R/MCMC_conjugacy.R +++ b/packages/nimble/R/MCMC_conjugacy.R @@ -133,7 +133,8 @@ conjugacyRelationshipsInputList <- list( link = 'linear', dependents = list( ##dmnorm = list(param = 'mean', contribution_mean = '(t(coeff) %*% prec %*% asCol(value-offset))[,1]', contribution_prec = 't(coeff) %*% prec %*% coeff')), - dmnorm = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)')), + dmnorm = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)'), + dmnormAD = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)')), ## LINK will be replaced with appropriate link via code processing ## original less efficient posterior definition: ## posterior = 'dmnorm_chol(mean = (inverse(prior_prec + contribution_prec) %*% (prior_prec %*% asCol(prior_mean) + asCol(contribution_mean)))[,1], @@ -144,6 +145,21 @@ conjugacyRelationshipsInputList <- list( mu <- backsolve(R, forwardsolve(t(R), A))[,1] dmnorm_chol(mean = mu, cholesky = R, prec_param = 1) }'), + list(prior = 'dmnormAD', + link = 'linear', + dependents = list( + ##dmnorm = list(param = 'mean', contribution_mean = '(t(coeff) %*% prec %*% asCol(value-offset))[,1]', contribution_prec = 't(coeff) %*% prec %*% coeff')), + dmnorm = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)'), + dmnormAD = list(param = 'mean', contribution_mean = '(calc_dmnormConjugacyContributions(coeff, prec, value-offset, 1, 0))[,1]', contribution_prec = 'calc_dmnormConjugacyContributions(coeff, prec, value-offset, 2, 0)')), + ## LINK will be replaced with appropriate link via code processing + ## original less efficient posterior definition: + ## posterior = 'dmnorm_chol(mean = (inverse(prior_prec + contribution_prec) %*% (prior_prec %*% asCol(prior_mean) + asCol(contribution_mean)))[,1], + ## cholesky = chol(prior_prec + contribution_prec), + ## prec_param = 1)'), + posterior = '{ R <- chol(prior_prec + contribution_prec) + A <- prior_prec %*% asCol(prior_mean) + asCol(contribution_mean) + mu <- backsolve(R, forwardsolve(t(R), A))[,1] + dmnorm_chol(mean = mu, cholesky = R, prec_param = 1) }'), ## wishart list(prior = 'dwish', @@ -159,7 +175,8 @@ conjugacyRelationshipsInputList <- list( ## changing to only use link='identity' case, since the link='linear' case was not correct ## -DT March 2017 ## dmnorm = list(param = 'prec', contribution_R = 'asCol(value-mean) %*% (asRow(value-mean) %*% coeff)', contribution_df = '1')), - dmnorm = list(param = 'prec', contribution_R = 'coeff * asCol(value-mean) %*% asRow(value-mean)', contribution_df = '1')), + dmnorm = list(param = 'prec', contribution_R = 'coeff * asCol(value-mean) %*% asRow(value-mean)', contribution_df = '1'), + dmnormAD = list(param = 'prec', contribution_R = 'coeff * asCol(value-mean) %*% asRow(value-mean)', contribution_df = '1')), posterior = 'dwish_chol(cholesky = chol(prior_R + contribution_R), df = prior_df + contribution_df, scale_param = 0)'), @@ -168,7 +185,8 @@ conjugacyRelationshipsInputList <- list( list(prior = 'dinvwish', link = 'multiplicativeScalar', # we only handle scalar 'coeff'; this naming is slightly awkward since for univar dists, link is of course scalar dependents = list( - dmnorm = list(param = 'cov', contribution_S = 'asCol(value-mean) %*% asRow(value-mean) / coeff', contribution_df = '1')), + dmnorm = list(param = 'cov', contribution_S = 'asCol(value-mean) %*% asRow(value-mean) / coeff', contribution_df = '1'), + dmnormAD = list(param = 'cov', contribution_S = 'asCol(value-mean) %*% asRow(value-mean) / coeff', contribution_df = '1')), posterior = 'dinvwish_chol(cholesky = chol(prior_S + contribution_S), df = prior_df + contribution_df, scale_param = 1)') @@ -999,7 +1017,7 @@ conjugacyClass <- setRefClass( for(contributionName in posteriorObject$neededContributionNames) { if(!(contributionName %in% dependents[[distName]]$contributionNames)) next contributionExpr <- dependents[[distName]]$contributionExprs[[contributionName]] - if(distName == 'dmnorm' && prior == 'dmnorm') { + if(distName %in% c('dmnorm','dmnormAD') && prior %in% c('dmnorm','dmnormAD')) { ## need to deal with [,1] in contribution_mean if(contributionName == 'contribution_mean') tmpExpr <- contributionExpr[[2]][[2]] else tmpExpr <- contributionExpr if(getNimbleOption('allowDynamicIndexing') && doDependentScreen) { @@ -1348,8 +1366,11 @@ cc_otherParamsCheck <- function(model, depNode, targetNode, skipExpansionsNode = if(!missing(depParamNodeName) && (names(paramsList)[i] == depParamNodeName)) { expr <- depNodeExprExpanded } else { expr <- cc_expandDetermNodesInExpr(model, paramsList[[i]], targetNode, skipExpansionsNode) } - if(cc_vectorizedComponentCheck(targetNode, expr)) return(FALSE) - if(cc_nodeInExpr(targetNode, expr)) { timesFound <- timesFound + 1 } ## we found 'targetNode' + ## We expect to find target in PDinverse_logdet() one extra time when dmnormAD used. + if(!(length(expr) > 1 && expr[[1]] == "PDinverse_logdet")) { + if(cc_vectorizedComponentCheck(targetNode, expr)) return(FALSE) + if(cc_nodeInExpr(targetNode, expr)) { timesFound <- timesFound + 1 } ## we found 'targetNode' + } } if(timesFound == 0) stop('something went wrong; targetNode not found in any parameter expressions') if(timesFound == 1) return(TRUE) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 9a548ca03..ee5410cb2 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -1170,7 +1170,7 @@ essNF_multivariate <- nimbleFunction( name = 'essNF_multivariate', contains = essNFList_virtual, setup = function(model, node) { - if(!(model$getDistribution(node) == 'dmnorm')) stop('something went wrong') + if(!(model$getDistribution(node) %in% c('dmnorm', 'dmnormAD'))) stop('sampler_ess: node `', node, '` does not have a dmnorm distribution') }, run = function() { mean <- model$getParam(node, 'mean') @@ -1201,10 +1201,10 @@ sampler_ess <- nimbleFunction( ## nested function and function list definitions essNFList <- nimbleFunctionList(essNFList_virtual) if(model$getDistribution(target) == 'dnorm') essNFList[[1]] <- essNF_univariate(model, target) - if(model$getDistribution(target) == 'dmnorm') essNFList[[1]] <- essNF_multivariate(model, target) + if(model$getDistribution(target) %in% c('dmnorm','dmnormAD')) essNFList[[1]] <- essNF_multivariate(model, target) ## checks if(length(target) > 1) stop('elliptical slice sampler only applies to one target node') - if(!(model$getDistribution(target) %in% c('dnorm', 'dmnorm'))) stop('elliptical slice sampler only applies to normal distributions') + if(!(model$getDistribution(target) %in% c('dnorm', 'dmnorm', 'dmnormAD'))) stop('elliptical slice sampler only applies to normal distributions') targetNames <- createNamesString(target) }, run = function() { @@ -2862,7 +2862,7 @@ sampler_polyagamma <- nimbleFunction( ## Conjugacy checking, part 1. if(check) { - if(!all(targetDists %in% c("dnorm", "dmnorm"))) + if(!all(targetDists %in% c("dnorm", "dmnorm", "dmnormAD"))) stop("polyagamma sampler: all target nodes must have `dnorm` or `dmnorm` priors. ", checkMessage) if(!all(model$getDistribution(yNodes) %in% c("dbern", "dbin")) ) stop("polyagamma sampler: response nodes must be distributed `dbern` or `dbin`. ", checkMessage) @@ -2981,7 +2981,7 @@ sampler_polyagamma <- nimbleFunction( singleSize <- FALSE dnormNodes <- targetDists == "dnorm" - dmnormNodes <- targetDists == "dmnorm" + dmnormNodes <- targetDists %in% c("dmnorm", "dmnormAD") n_dnorm <- sum(dnormNodes) n_dmnorm <- sum(dmnormNodes) diff --git a/packages/nimble/R/cppDefs_ADtools.R b/packages/nimble/R/cppDefs_ADtools.R index 6d6c4e367..1392efe5b 100644 --- a/packages/nimble/R/cppDefs_ADtools.R +++ b/packages/nimble/R/cppDefs_ADtools.R @@ -198,6 +198,20 @@ make_deriv_function <- function(origFun, nDim = 1, type = 'double', ref = TRUE, const = TRUE)) newFun$args$addSymbol(cppVar(baseType = 'bool', name = "DO_UPDATE_")) newFun$args$addSymbol(cppVar(baseType = 'bool', name = "RESET_")) + newFun$args$addSymbol(cppNimArr(name = 'ARGZ_outInds__', + nDim = 1, type = 'double', ref = TRUE, const = TRUE)) + + inDirSym <- cppNimArr(name = 'ARGZ_inDir__', + nDim = 1, type = 'double', ref = TRUE, const = TRUE) + outDirSym <- cppNimArr(name = 'ARGZ_outDir__', + nDim = 1, type = 'double', ref = TRUE, const = TRUE) + if(meta) { + inDirSym <- cppVarSym2templateTypeCppVarSym(inDirSym) + outDirSym <- cppVarSym2templateTypeCppVarSym(outDirSym) + } + newFun$args$addSymbol(inDirSym) + newFun$args$addSymbol(outDirSym) + newFun$args$addSymbol(cppVar(name = 'ARGZ_ADinfo_', ref = TRUE, baseType = "nimbleCppADinfoClass")) @@ -233,6 +247,9 @@ make_deriv_function <- function(origFun, getDerivsRcall <- substitute(returnList_ <- GETDERIVS_WRAPPER( INNERCALL, ARGZ_nimDerivsOrders_, ARGZ_wrtVector_ , + ARGZ_outInds__, + ARGZ_inDir__, + ARGZ_outDir__, recordingInfo_), list(INNERCALL = innerRcall, GETDERIVS_WRAPPER = as.name(getDerivs_wrapper))) @@ -240,7 +257,10 @@ make_deriv_function <- function(origFun, getDerivs_wrapper <- 'getDerivs_wrapper' getDerivsRcall <- substitute(returnList_ <- GETDERIVS_WRAPPER( INNERCALL, ARGZ_nimDerivsOrders_, - ARGZ_wrtVector_ ), + ARGZ_wrtVector_ , + ARGZ_outInds__, + ARGZ_inDir__, + ARGZ_outDir__), list(INNERCALL = innerRcall, GETDERIVS_WRAPPER = as.name(getDerivs_wrapper))) } diff --git a/packages/nimble/R/cppDefs_nimbleFunction.R b/packages/nimble/R/cppDefs_nimbleFunction.R index 73b228eaf..32a13d9f0 100644 --- a/packages/nimble/R/cppDefs_nimbleFunction.R +++ b/packages/nimble/R/cppDefs_nimbleFunction.R @@ -52,7 +52,7 @@ cppNimbleClassClass <- setRefClass('cppNimbleClassClass', loaded = 'ANY', Cwritten = 'ANY', RCfunDefs = 'ANY' - + ), methods = list( getDefs = function() { @@ -183,7 +183,7 @@ cppNimbleClassClass <- setRefClass('cppNimbleClassClass', pointerLine <- list() if(!(is.null(nimCompProc[['nimbleListObj']]))){ flagText <- paste0('RCopiedFlag = false') - flagLine[[1]] <- substitute(FLAGTEXT, + flagLine[[1]] <- substitute(FLAGTEXT, list(FLAGTEXT = as.name(flagText))) pointerText <- paste0('RObjectPointer = NULL') pointerLine[[1]] <- substitute(POINTERTEXT, @@ -214,7 +214,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', ), methods = list( getDefs = function() { - c(callSuper(), SEXPmemberInterfaceFuns) + c(callSuper(), SEXPmemberInterfaceFuns) }, getHincludes = function() { c(callSuper(), unlist(lapply(SEXPmemberInterfaceFuns, function(x) x$getHincludes()), recursive = FALSE)) @@ -236,7 +236,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', inheritance <<- c(inheritance, 'nodeFun') parentsSizeAndDims <<- environment(nfProc$nfGenerator)$parentsSizeAndDims ADconstantsInfo <<- environment(nfProc$nfGenerator)$ADconstantsInfo - + } } }, @@ -260,7 +260,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', ) SEXPmemberInterfaceFuns <<- SEXPmemberInterfaceFuns[ !unlist(lapply(SEXPmemberInterfaceFuns, is.null)) ] - + nimCompProc <<- nfProc }, buildFunctionDefs = function() { @@ -289,7 +289,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', ## functionDefs[[newFunName]] <<- res$fun ## useModelInfo$nodeFxnVector_name <- res[['nodeFxnVector_name']] ## } else { - ## We need a separate one for the "reconfigure" system that is non-static + ## We need a separate one for the "reconfigure" system that is non-static ## if(isTRUE(nimbleOptions('useADreconfigure'))) { newFunName <- paste0(funName, '_AD2_') res <- makeTypeTemplateFunction(newName = newFunName, @@ -299,7 +299,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', derivControl = derivControl) functionDefs[[newFunName]] <<- res$fun useModelInfo$nodeFxnVector_name <- res[['nodeFxnVector_name']] - + ## } ## } useModelInfo @@ -323,7 +323,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', ## Add meta-taping type template function of the new function newMetaFunName <- paste0(funName, '_AD2__deriv_') argTransName <- paste0(funName, '_ADargumentTransfer2__AD2_') - + functionDefs[[newMetaFunName]] <<- make_deriv_function(regularFun, newMetaFunName, independentVarNames, @@ -475,7 +475,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', add_deriv_function(funName, independentVarNames, derivControl) - + funIndex + 1 ## function return value increments by one in non-meta case } else { funIndex ## but does not increment in meta case @@ -531,7 +531,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', addInheritance("nimbleFunctionCppADbase") addADinfoObjects(.self) } - constructorCode + constructorCode }, buildCmultiInterface = function(dll = NULL) { sym <- if(!is.null(dll)) @@ -577,7 +577,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', }, buildAll = function(where = where) { baseClassObj <- environment(nfProc$nfGenerator)$contains - + if(!is.null(baseClassObj)) { inheritance <<- inheritance[inheritance != 'NamedObjects'] baseClassName <- environment(baseClassObj)$name @@ -592,9 +592,9 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', updateADproxyModelMethods(.self) } } - + addCopyFromRobject() - + callSuper(where) if(handleDerivs) { if(!is.null(constructorCode)) { ## insert AD-related code to constructor @@ -645,6 +645,7 @@ cppNimbleFunctionClass <- setRefClass('cppNimbleFunctionClass', ## The next block of code has the initial setup for an AST processing stage ## to make modifications for AD based on context etc. modifyForAD_handlers <- c(list( + nimDerivs_dummy = 'modifyForAD_nimDerivs_dummy', `[` = 'modifyForAD_indexingBracket', # pow = 'modifyForAD_issuePowWarning', eigenBlock = 'modifyForAD_eigenBlock', @@ -715,7 +716,7 @@ exprClasses_modifyForAD <- function(code, symTab, workEnv[['nodeFxnVector_name']], ", recordingInfo_.tape_id(), recordingInfo_.tape_handle());"))) ## ", TYPE_::get_tape_id_nimble(), TYPE_::get_tape_handle_nimble());"))) - + set_atomic_info_line <- RparseTree2ExprClasses(cppLiteral( paste0("set_CppAD_atomic_info_for_model(", @@ -734,7 +735,7 @@ exprClasses_modifyForAD <- function(code, symTab, setArg(code, 1, newExpr) } } - } + } invisible(NULL) } @@ -758,6 +759,23 @@ recurse_modifyForAD <- function(code, symTab, workEnv) { invisible(NULL) } +modifyForAD_nimDerivs_dummy <- function(code, symTab, workEnv) { + # nimDerivs_dummy is a placeholder in the AST that does not get output in the genCpp step. + # In this handler, we need to look for make_vector_if_needed for numeric literals for the inDir and outDir arguments + # and wrap them in TYPE_. + for(argName in c('inDir', 'outDir')) { + if(!is.null(code$args[[argName]])) { + arg <- code$args[[argName]] + if(!inherits(arg, 'exprClass')) next + if(arg$isCall && arg$name == 'make_vector_if_necessary') { + if(length(arg$args) > 0 && is.numeric(arg$args[[1]])) + insertExprClassLayer(arg, 1, 'TYPE_') + } + } + } + invisible(NULL) +} + modifyForAD_indexingBracket <- function(code, symTab, workEnv) { if(!isTRUE(workEnv$.anyAD)) return(invisible(NULL)) if(is.null(code$type)) return(invisible(NULL)) # arises from constructions like setMap that lack type annotation @@ -821,7 +839,7 @@ modifyForAD_AssignEigenMap <- function(code, symTab, workEnv) { ## We extract the variable name by undoing the construction of the Eigen Map variable name. ## This would be nicer to do with generic functions that undo the label creation. ## At this time of this writing, hacking the solution here was more practical. - ## + ## ## This method of extracting the name of the variable we are making a map into ## attempted to invert the name-generation and name-mangling steps that led to it. ## However, it is not fully general. In addition, it is not clear if this would @@ -1135,7 +1153,7 @@ makeSingleCopyCall <- function(varName, cppCopyType, buildDerivs = FALSE) { switch(cppCopyType, 'nimbleFunction' = { - cppLiteral(paste0("COPY_NIMBLE_FXN_FROM_R_OBJECT(\"", varName, "\");")) + cppLiteral(paste0("COPY_NIMBLE_FXN_FROM_R_OBJECT(\"", varName, "\");")) }, 'numericVector' = { cppLiteral(paste0("COPY_NUMERIC_VECTOR_FROM_R_OBJECT(\"", varName, "\");")) @@ -1200,7 +1218,7 @@ makeCopyFromRobjectDef <- function(cppCopyTypes, objectDefs = localVars) copyFromRobjectInterfaceDef <- NULL - + list(copyFromRobjectDef = copyFromRobjectDef, copyFromRobjectInterfaceDef = copyFromRobjectInterfaceDef) } diff --git a/packages/nimble/R/distributions_implementations.R b/packages/nimble/R/distributions_implementations.R index 1b3a6a7e3..78da6429f 100644 --- a/packages/nimble/R/distributions_implementations.R +++ b/packages/nimble/R/distributions_implementations.R @@ -545,6 +545,36 @@ qdexp <- function(p, location = 0, scale = 1, rate = 1/scale, lower.tail = TRUE, #' dmnorm_chol(x, mean, ch, prec_param = FALSE) NULL +PDinverse_logdet <- function(mat) { + if(storage.mode(mat) != 'double') + storage.mode(mat) <- 'double' + out <- .Call(C_PDinverse_logdet, mat) + return(out) +} + +#' @rdname MultivariateNormal +#' @export +dmnorm_inv_ld <- function(x, mean, mat, inv_ld, prec_param = TRUE, log = FALSE) { + if(storage.mode(inv_ld) != 'double') + storage.mode(inv_ld) <- 'double' + if(storage.mode(mat) != 'double') + storage.mode(mat) <- 'double' + .Call(C_dmnorm_inv_ld, as.double(x), as.double(mean), mat, inv_ld, + as.double(prec_param), as.logical(log)) +} + +#' @rdname MultivariateNormal +#' @export +rmnorm_inv_ld <- function(n=1, mean, mat, inv_ld, prec_param = TRUE) { + if(n != 1) warning('rmnorm_inv_ld only handles n = 1 at the moment') + if(storage.mode(inv_ld) != 'double') + storage.mode(inv_ld) <- 'double' + if(storage.mode(mat) != 'double') + storage.mode(mat) <- 'double' + .Call(C_rmnorm_inv_ld, as.double(mean), mat, inv_ld, + as.double(prec_param)) +} + #' @rdname MultivariateNormal #' @export dmnorm_chol <- function(x, mean, cholesky, prec_param = TRUE, log = FALSE) { diff --git a/packages/nimble/R/distributions_inputList.R b/packages/nimble/R/distributions_inputList.R index a8d1926de..0de248cfc 100644 --- a/packages/nimble/R/distributions_inputList.R +++ b/packages/nimble/R/distributions_inputList.R @@ -210,6 +210,15 @@ distributionsInputList <- list( 'cov = calc_dmnormAltParams(cholesky, prec_param, 0)'), types = c('value = double(1)', 'mean = double(1)', 'cholesky = double(2)', 'prec = double(2)', 'cov = double(2)')), + dmnormAD = list(BUGSdist = 'dmnormAD(mean, prec, cov, inv_ld, prec_param)', + Rdist = c('dmnorm_inv_ld(mean, mat = cov, inv_ld = PDinverse_logdet(cov), prec_param = 0)', + 'dmnorm_inv_ld(mean, mat = prec, inv_ld = PDinverse_logdet(prec), prec_param = 1)'), + altParams= c('prec = calc_dmnorm_inv_ld_AltParams(mat, inv_ld, prec_param, 1)', + 'cov = calc_dmnorm_inv_ld_AltParams(mat, inv_ld, prec_param, 0)'), + mixedSizes = TRUE, + types = c('value = double(1)', 'mean = double(1)', 'inv_ld = double(1)', + 'prec = double(2)', 'cov = double(2)', 'mat = double(2)')), + dmulti = list(BUGSdist = 'dmulti(prob, size)', Rdist = 'dmulti(size, prob)', types = c('value = double(1)', 'prob = double(1)'), diff --git a/packages/nimble/R/genCpp_eigenization.R b/packages/nimble/R/genCpp_eigenization.R index 636cad2d8..07cd35df4 100644 --- a/packages/nimble/R/genCpp_eigenization.R +++ b/packages/nimble/R/genCpp_eigenization.R @@ -110,7 +110,7 @@ eigenizeCalls <- c( ## component-wise unarys valid for either Eigen array or mat ) eigenizeCallsBeforeRecursing <- c( ## These cannot be calls that trigger aliasRisk. getParam always triggers an intermediate so it should never need handling here - makeCallList(c('size','nimArr_dmnorm_chol', 'nimArr_dmvt_chol', 'nimArr_dlkj_corr_cholesky', 'nimArr_dwish_chol', 'nimArr_dinvwish_chol', 'nimArr_dcar_normal', 'nimArr_dcar_proper', 'nimArr_ddirch','calculate','calculateDiff','getLogProb', 'getParam', 'getBound', 'getNodeFunctionIndexedInfo', 'concatenateTemp', 'MAKE_FIXED_VECTOR', 'hardCodedVectorInitializer'), 'eigenize_doNotRecurse'), + makeCallList(c('size','nimArr_dmnorm_chol', 'nimArr_dmnorm_inv_ld', 'nimArr_dmvt_chol', 'nimArr_dlkj_corr_cholesky', 'nimArr_dwish_chol', 'nimArr_dinvwish_chol', 'nimArr_dcar_normal', 'nimArr_dcar_proper', 'nimArr_ddirch','calculate','calculateDiff','getLogProb', 'getParam', 'getBound', 'getNodeFunctionIndexedInfo', 'concatenateTemp', 'MAKE_FIXED_VECTOR', 'hardCodedVectorInitializer'), 'eigenize_doNotRecurse'), list(coeffSetter = 'eigenize_coeffSetter', nfVar = 'eigenize_nfVar', chainedCall = 'eigenize_chainedCall'), diff --git a/packages/nimble/R/genCpp_operatorLists.R b/packages/nimble/R/genCpp_operatorLists.R index 5f37d85e0..194101559 100644 --- a/packages/nimble/R/genCpp_operatorLists.R +++ b/packages/nimble/R/genCpp_operatorLists.R @@ -194,6 +194,7 @@ callToSkipInEigenization <- c('copy', 'blank', 'rankSample', 'nimArr_dmnorm_chol', + 'nimArr_dmnorm_inv_ld', 'nimArr_dmvt_chol', 'nimArr_dlkj_corr_cholesky', 'nimArr_dwish_chol', @@ -205,6 +206,7 @@ callToSkipInEigenization <- c('copy', 'nimArr_dinterval', 'nimArr_ddirch', 'nimArr_rmnorm_chol', + 'nimArr_rmnorm_inv_ld', 'nimArr_rmvt_chol', 'nimArr_rlkj_corr_cholesky', 'nimArr_rwish_chol', @@ -340,7 +342,8 @@ nimDerivsPrependTypeOperators <- c("dnorm", "dpois", "dgamma", "dinvgamma", "dsq "dexp_nimble", "dexp", "ddexp", "dlnorm", "dweibull", "dbinom", "dbeta", "dchisq", "dlogis", "dt", "dt_nonstandard", "nimArr_dmulti", "nimArr_dcat", "dnbinom", "dunif", "pairmax", "pairmin", - "nimArr_ddirch", "nimArr_dmvt_chol", "nimArr_dmnorm_chol", + "nimArr_ddirch", "nimArr_dmvt_chol", "nimArr_dmnorm_chol", "nimArr_dmnorm_inv_ld", + "PDinverse_logdet", "nimArr_dwish_chol", "nimArr_dinvwish_chol", "nimArr_dlkj_corr_cholesky", "nimArr_dcar_normal", "nimArr_dcar_proper", "nimStep", 'ilogit', 'icloglog', 'iprobit', 'probit', 'cloglog', diff --git a/packages/nimble/R/genCpp_processSpecificCalls.R b/packages/nimble/R/genCpp_processSpecificCalls.R index 1d01a8c52..911d5ccb3 100644 --- a/packages/nimble/R/genCpp_processSpecificCalls.R +++ b/packages/nimble/R/genCpp_processSpecificCalls.R @@ -41,9 +41,10 @@ specificCallHandlers = c( '^' = "powHandler"), makeCallList(names(specificCallReplacements), 'replacementHandler'), makeCallList(c('nimNumeric', 'nimLogical', 'nimInteger', 'nimMatrix', 'nimArray'), 'nimArrayGeneralHandler' ), - makeCallList(c('dmnorm_chol', 'dmvt_chol', 'dlkj_corr_cholesky', 'dwish_chol', 'dinvwish_chol', 'dcar_normal', 'dcar_proper', 'dmulti', 'dcat', 'dinterval', 'ddirch'), 'dmFunHandler') + makeCallList(c('dmnorm_chol', 'dmnorm_inv_ld', 'dmvt_chol', 'dlkj_corr_cholesky', 'dwish_chol', 'dinvwish_chol', 'dcar_normal', 'dcar_proper', 'dmulti', 'dcat', 'dinterval', 'ddirch'), 'dmFunHandler') ) specificCallHandlers[['rmnorm_chol']] <- 'rmFunHandler' +specificCallHandlers[['rmnorm_inv_ld']] <- 'rmFunHandler' specificCallHandlers[['rmvt_chol']] <- 'rmFunHandler' specificCallHandlers[['rlkj_corr_cholesky']] <- 'rmFunHandler' specificCallHandlers[['rwish_chol']] <- 'rmFunHandler' @@ -262,7 +263,7 @@ replacementHandler <- function(code, symTab) { } dmFunHandler <- function(code, symTab) { - if(code$name %in% c('dwish_chol', 'dinvwish_chol', 'dmnorm_chol', 'dmvt_chol')) + if(code$name %in% c('dwish_chol', 'dinvwish_chol', 'dmnorm_chol', 'dmnorm_inv_ld', 'dmvt_chol')) code$args$overwrite_inputs <- 0 code$name <- paste0('nimArr_', code$name) } diff --git a/packages/nimble/R/genCpp_sizeProcessing.R b/packages/nimble/R/genCpp_sizeProcessing.R index ba1cf9746..5c2870989 100644 --- a/packages/nimble/R/genCpp_sizeProcessing.R +++ b/packages/nimble/R/genCpp_sizeProcessing.R @@ -6,6 +6,7 @@ sizeProc_storage_mode <- function(x) { } assignmentAsFirstArgFuns <- c('nimArr_rmnorm_chol', + 'nimArr_rmnorm_inv_ld', 'nimArr_rmvt_chol', 'nimArr_rlkj_corr_cholesky', 'nimArr_rwish_chol', @@ -101,7 +102,8 @@ sizeCalls <- c( Rf_eval = 'sizeReval', nimbleConvert = 'sizeNimbleConvert', nimbleUnconvert = 'sizeNimbleUnconvert', - asReturnSymbol = 'sizeAsReturnSymbol'), + asReturnSymbol = 'sizeAsReturnSymbol', + PDinverse_logdet = 'sizePDinverse_logdet'), makeCallList(scalar_distribution_dFuns, 'sizeRecyclingRule'), makeCallList(scalar_distribution_pFuns, 'sizeRecyclingRule'), makeCallList(scalar_distribution_qFuns, 'sizeRecyclingRule'), @@ -119,6 +121,7 @@ sizeCalls <- c( rexp = 'sizeRecyclingRuleRfunction', makeCallList(c('nimAnyNA','nimAnyNaN'), 'sizeScalarRecurse'), makeCallList(c('nimArr_dmnorm_chol', + 'nimArr_dmnorm_inv_ld', 'nimArr_dmvt_chol', 'nimArr_dlkj_corr_cholesky', 'nimArr_dwish_chol', @@ -130,6 +133,7 @@ sizeCalls <- c( 'nimArr_dinterval', 'nimArr_ddirch'), 'sizeScalarRecurseAllowMaps'), makeCallList(c('nimArr_rmnorm_chol', + 'nimArr_rmnorm_inv_ld', 'nimArr_rmvt_chol', 'nimArr_rlkj_corr_cholesky', 'nimArr_rwish_chol', @@ -1279,11 +1283,13 @@ sizeNimDerivs <- function(code, symTab, typeEnv){ ## asserts <- sizeNimbleListReturningFunction(code, symTab, typeEnv) ## lift wrt if needed. I'm not sure why sizeNimbleListReturningFunction doesn't handle lifting - if(inherits(code$args[['wrt']], 'exprClass')) { - if(!code$args[['wrt']]$isName) { - iWrt <- which(names(code$args) == 'wrt') - if(length(iWrt) != 1) stop("problem working on wrt argument to nimDerivs") - asserts <- c(asserts, sizeInsertIntermediate(code, iWrt, symTab, typeEnv) ) + for(liftArgName in c("wrt", "outInds", "inDir", "outDir")) { + if(inherits(code$args[[liftArgName]], 'exprClass')) { + if(!code$args[[liftArgName]]$isName) { + iLift <- which(names(code$args) == liftArgName) + if(length(iLift) != 1) stop("problem working on ", liftArgName, " argument to nimDerivs") + asserts <- c(asserts, sizeInsertIntermediate(code, iLift, symTab, typeEnv)) + } } } if(inherits(code$args[['order']], 'exprClass')) { @@ -1302,16 +1308,23 @@ sizeNimDerivs <- function(code, symTab, typeEnv){ ## asserts <- c(asserts, sizeInsertIntermediate(code, iOrder, symTab, typeEnv) ) } } + for(makeVecName in c('wrt', 'order', 'outInds', 'inDir', 'outDir')) { + iMakeVec <- which(names(code$args) == makeVecName) + if(length(iMakeVec) != 1) stop("problem working on ", makeVecName, " argument to nimDerivs") + insertExprClassLayer(code, iMakeVec, 'make_vector_if_necessary', + type = 'double', + nDim = 1, + sizeExprs = list()) + } +# insertExprClassLayer(code, which(names(code$args)=='wrt'), 'make_vector_if_necessary', +# type = 'double', +# nDim = 1, +# sizeExprs = list()) - insertExprClassLayer(code, which(names(code$args)=='wrt'), 'make_vector_if_necessary', - type = 'double', - nDim = 1, - sizeExprs = list()) - - a1 <- insertExprClassLayer(code, which(names(code$args)=='order'), 'make_vector_if_necessary', - type = 'double', - nDim = 1, - sizeExprs = list()) +# a1 <- insertExprClassLayer(code, which(names(code$args)=='order'), 'make_vector_if_necessary', +# type = 'double', +# nDim = 1, +# sizeExprs = list()) newADinfoName <- ADinfoLabel() ## symTab$addSymbol(symbolADinfo$new(name = newADinfoName)) if(!is.list(code$aux)) @@ -3285,6 +3298,65 @@ sizeMatrixSquareReduction <- function(code, symTab, typeEnv) { if(length(asserts) == 0) NULL else asserts } +sizePDinverse_logdet <- function(code, symTab, typeEnv) { + # Modeled after a combination of generalFunSize and sizeUnaryCwiseSquare + if(length(code$args) != 1) { + stop(exprClassProcessingErrorMsg(code, 'sizePDinverse_logdet called with argument length != 1.'), call. = FALSE) + } + a1 <- code$args[[1]] + if(!inherits(a1, 'exprClass')) + stop(exprClassProcessingErrorMsg(code, 'sizePDinverse_logdet called with argument that is not an expression.'), call. = FALSE) + ## Ensure that simple maps being passed will be passed without extra + ## copy that would occur from lifting an Eigen expression. + ## (Relevant primarily for first argument, but no harm in applying to second.) + for(i in seq_along(code$args)) { + if(inherits(code$args[[i]], 'exprClass')) { + if(code$args[[i]]$name == "[") { + if(inherits(code$args[[i]]$args[[1]], + 'exprClass')) { ## must be true, but I'm being defensive + if(code$args[[i]]$args[[1]]$isName) { + insertExprClassLayer(code, i, 'passByMap') + } + } + } + } + } + asserts <- recurseSetSizes(code, symTab, typeEnv) + ## lift any argument that is an expression + for(i in seq_along(code$args)) { + if(inherits(code$args[[i]], 'exprClass')) { + if(!code$args[[i]]$isName) { + forceType <- NULL + if(i==2) forceType <- 'double' ## second argument is always double + asserts <- c(asserts, sizeInsertIntermediate(code, i, symTab, typeEnv, forceType = forceType) ) + } + } + } + a1 <- code$args[[1]] + if(a1$nDim != 2) + stop(exprClassProcessingErrorMsg(code, 'sizePDinverse_logdet called with argument that is not a matrix.'), call. = FALSE) + a1SizeExprs <- a1$sizeExprs + if(is.integer(a1SizeExprs[[1]]) && is.integer(a1SizeExprs[[2]])) { + newSizeExpr <- as.integer(a1SizeExprs[[1]] * a1SizeExprs[[2]] + 1) + } else { + newSizeExpr <- substitute(((A) * (B)) + 1, + list(A = a1SizeExprs[[1]], B = a1SizeExprs[[2]])) + } + code$sizeExprs <- list(newSizeExpr) + code$nDim <- 1 + code$toEigenize <- 'no' + code$type <- 'double' + + ## Self-lift if this expression is amid a larger expression. + if(!(code$caller$name %in% c('{','<-','<<-','='))) { + asserts <- c(asserts, sizeInsertIntermediate(code$caller, code$callerArgID, symTab, typeEnv)) + } else + typeEnv$.ensureNimbleBlocks <- TRUE + + + invisible(asserts) +} + sizeUnaryCwiseSquare <- function(code, symTab, typeEnv) { if(length(code$args) != 1){ stop(exprClassProcessingErrorMsg(code, 'sizeUnaryCwiseSquare called with argument length != 1.'), call. = FALSE) @@ -3779,6 +3851,8 @@ sizeBinaryCwise <- function(code, symTab, typeEnv) { mvFirstArgCheckLists <- list(nimArr_rmnorm_chol = list(c(1, 2, 0), ## dimensionality of ordered arguments AFTER the first, which is for the return value. e.g. mean (1D), chol(2D), prec_param(scalar) 1, 'double'), ## 1 = argument from which to take answer size, double = answer type + nimArr_rmnorm_inv_ld = list(c(1, 2, 1, 0), + 1, 'double'), nimArr_rmvt_chol = list(c(1, 2, 0, 0), ## dimensionality of ordered arguments AFTER the first, which is for the return value. e.g. mean (1D), chol(2D), df(scalar), prec_param(scalar) 1, 'double'), ## 1 = argument from which to take answer size, double = answer type nimArr_rlkj_corr_cholesky = list(c(0, 0), ## eta, p diff --git a/packages/nimble/R/genCpp_toEigenize.R b/packages/nimble/R/genCpp_toEigenize.R index 0994d7c0a..7ec659003 100644 --- a/packages/nimble/R/genCpp_toEigenize.R +++ b/packages/nimble/R/genCpp_toEigenize.R @@ -91,6 +91,7 @@ toEigenCalls <- c( makeCallList(c('nim_IsNA', 'nim_IsNaN'), 'toEigenScalarRecurse'), makeCallList(c('nimArr_dmnorm_chol', + 'nimArr_dmnorm_inv_ld', 'nimArr_dmvt_chol', 'nimArr_dlkj_corr_cholesky', 'nimArr_dwish_chol', @@ -100,6 +101,7 @@ toEigenCalls <- c( 'nimArr_dinterval', 'nimArr_ddirch'), 'toEigenScalarRecurse'), makeCallList(c('nimArr_rmnorm_chol', + 'nimArr_rmnorm_inv_ld', 'nimArr_rmvt_chol', 'nimArr_rlkj_corr_cholesky', 'nimArr_rwish_chol', diff --git a/packages/nimble/R/miscFunctions.R b/packages/nimble/R/miscFunctions.R index 05ee6715f..acbb48d8a 100644 --- a/packages/nimble/R/miscFunctions.R +++ b/packages/nimble/R/miscFunctions.R @@ -63,15 +63,33 @@ calc_dmnormAltParams <- nimbleFunction( } else { I <- diag(dim(cholesky)[1]) ans <- backsolve(cholesky, forwardsolve(t(cholesky), I)) - ## Chris suggests: - ## tmp <- forwardsolve(L, I) - ## ans <- crossprod(tmp) } returnType(double(2)) return(ans) } ) +calc_dmnorm_inv_ld_AltParams <- nimbleFunction( + name = 'calc_dmnorm_inv_ld_AltParams', + run = function(mat = double(2), inv_ld = double(1), prec_param = double(), return_prec = double()) { + ## No need for inversion as desired result is either in `mat` from + ## original parameter provided or inverse is in `inv_ld`. + if(prec_param == return_prec) + return(mat) + + n <- sqrt(length(inv_ld)-1) + nsq <- n * n + ans <- matrix(inv_ld[1:nsq], nrow = n, ncol = n) + if(n > 1) { # Fill lower triangle as PDinverse_logdet only guarantees upper. + for(i in 2:n) + for(j in 1:(i-1)) + ans[i,j] <- ans[j,i] + } + return(ans) + returnType(double(2)) + } +) + ## This is used in conjugacy definition for ddirch, to calculate 'contribution' ## terms from dcat dependents. calc_dcatConjugacyContributions <- nimbleFunction( @@ -113,9 +131,6 @@ calc_dwishAltParams <- nimbleFunction( } else { I <- diag(dim(cholesky)[1]) ans <- backsolve(cholesky, forwardsolve(t(cholesky), I)) - ## Chris suggests: - ## tmp <- forwardsolve(L, I) - ## ans <- crossprod(tmp) } returnType(double(2)) return(ans) diff --git a/packages/nimble/R/nimbleFunction_keywordProcessing.R b/packages/nimble/R/nimbleFunction_keywordProcessing.R index dd4bed90a..d67d42f63 100644 --- a/packages/nimble/R/nimbleFunction_keywordProcessing.R +++ b/packages/nimble/R/nimbleFunction_keywordProcessing.R @@ -818,6 +818,9 @@ nimDerivs_keywordInfo <- keywordInfoClass( ## First check to see if nimFxn argument is a method. fxnCall <- code[[2]][[1]] order <- code[['order']] + outInds <- code[['outInds']] + inDir <- code[['inDir']] + outDir <- code[['outDir']] calculateCase <- FALSE if(deparse(fxnCall) == 'calculate'){ @@ -831,7 +834,15 @@ nimDerivs_keywordInfo <- keywordInfoClass( code[[2]] <- matchKeywordCode(code[[2]], nfProc) calculateCase <- TRUE } + non_NA_arg <- function(arg) { + if(is.call(arg)) return(TRUE) + if(is.name(arg)) return(TRUE) + !isTRUE(all(is.na(arg))) + } if(calculateCase) { + if(non_NA_arg(outInds) || non_NA_arg(inDir) || non_NA_arg(outDir)) { + warning("nimDerivs with calculate() does not support outInds, inDir, or outDir arguments. They will be ignored.") + } innerCode <- code[[2]] if(length(code$wrt)==1) if(!is.name(code$wrt)) @@ -852,6 +863,28 @@ nimDerivs_keywordInfo <- keywordInfoClass( ## Only make the wrt substitution if the names are baked in as character wrtArg <- code$wrt + + # For extending to outInds, inDir, and outDir, we will + # 1. Give warnings if either wrt and inDir are given + # or outInds and outDir are given. + # Treat wrt and outInds as over-riding the others. + # 2. I think we can leave the rest as is, except set NAs to -1s for an easier flag inside C++. + # wrt gets more handling here, resulting in default (-1, -1) if wrt is NULL. + # But sizeProcessing with do makeVectorIfNeeded so the others can be left as is. + + ## We thought about some argument warnings here, but any warnings need to be run-time. + ## wrt and inDir don't make sense in one call, but they can both be run-time so can change betwen calls. + ## outDir and outInds are similar and are even both allowed because in the (run-time) case or order=2 + ## the outInds are used for the Jacobian and the outDir for the Hessian reverse. + for(nameToCheck in c('outInds')) { ## Oops, general loop but only one case. + if(!non_NA_arg(code[[nameToCheck]])) + code[[nameToCheck]] <- substitute(A, list(A=-1)) # get -1 as a number, not parse tree of -(1) + } + for(nameToCheck in c('inDir', 'outDir')) { + if(!non_NA_arg(code[[nameToCheck]])) + code[[nameToCheck]] <- as.numeric(NA) # ensure it is a numeric NA. + } + 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 @@ -1024,7 +1057,8 @@ matchFunctions[['nimDerivs']] <- function(call = NA, order = nimC(0,1,2), wrt = NA, model = NA, updateNodes = NA, constantNodes = NA, - do_update = TRUE, reset = FALSE) {} + do_update = TRUE, reset = FALSE, + outInds=NA, inDir=NA, outDir=NA) {} matchFunctions[['derivInfo']] <- derivInfo matchFunctions[['besselK']] <- function(x, nu, expon.scaled = FALSE){} matchFunctions[['dgamma']] <- function(x, shape, rate = 1, scale, log = FALSE){} diff --git a/packages/nimble/R/options.R b/packages/nimble/R/options.R index e6388c584..ec362b13e 100644 --- a/packages/nimble/R/options.R +++ b/packages/nimble/R/options.R @@ -227,7 +227,8 @@ nimOptimMethod("bobyqa", stripUnusedTypeDefs = TRUE, digits = NULL, enableVirtualNodeFunctionDefs = FALSE, - checkDerivsArgs = TRUE + checkDerivsArgs = TRUE, + useADdmnorm = TRUE ) ) diff --git a/packages/nimble/R/parameterTransform.R b/packages/nimble/R/parameterTransform.R index fe7eb8ee7..697bc38fd 100644 --- a/packages/nimble/R/parameterTransform.R +++ b/packages/nimble/R/parameterTransform.R @@ -181,7 +181,7 @@ parameterTransform <- nimbleFunction( next } stop('`parameterTransform` system doesn\'t have a transformation for the bounds of node: ', node, ', which are (', bounds[1], ', ', bounds[2], ')') } else { ## multivariate - if(dist %in% c('dmnorm', 'dmvt', 'dcar_normal', 'dcar_proper') || ## 6: multivariate {normal, t, CAR}, + if(dist %in% c('dmnorm', 'dmnormAD', 'dmvt', 'dcar_normal', 'dcar_proper') || ## 6: multivariate {normal, t, CAR}, isUserDefined(dist)) ## all multivariate user-defined distributions, { ## and non-scalar determ nodes when allowDeterm is TRUE if(isUserDefined(dist) && getNimbleOption('parameterTransformWarnUserDists')) diff --git a/packages/nimble/inst/CppCode/GNUmakefile b/packages/nimble/inst/CppCode/GNUmakefile index a9ca46f40..e4d889a06 100644 --- a/packages/nimble/inst/CppCode/GNUmakefile +++ b/packages/nimble/inst/CppCode/GNUmakefile @@ -1,4 +1,4 @@ -SRC=nimOptim.cpp nimDerivs_atomic_dyn_ind.cpp nodeFun.cpp accessorClasses.cpp nimbleCppAD.cpp RcppUtils.cpp EigenTypedefs.cpp eigenUsingClasses.cpp RcppNimbleUtils.cpp Utils.cpp NamedObjects.cpp ModelClassUtils.cpp dists.cpp nimDists.cpp dllFinalizer.cpp nimbleGraph.cpp smartPtrs.cpp predefinedNimbleLists.cpp nimIntegrate.cpp nimDerivs_atomic_classes.cpp nimDerivs_atomic_matmult.cpp nimDerivs_atomic_matinverse.cpp nimDerivs_atomic_backsolve.cpp nimDerivs_atomic_forwardsolve.cpp nimDerivs_atomic_cholesky.cpp nimDerivs_atomic_solve_base.cpp nimDerivs_atomic_pow_int.cpp nimDerivs_atomic_zround.cpp nimDerivs_atomic_log_pow_int.cpp nimDerivs_atomic_probit.cpp +SRC=nimDerivs_atomic_PDinverse_logdet.cpp nimOptim.cpp nimDerivs_atomic_dyn_ind.cpp nodeFun.cpp accessorClasses.cpp nimbleCppAD.cpp RcppUtils.cpp EigenTypedefs.cpp eigenUsingClasses.cpp RcppNimbleUtils.cpp Utils.cpp NamedObjects.cpp ModelClassUtils.cpp dists.cpp nimDists.cpp dllFinalizer.cpp nimbleGraph.cpp smartPtrs.cpp predefinedNimbleLists.cpp nimIntegrate.cpp nimDerivs_atomic_classes.cpp nimDerivs_atomic_matmult.cpp nimDerivs_atomic_matinverse.cpp nimDerivs_atomic_backsolve.cpp nimDerivs_atomic_forwardsolve.cpp nimDerivs_atomic_cholesky.cpp nimDerivs_atomic_solve_base.cpp nimDerivs_atomic_pow_int.cpp nimDerivs_atomic_zround.cpp nimDerivs_atomic_log_pow_int.cpp nimDerivs_atomic_probit.cpp OBJECTS=$(SRC:%.cpp=%.o) # nimbleGraph.cpp is only needed in SRC for --enable-dylib=true, but it's not a big deal to leave include it either way diff --git a/packages/nimble/inst/CppCode/dists.cpp b/packages/nimble/inst/CppCode/dists.cpp index 8c8b87f75..f43794e8a 100644 --- a/packages/nimble/inst/CppCode/dists.cpp +++ b/packages/nimble/inst/CppCode/dists.cpp @@ -28,6 +28,7 @@ //#include "Utils.h" // moved to dists.h #include "nimble/dists.h" #include +#include // Detects NA bool R_IsNA_ANY(double* P, int s) { @@ -1211,6 +1212,300 @@ SEXP C_rmnorm_chol(SEXP mean, SEXP chol, SEXP prec_param) return ans; } +// Drafted by GitHub copilot, modified by NIMBLE team. +// This function computes the inverse and log-determinant of a positive-definite matrix. +// Like chol, only the upper triangular part of matPtr is used. +// The first n*n elements of ans will be the inverse, but also only the upper triangular part. +// This means there will be many zeros, which for now we retain because then the vector can be used +// as a matrix within further copying. +// Variables are named as if the matPtr input is a covariance and the output has the precision elements, +// but it can be the other way around (if prec_param == 1 for dmnormAD) +void PDinverse_logdet_internal(double *matPtr, double *ans, int n) { + Eigen::Map> + A(matPtr, n, n); + + double logdet_cov = 0.0; + + //Eigen::MatrixXd chol; + auto chol = (A).template selfadjointView().llt(); + // Eigen::LLT llt(A); // could be used instead of chol below + // chol = EIGEN_CHOL(A); // note: llt could be used directly for solve below. + double half_logdet = 0.0; + for(int i = 0; i < n; ++i) { + half_logdet += std::log(chol.matrixU()(i, i)); + } + logdet_cov += 2.0 * half_logdet; + // Compute precision matrix + Eigen::MatrixXd Prec(n, n); + Prec = chol.solve(Eigen::MatrixXd::Identity(n, n)).template triangularView(); + for(size_t j = 0; j < n; ++j) { + for(size_t i = 0; i < n; ++i) + ans[j * n + i] = Prec(i, j); + } + + // Last element is log(det(cov)) + ans[n * n] = logdet_cov; +} + +// Drafted by GitHub copilot, modified by NIMBLE team. +// This function is the R interface for PDinverse_logdet. +SEXP C_PDinverse_logdet(SEXP mat) { + if(!Rf_isMatrix(mat) || !Rf_isReal(mat)) + RBREAK("Error (C_PDinverse_logdet): 'mat' must be a real matrix.\n"); + + int* dims = INTEGER(Rf_getAttrib(mat, R_DimSymbol)); + if(dims[0] != dims[1]) + RBREAK("Error (C_PDinverse_logdet): 'mat' must be a square matrix.\n"); + int n = dims[0]; + + double* c_mat = REAL(mat); + + SEXP ans; + PROTECT(ans = Rf_allocVector(REALSXP, n * n + 1)); + double* c_ans = REAL(ans); + + PDinverse_logdet_internal(c_mat, c_ans, n); + + UNPROTECT(1); + return ans; +} + +// ...existing code... +// Begin multivariate normal version that uses pre-computed precision with (log(det(cov))) tacked on. +// Drafted by GitHub copoilot, modified by NIMBLE team. +double dmnorm_inv_ld(double* x, double* mean, + double* mat, double* PDinv_ldet, int n, + int prec_param, int give_log, int overwrite_inputs) { + // inv_ld: length n*n+1, first n*n elements are precision matrix (column-major), last is log(det(cov)) + double* xCopy; + double dens = -n * M_LN_SQRT_2PI; + int i, j; + + if (R_IsNA_ANY(x, n) || R_IsNA_ANY(mean, n)) + return NA_REAL; + if (R_IsNaN_ANY(x, n) || R_IsNaN_ANY(mean, n)) + return R_NaN; + if(!prec_param) { + if(R_IsNA_ANY(PDinv_ldet, n*n+1)) + return R_NaN; + if(R_IsNaN_ANY(PDinv_ldet, n*n+1)) + return R_NaN; + if(!R_FINITE_ANY(PDinv_ldet, n*n+1)) return R_D__0; + } else { + if(R_IsNA_ANY(mat, n*n)) + return R_NaN; + if(R_IsNaN_ANY(mat, n*n)) + return R_NaN; + if(R_IsNA(PDinv_ldet[n*n])) + return R_NaN; + if(R_IsNaN(PDinv_ldet[n*n])) + return R_NaN; + if(!R_FINITE_ANY(mat, n*n) || !R_FINITE(PDinv_ldet[n*n])) return R_D__0; + } + if(!R_FINITE_ANY(x, n) || !R_FINITE_ANY(mean, n)) return R_D__0; + + // Add log(det(cov)) term (note: it's log(det(cov)), not log(det(prec))) + double pmhalf = prec_param ? 0.5 : -0.5; + dens += pmhalf * PDinv_ldet[n*n]; + + // Compute (x - mean) + if(overwrite_inputs) { + xCopy = x; + for(i = 0; i < n; i++) + xCopy[i] -= mean[i]; + } else { + xCopy = new double[n]; + for(i = 0; i < n; i++) + xCopy[i] = x[i] - mean[i]; + } + + // First version from copilot did not use BLAS: + // Compute quadratic form: xCopy' * Prec * xCopy + // double quad = 0.0; + // for(i = 0; i < n; i++) { + // double tmp = 0.0; + // for(j = 0; j < n; j++) { + // tmp += inv_ld[j*n + i] * xCopy[j]; // column-major + // } + // quad += xCopy[i] * tmp; + // } + + // I asked for a second version using BLAS: + // Use BLAS for quadratic form: y = Prec * xCopy, quad = xCopy' * y + double* y = new double[n]; + double alpha = 1.0, beta = 0.0; + int inc1 = 1; + double *inv_ld = prec_param ? mat : PDinv_ldet; + //F77_CALL(dgemv)("N", &n, &n, &alpha, inv_ld, &n, xCopy, &inc1, &beta, y, &inc1 FCONE); + F77_CALL(dsymv)("U", &n, &alpha, inv_ld, &n, xCopy, &inc1, &beta, y, &inc1 FCONE); + + double quad = F77_CALL(ddot)(&n, xCopy, &inc1, y, &inc1); + delete [] y; + + dens += -0.5 * quad; + + if(!overwrite_inputs) + delete [] xCopy; + + return give_log ? dens : exp(dens); +} + +// This one needed more human editing for the recycling rule portion. +SEXP C_dmnorm_inv_ld(SEXP x, SEXP mean, SEXP mat, SEXP inv_ld, SEXP prec_param, SEXP return_log) +// calculates mv normal density given precision matrix and log(det(cov)) +// current prec_param is ignored and assume FALSE. +// inv_ld should be a numeric vector of length n*n+1: first n*n elements are precision matrix (column-major), last is log(det(cov)) +{ + if(!Rf_isReal(x) || !Rf_isReal(mean)) + RBREAK("Error (C_dmnorm_inv_ld): 'x' and 'mean' should be real valued.\n"); + if(!Rf_isReal(inv_ld) || !Rf_isLogical(return_log)) + RBREAK("Error (C_dmnorm_inv_ld): invalid input type for one of the arguments.\n"); + if(!Rf_isMatrix(mat) || !Rf_isReal(mat)) + RBREAK("Error (C_dmnorm_inv_ld): 'mat' must be a real matrix.\n"); + if(!Rf_isReal(prec_param)) + RBREAK("Error (C_dmnorm_inv_ld): invalid input type for prec_param.\n"); + int* dims = INTEGER(Rf_getAttrib(mat, R_DimSymbol)); + if(dims[0] != dims[1]) + RBREAK("Error (C_dmnorm_inv_ld): 'mat' must be a square matrix.\n"); + int p = dims[0]; + + int n_x = LENGTH(x); + if(n_x != p) + RBREAK("Error (C_dmnorm_inv_ld): 'x' and 'mat' are not of compatible sizes.\n"); + + if(LENGTH(inv_ld) != n_x*n_x + 1) + RBREAK("Error (C_dmnorm_inv_ld): 'inv_ld' must be length n*n+1.\n"); + + int give_log = (int) LOGICAL(return_log)[0]; + double* c_x = REAL(x); + double* c_mean = REAL(mean); + double * c_mat = REAL(mat); + double* c_inv_ld = REAL(inv_ld); + double c_prec_param = REAL(prec_param)[0]; + + int n_mean = LENGTH(mean); + + // recycling rule for the mean. + double* full_mean; + if(n_mean < n_x) { + full_mean = new double[n_x]; + int i_mean = 0; + for(int i = 0; i < n_x; i++) { + full_mean[i] = c_mean[i_mean++]; + if(i_mean == n_mean) i_mean = 0; + } + } else full_mean = c_mean; + + SEXP ans; + PROTECT(ans = Rf_allocVector(REALSXP, 1)); + REAL(ans)[0] = dmnorm_inv_ld(c_x, full_mean, c_mat, c_inv_ld, n_x, c_prec_param, give_log, 0); + if(n_mean < n_x) + delete [] full_mean; + UNPROTECT(1); + return ans; +} + +// This one needed a fair bit of editing. +void rmnorm_inv_ld(double *ans, double* mean, + double *mat, double* inv_ld, int n, int prec_param) { + // inv_ld: length n*n+1, first n*n elements are precision matrix (column-major), last is log(det(cov)) + // last element of inv_ld is not used. + // Generate standard normal + bool return_NaN = ISNAN_ANY(mean, n); + if(!prec_param) return_NaN |= ISNAN_ANY(inv_ld, n*n) || (!R_FINITE_ANY(inv_ld, n*n)); + else return_NaN |= ISNAN(inv_ld[n*n] || ISNAN_ANY(mat, n*n)) || + (!R_FINITE(inv_ld[n*n]) || !R_FINITE_ANY(mat, n*n)); + if (return_NaN) { + for(int j = 0; j < n; j++) + ans[j] = R_NaN; + return; + } + + for(int i = 0; i < n; i++) + ans[i] = norm_rand(); + + // Cholesky decomposition of precision matrix + double* chol_prec = new double[n*n]; + if(!prec_param) { + for(int i = 0; i < n*n; i++) + chol_prec[i] = inv_ld[i]; + } else { + for(int i = 0; i < n*n; i++) + chol_prec[i] = mat[i]; + } + char uplo = 'U'; + int info(0); + F77_CALL(dpotrf)(&uplo, &n, chol_prec, &n, &info FCONE); + if (info != 0) { + Rf_error("Error in Cholesky decomposition in rnorm_inv_ld: dpotrf returned info != %d", info); + } + // Solve U^T y = z (z = ans), then U x = y + char transPrec = 'N'; + char diag = 'N'; + int incx = 1; + int lda(n); + F77_CALL(dtrsv)(&uplo, &transPrec, &diag, &n, + chol_prec, &lda, ans, &incx FCONE FCONE FCONE); + + // Add mean + for(int i = 0; i < n; i++) + ans[i] += mean[i]; + + delete [] chol_prec; +} + +// also substantial editing +SEXP C_rmnorm_inv_ld(SEXP mean, SEXP mat, SEXP inv_ld, SEXP prec_param) +// generates single mv normal draw given precision matrix and log(det(cov)) +// inv_ld should be a numeric vector of length n*n+1: first n*n elements are precision matrix (column-major), last is log(det(cov)) +{ +if(!Rf_isReal(mean)) + RBREAK("Error (C_rmnorm_inv_ld): 'mean' should be real-valued\n"); + if(!Rf_isReal(inv_ld)) + RBREAK("Error (C_rmnorm_inv_ld): invalid input type for one of the arguments.\n"); + if(!Rf_isMatrix(mat) || !Rf_isReal(mat)) + RBREAK("Error (C_rmnorm_inv_ld): 'mat' must be a real matrix.\n"); + if(!Rf_isReal(prec_param)) + RBREAK("Error (C_rmnorm_inv_ld): invalid input type for prec_param.\n"); + + int* dims = INTEGER(Rf_getAttrib(mat, R_DimSymbol)); + if(dims[0] != dims[1]) + RBREAK("Error (C_rmnorm_inv_ld): 'mat' must be a square matrix.\n"); + int n_mean = LENGTH(mean); + int n_mat = dims[0]; + if(LENGTH(inv_ld) != n_mat*n_mat + 1) + RBREAK("Error (C_rmnorm_inv_ld): 'inv_ld' must be length n*n+1.\n"); + int n_values = n_mat; //sqrt(LENGTH(inv_ld) - 1); + + double* c_mean = REAL(mean); + double c_prec_param = REAL(prec_param)[0]; + double* c_inv_ld = REAL(inv_ld); + double* c_mat = REAL(mat); + double* full_mean; + + // recycling rule for the mean. + if(n_mean < n_values) { + full_mean = new double[n_values]; + int i_mean = 0; + for(int i = 0; i < n_values; i++) { + full_mean[i] = c_mean[i_mean++]; + if(i_mean == n_mean) i_mean = 0; + } + } else full_mean = c_mean; + + GetRNGstate(); + + SEXP ans; + PROTECT(ans = Rf_allocVector(REALSXP, n_mean)); + rmnorm_inv_ld(REAL(ans), full_mean, c_mat, c_inv_ld, n_mean, c_prec_param); + + PutRNGstate(); + if(n_mean < n_values) + delete [] full_mean; + UNPROTECT(1); + return ans; +} +// ...existing code... // Begin multivariate t double dmvt_chol(double* x, double* mu, double* chol, double df, int n, double prec_param, int give_log, int overwrite_inputs) { diff --git a/packages/nimble/inst/CppCode/nimDerivs_atomic_PDinverse_logdet.cpp b/packages/nimble/inst/CppCode/nimDerivs_atomic_PDinverse_logdet.cpp new file mode 100644 index 000000000..859ee25dd --- /dev/null +++ b/packages/nimble/inst/CppCode/nimDerivs_atomic_PDinverse_logdet.cpp @@ -0,0 +1,487 @@ +#include +#include + +/* +Atomic class for matrix inverse and log determinant. +See PDinverse_logdet. +This implementation follows the successful pattern of TMB. +The motivating use case is the multivariate normal distribution. + +See the matinverse atomic notes as well. + +Input X is a positive definite (PD) matrix. +Output Y is a vector of precision elements followed +by a scalar that is the log determinant of the covariance. +These can internally be computed with Cholesky decomposition, +but Cholesky is not placed on the AD tape. Rather, derivatives +for the matrix inverse use the rules for matrix inverse +and derivatives for the log determinant use the rules for +log determinant (Jacobi's formula). The latter uses the +matrix inverse, giving a second reason (beyond shared use +of Cholesky) to do these steps together. + +X could represent the precision or the covariance matrix. +To-Do: Figure out how to handle prec vs. cov. +---- +For covariance matrix X: +Value: +Y = c(X^{-1}, log(det(X))) = c(Y1, Y2) +X is n-x-n. +Y is n*n + 1 vector. Y1 are the n*n precision elements, Y2 is log(det(X)). + +---- +Forward first order +dY1 = -Y1 * dX * Y1 (see atomic_matinverse notes) + +// Extend to the case of symmetric matrices represented by triangular parts. + +dY2 = trace(Y1 * dX) // Jacobi's formula for log determinant + +---- +Reverse first order +From Y1: +Xadjoint = -t(Y1) %*% Yadjoint1 %*% t(Y1); + + +From Y2: + = + = + = sum_i + = sum_i + = sum_i +Xadjoint[, i] += Y1[i,]^T Yadjoint2 +i.e. +Xadjoint += t(Y1) * Yadjoint2 + +---- +Reverse second order +// I am going to wait on this. +*/ + +/* +Storage and usage of diagonal half-matrices: +- +*/ + +atomic_PDinverse_logdet_class::atomic_PDinverse_logdet_class(const std::string& name) : + CppAD::atomic_three(name) {}; + +bool atomic_PDinverse_logdet_class::for_type( + const CppAD::vector& parameter_x , + const CppAD::vector& type_x , + CppAD::vector& type_y ) +{ + // printf("In PDinverse_logdet for_type\n"); + size_t n = type_x.size(); // y is one longer than x + // All types must be the same. + CppAD::ad_type_enum final_type(CppAD::constant_enum); + CppAD::ad_type_enum this_type; + for(size_t i = 0; i < n; ++i) { + this_type = type_x[i]; + if(this_type == CppAD::dynamic_enum) + final_type = CppAD::dynamic_enum; + if(this_type == CppAD::variable_enum) { + final_type = CppAD::variable_enum; + break; + } + } + for(size_t i = 0; i < n+1; ++i) type_y[i] = final_type; + return true; +} + +bool atomic_PDinverse_logdet_class::rev_depend( + const CppAD::vector& parameter_x , + const CppAD::vector& type_x , + CppAD::vector& depend_x , + const CppAD::vector& depend_y + ) { + // printf("In PDinverse_logdet reverse_depend\n"); + bool any_depend_y(false); + size_t ny = depend_y.size(); + for(size_t i = 0; i < ny; ++i) { + if(depend_y[i]) { + any_depend_y = true; + break; + } + } + if(depend_x.size() != ny-1) { + printf("In PDinverse_logdet rev_depend, somehow size of depend_x (n*n) does not match size of depend_y (n*n+1)). That should never happen.\n"); + } + size_t nx = depend_x.size(); + for(size_t i = 0; i < nx; ++i) { + depend_x[i] = any_depend_y; + } + return false; +} + +// Forward mode +inline bool atomic_PDinverse_logdet_class::forward( + const CppAD::vector& parameter_x, + const CppAD::vector& type_x, + size_t need_y, + size_t order_low, + size_t order_up, + const CppAD::vector& taylor_x, + CppAD::vector& taylor_y +) { + //forward mode + // printf("In PDinverse_logdet forward\n"); + size_t nrow = order_up + 1; + size_t n = static_cast< size_t>(sqrt(static_cast(taylor_x.size()/nrow))); + EigenConstMap Xmap(&taylor_x[0], n, n, EigStrDyn(nrow*n, nrow) ); + if((order_low <= 0) && (order_up >= 0)) { // value + EigenMap Ymap(&taylor_y[0], n, n, EigStrDyn(nrow*n, nrow ) ); + + // Assisted by copilot (begin) + // (i) Cholesky decomposition (upper-triangular, as in NIMBLE convention) + auto chol = Xmap.template selfadjointView().llt(); + //Eigen::MatrixXd chol = llt.matrixU(); // Upper-triangular Cholesky factor + + // (ii) Matrix inverse from the Cholesky + // Solve U^T U = X, so X^{-1} = U^{-1} (U^T)^{-1} + Ymap = chol.solve(Eigen::MatrixXd::Identity(n, n)).template triangularView(); + + // (iii) Log determinant from the Cholesky + // log(det(X)) = 2 * sum(log(diag(U))) + double logdet = 0.0; + for(int i = 0; i < n; ++i) + logdet += 2.0 * std::log(chol.matrixU()(i, i)); + // Assisted by copilot (end) + + taylor_y[n * n * nrow + 0] = logdet; // Last element is log(det(X)) + double_cache.set_cache( 0, 0, order_up, taylor_x, taylor_y ); + } + if((order_low <= 1) && (order_up >= 1)) { + // printf("In forward >1\n"); + double_cache.check_and_set_cache(this, + parameter_x, + type_x, + 0, + order_up, + taylor_x, + taylor_y.size()); + size_t cache_nrow = double_cache.nrow(); + EigenMap Ymap(double_cache.taylor_y_ptr(), n, n, EigStrDyn(cache_nrow*n, cache_nrow ) ); + EigenMap dYmap(&taylor_y[1], n, n, EigStrDyn(nrow*n, nrow ) ); + EigenConstMap dXmap(&taylor_x[1], n, n, EigStrDyn(nrow*n, nrow)); + // multiplying two selfadjointViews does not seem to work. + // So, I will materialize at least one of them. If that works, I'll make it a member variable. + MatrixXd Ysym = Ymap.template selfadjointView(); + MatrixXd dXsym = dXmap.template selfadjointView(); +// dYmap = -(Ymap.template selfadjointView() * dXmap.template selfadjointView() * Ymap.template selfadjointView()).template triangularView(); + dYmap = -(Ysym * dXsym * Ysym); + dYmap = dYmap.template triangularView(); + +// std::cout<<"dXmap "< >& parameter_x, + const CppAD::vector& type_x, + size_t need_y, + size_t order_low, + size_t order_up, + const CppAD::vector >& taylor_x, + CppAD::vector >& taylor_y +) { + //forward mode + // printf("In PDinverse_logdet forward\n"); + size_t nrow = order_up + 1; + size_t n = static_cast< size_t>(sqrt(static_cast(taylor_x.size()/nrow))); + metaEigenConstMap Xmap(&taylor_x[0], n, n, EigStrDyn(nrow*n, nrow) ); + if((order_low <= 0) && (order_up >= 0)) { // value + // metaEigenMap Ymap(&taylor_y[0], n, n, EigStrDyn(nrow*n, nrow ) ); + NimArr<2, CppAD::AD > NimArrX; + NimArr<1, CppAD::AD > NimArrY; + NimArrX.setSize(n, n); + for (size_t i = 0; i < n; ++i) + for (size_t j = 0; j < n; ++j) + NimArrX(i, j) = Xmap(i, j); + NimArrY = nimDerivs_PDinverse_logdet(NimArrX); + for (size_t i = 0; i < n*n+1; ++i) + taylor_y[i*nrow] = NimArrY[i]; + CppADdouble_cache.set_cache( 0, 0, order_up, taylor_x, taylor_y ); + } + if((order_low <= 1) && (order_up >= 1)) { + // printf("In meta-forward >1\n"); + CppADdouble_cache.check_and_set_cache(this, + parameter_x, + type_x, + 0, + order_up, + taylor_x, + taylor_y.size()); + size_t cache_nrow = CppADdouble_cache.nrow(); + metaEigenMap Ymap(CppADdouble_cache.taylor_y_ptr(), n, n, EigStrDyn(cache_nrow*n, cache_nrow ) ); + metaEigenMap dYmap(&taylor_y[1], n, n, EigStrDyn(nrow*n, nrow ) ); + metaEigenConstMap dXmap(&taylor_x[1], n, n, EigStrDyn(nrow*n, nrow)); + MatrixXd_CppAD Ysym = Ymap.template selfadjointView(); + MatrixXd_CppAD dXsym = dXmap.template selfadjointView(); +// dYmap = -(Ysym * dXsym * Ysym); + dYmap = nimDerivs_matmult(-Ysym, nimDerivs_matmult( dXsym, Ysym ) ); + dYmap = dYmap.template triangularView(); + // Assisted by copilot (begin) + // Jacobi's formula for log determinant: trace(Y1 * dX) + // Since trace(A*B) = sum_ij A_ij * B_ji, so sum of element-wise product of Ymap and dXmap.transpose() + taylor_y[n * n * nrow + 1] = (Ysym.array() * dXsym.transpose().array()).sum(); + // Assisted by copilot (end) + } + return true; +} + +// Reverse mode +inline bool atomic_PDinverse_logdet_class::reverse( + const CppAD::vector& parameter_x, + const CppAD::vector& type_x, + size_t order_up, + const CppAD::vector& taylor_x, + const CppAD::vector& taylor_y, + CppAD::vector& partial_x, + const CppAD::vector& partial_y +) { + size_t nrow = order_up + 1; + size_t n = static_cast(sqrt(static_cast(taylor_x.size()/nrow))); + + // if(order_up >= 1) { + // EigenConstMap dYmap(&taylor_y[1], n, n, EigStrDyn(nrow*n, nrow ) ); + // std::cout << "Ymap as a matrix starting reverse 2:\n"; + // // Output dYmap as a matrix for readability + // for(int i = 0; i < dYmap.rows(); ++i) { + // for(int j = 0; j < dYmap.cols(); ++j) { + // std::cout << dYmap(i, j) << " "; + // } + // std::cout << std::endl; + // } + // } + + double_cache.check_and_set_cache(this, + parameter_x, + type_x, + 0, // only use cached values up to order 0 + order_up, + taylor_x, + taylor_y.size()); + + size_t cache_nrow = double_cache.nrow(); + EigenConstMap Ymap(double_cache.taylor_y_ptr(), n, n, EigStrDyn(cache_nrow*n, cache_nrow ) ); + MatrixXd Ysym = Ymap.template selfadjointView(); + EigenConstMap Xmap(&taylor_x[0], n, n, EigStrDyn(nrow*n, nrow) ); + MatrixXd Xsym = Xmap.template selfadjointView(); + EigenMap Xadjoint_map(&partial_x[0], n, n, EigStrDyn(nrow*n, nrow) ); + if(order_up >= 0) { + EigenConstMap Yadjoint_map(&partial_y[0], n, n, EigStrDyn(nrow*n, nrow ) ); + Xadjoint_map = (-Ysym.transpose() * Yadjoint_map.template triangularView() * Ysym.transpose()); + for(size_t i = 0; i < n; ++i) { + Xadjoint_map(i, i) += Ysym(i, i) * partial_y[n*n*nrow+0]; // Add the contribution from logdet + for(size_t j = i+1; j < n; ++j) { + Xadjoint_map(i, j) += Xadjoint_map(j, i) + // complete the contribution from inverse + 2 * Ysym(i, j) * partial_y[n*n*nrow+0]; // Add the contribution from logdet + Xadjoint_map(j, i) = 0; + } + } + // Xadjoint_map += partial_y[n*n*nrow+0] * Ysym.transpose(); // Add the contribution from logdet + } + // TO-DO: Make some of the Eigen matrix objects member variables. + if(order_up >= 1) { + EigenConstMap Xdot_map(&taylor_x[1], n, n, EigStrDyn(nrow*n, nrow ) ); // K-dot + MatrixXd Xdot_sym = Xdot_map.template selfadjointView(); //Sigma-dot + + EigenConstMap Ydot_adjoint_map(&partial_y[1], n, n, EigStrDyn(nrow*n, nrow ) ); // J-dot adjoint + MatrixXd Ydot_adjoint_sym = Ydot_adjoint_map.template selfadjointView(); // P-dot adjoint + + EigenMap Xdot_adjoint_map(&partial_x[1], n, n, EigStrDyn(nrow*n, nrow) ); // K-dot adjoint + + //Eigen::MatrixXd Y_Xdot_Y_transpose = (Ysym * Xdot_sym * Ysym).transpose(); + //Eigen::MatrixXd Ydot_adjoint_Ytranspose = Ydot_adjoint_sym * Ysym.transpose(); + Eigen::MatrixXd Y_Jdotadjoint_Y = Ysym * Ydot_adjoint_map.template triangularView() * Ysym; + Eigen::MatrixXd Y_Xdot = Ysym * Xdot_sym; + Eigen::MatrixXd new_Xadjoint_terms = (Y_Jdotadjoint_Y * Y_Xdot.transpose() + + Y_Xdot * Y_Jdotadjoint_Y); + Eigen::MatrixXd new_Xadjoint_terms_from_logdet = -Y_Xdot * Ysym; + for(size_t i = 0; i < n; ++i) { + Xadjoint_map(i, i) += new_Xadjoint_terms(i, i); // from inverse to diag adj of X + Xdot_adjoint_map(i, i) = -Y_Jdotadjoint_Y(i, i); // from inverse to diag adj of Xdot + Xadjoint_map(i, i) += new_Xadjoint_terms_from_logdet(i, i) * partial_y[n*n*nrow+1]; // from logdet to diag adj of X + Xdot_adjoint_map(i, i) += Ysym(i, i) * partial_y[n*n*nrow+1]; // from logdet to diag adj of Xdot + + for(size_t j = i+1; j < n; ++j) { + Xadjoint_map(i, j) += new_Xadjoint_terms(i, j) + new_Xadjoint_terms(j, i); // from inverse to adj of X + Xdot_adjoint_map(i, j) = -Y_Jdotadjoint_Y(i, j)-Y_Jdotadjoint_Y(j, i); // from inverse to adj of Xdot + Xadjoint_map(i, j) += (new_Xadjoint_terms_from_logdet(i,j) + + new_Xadjoint_terms_from_logdet(j,i)) * partial_y[n*n*nrow+1]; // from logdet to adj of X + Xdot_adjoint_map(i, j) += 2* Ysym(i, j) * partial_y[n*n*nrow+1]; // from logdet to adj of Xdot + Xdot_adjoint_map(j, i) = 0; + } + } + // Xdot_adjoint_map = -(Y_Jdotadjoint_Y.template triangularView()); + } +// Xadjoint_map = Xadjoint_map.template triangularView(); + return true; +} + +inline bool atomic_PDinverse_logdet_class::reverse( + const CppAD::vector >& parameter_x , + const CppAD::vector& type_x , + size_t order_up , + const CppAD::vector >& taylor_x , + const CppAD::vector >& taylor_y , + CppAD::vector >& partial_x , + const CppAD::vector >& partial_y ) +{ + size_t nrow = order_up + 1; + size_t n = static_cast(sqrt(static_cast(taylor_x.size()/nrow))); + + CppADdouble_cache.check_and_set_cache(this, + parameter_x, + type_x, + 0, // only use cached values up to order 0 + order_up, + taylor_x, + taylor_y.size()); + size_t cache_nrow = CppADdouble_cache.nrow(); + metaEigenConstMap Ymap(CppADdouble_cache.taylor_y_ptr(), n, n, EigStrDyn(cache_nrow*n, cache_nrow ) ); + MatrixXd_CppAD Ysym = Ymap.template selfadjointView(); + metaEigenConstMap Xmap(&taylor_x[0], n, n, EigStrDyn(nrow*n, nrow) ); + MatrixXd_CppAD Xsym = Xmap.template selfadjointView(); + metaEigenMap Xadjoint_map(&partial_x[0], n, n, EigStrDyn(nrow*n, nrow) ); + // QUESTION: Save recorded operations by hand-coding due to the triangular stuff. + if(order_up >= 0) { + metaEigenConstMap Yadjoint_map(&partial_y[0], n, n, EigStrDyn(nrow*n, nrow ) ); + //Xadjoint_map = (-Ysym.transpose() * Yadjoint_map.template triangularView() * Ysym.transpose()); + Xadjoint_map = nimDerivs_matmult(-Ysym.transpose(), + nimDerivs_matmult(Yadjoint_map.template triangularView(), + Ysym.transpose())); + for(size_t i = 0; i < n; ++i) { + Xadjoint_map(i, i) += Ysym(i, i) * partial_y[n*n*nrow+0]; // Add the contribution from logdet + for(size_t j = i+1; j < n; ++j) { + Xadjoint_map(i, j) += Xadjoint_map(j, i) + // complete the contribution from inverse + 2 * Ysym(i, j) * partial_y[n*n*nrow+0]; // Add the contribution from logdet + Xadjoint_map(j, i) = 0; + } + } + } + + if(order_up >= 1) { + metaEigenConstMap Xdot_map(&taylor_x[1], n, n, EigStrDyn(nrow*n, nrow ) ); + MatrixXd_CppAD Xdot_sym = Xdot_map.template selfadjointView(); //Sigma-dot + metaEigenConstMap Ydot_adjoint_map(&partial_y[1], n, n, EigStrDyn(nrow*n, nrow ) ); + MatrixXd_CppAD Ydot_adjoint_sym = Ydot_adjoint_map.template selfadjointView(); // P-dot adjoint + metaEigenMap Xdot_adjoint_map(&partial_x[1], n, n, EigStrDyn(nrow*n, nrow) ); + //metaEigenMatrixXd Y_Xdot_Y_transpose = nimDerivs_matmult(Ymap, nimDerivs_matmult( Xdot_map, Ymap)).transpose(); + //metaEigenMatrixXd Ydot_adjoint_Ytranspose = nimDerivs_matmult(Ydot_adjoint_map, Ymap.transpose()); + MatrixXd_CppAD Y_Jdotadjoint_Y = nimDerivs_matmult(Ysym, + nimDerivs_matmult(Ydot_adjoint_map.template triangularView(), Ysym)); + MatrixXd_CppAD Y_Xdot = nimDerivs_matmult(Ysym, Xdot_sym); + MatrixXd_CppAD new_Xadjoint_terms = + nimDerivs_matmult(Y_Jdotadjoint_Y, Y_Xdot.transpose()) + + nimDerivs_matmult(Y_Xdot, Y_Jdotadjoint_Y); + MatrixXd_CppAD new_Xadjoint_terms_from_logdet = nimDerivs_matmult(-Y_Xdot, Ysym); + + for(size_t i = 0; i < n; ++i) { + Xadjoint_map(i, i) += new_Xadjoint_terms(i, i); + Xdot_adjoint_map(i, i) = -Y_Jdotadjoint_Y(i, i); + + Xadjoint_map(i, i) += new_Xadjoint_terms_from_logdet(i, i) * partial_y[n*n*nrow+1]; // from logdet to diag adj of X + + Xdot_adjoint_map(i, i) += Ysym(i, i) * partial_y[n*n*nrow+1]; + for(size_t j = i+1; j < n; ++j) { + Xadjoint_map(i, j) += new_Xadjoint_terms(i, j) + new_Xadjoint_terms(j, i); + Xdot_adjoint_map(i, j) = -Y_Jdotadjoint_Y(i, j)- Y_Jdotadjoint_Y(j, i); + Xadjoint_map(i, j) += (new_Xadjoint_terms_from_logdet(i,j) + + new_Xadjoint_terms_from_logdet(j,i)) * partial_y[n*n*nrow+1]; // from logdet to adj of X + Xdot_adjoint_map(i, j) += 2* Ysym(i, j) * partial_y[n*n*nrow+1]; + Xdot_adjoint_map(j, i) = 0; + } + } + } + return true; +} + +/* void PDinverse_logdet(const MatrixXd_CppAD &x, // This (non-template) type forces any incoming expression to be evaluated + MatrixXd_CppAD &y) { + // static PDinverse_logdet_class PDinverse_logdet("PDinverse_logdet"); // this has no state information so the same object can be used for all cases + atomic_PDinverse_logdet_class *PDinverse_logdet; + size_t n = x.rows(); + std::vector > xVec(n*n); + mat2vec(x, xVec); + std::vector > yVec(n*n+1); + bool recording = CppAD::AD::get_tape_handle_nimble() != nullptr; + if(!recording) { + PDinverse_logdet = new atomic_PDinverse_logdet_class("PDinverse_logdet"); + } else { + void *tape_mgr = CppAD::AD::get_tape_handle_nimble()->nimble_CppAD_tape_mgr_ptr(); + PDinverse_logdet = new_atomic_PDinverse_logdet(tape_mgr, "PDinverse_logdet"); + } + (*PDinverse_logdet)(xVec, yVec); + y.resize(n*n+1, 1); + for(size_t i = 0; i < n*n+1; ++i) { + y(i) = yVec[i]; + } + if(!recording) { + delete PDinverse_logdet; + } else { + track_nimble_atomic(PDinverse_logdet, + CppAD::AD::get_tape_handle_nimble()->nimble_CppAD_tape_mgr_ptr(), + CppAD::local::atomic_index_info_vec_manager_nimble::manage() ); + } +} */ + +/* +MatrixXd_CppAD nimDerivs_PDinverse_logdet(const MatrixXd_CppAD &x, const CppAD::AD &prec_param) { + // prec_param is ignored for now. + MatrixXd_CppAD ans; + PDinverse_logdet(x, ans); + return ans; +} +*/ + +NimArr<1, CppAD::AD > nimDerivs_PDinverse_logdet(const NimArr<2, CppAD::AD > &x) { + // prec_param was being ignored and has been removed for now. + atomic_PDinverse_logdet_class *PDinverse_logdet; + size_t n = x.dim()[0]; + std::vector > xVec(n*n); + for(size_t i = 0; i < n; ++i) { + for(size_t j = 0; j < n; ++j) { + xVec[j*n + i] = x(i, j); + } + } + std::vector > yVec(n*n+1); + bool recording = CppAD::AD::get_tape_handle_nimble() != nullptr; + if(!recording) { + PDinverse_logdet = new atomic_PDinverse_logdet_class("PDinverse_logdet"); + } else { + void *tape_mgr = CppAD::AD::get_tape_handle_nimble()->nimble_CppAD_tape_mgr_ptr(); + PDinverse_logdet = new_atomic_PDinverse_logdet(tape_mgr, "PDinverse_logdet"); + } + (*PDinverse_logdet)(xVec, yVec); + NimArr<1, CppAD::AD > y; + y.initialize(0, false, n*n+1); + for(size_t i = 0; i < n*n+1; ++i) { + y(i) = yVec[i]; + } + if(!recording) { + delete PDinverse_logdet; + } else { + track_nimble_atomic(PDinverse_logdet, + CppAD::AD::get_tape_handle_nimble()->nimble_CppAD_tape_mgr_ptr(), + CppAD::local::atomic_index_info_vec_manager_nimble::manage() ); + } + return y; +} diff --git a/packages/nimble/inst/CppCode/nimDerivs_atomic_classes.cpp b/packages/nimble/inst/CppCode/nimDerivs_atomic_classes.cpp index 5c3ea8a33..09a43eb5f 100644 --- a/packages/nimble/inst/CppCode/nimDerivs_atomic_classes.cpp +++ b/packages/nimble/inst/CppCode/nimDerivs_atomic_classes.cpp @@ -26,6 +26,7 @@ ATOMIC_NEW_DELETE_(backsolve) ATOMIC_NEW_DELETE_(cholesky) ATOMIC_NEW_DELETE_(forwardsolve) ATOMIC_NEW_DELETE_(matinverse) +ATOMIC_NEW_DELETE_(PDinverse_logdet) ATOMIC_NEW_DELETE_(matmult) ATOMIC_NEW_DELETE_(pow_int) ATOMIC_NEW_DELETE_(zround) @@ -33,6 +34,7 @@ ATOMIC_NEW_DELETE_(floor) ATOMIC_NEW_DELETE_(ceil) ATOMIC_NEW_DELETE_(ftrunc) ATOMIC_NEW_DELETE_(nimRound) +ATOMIC_NEW_DELETE_(nimStep) ATOMIC_NEW_DELETE_(log_pow_int) ATOMIC_NEW_DELETE_(zb_over_a) ATOMIC_NEW_DELETE_(probit) @@ -149,6 +151,19 @@ atomic_nimRound_class *nimble_CppAD_tape_mgr::get_atomic_nimRound(std::vector(atomic_ptrs[nimRound_index].first); } +atomic_nimStep_class *track_atomic_nimStep(void* tape_mgr_ptr, + std::vector* vec_ptr) { + return reinterpret_cast(tape_mgr_ptr)->get_atomic_nimStep(vec_ptr); +} +atomic_nimStep_class *nimble_CppAD_tape_mgr::get_atomic_nimStep(std::vector* vec_ptr) { + if(!nimStep_exists) { + nimStep_index = atomic_ptrs.size(); + atomic_ptrs.push_back(atomic_pair(new_atomic_nimStep("atomic_nimStep_managed"), vec_ptr) ); + nimStep_exists = true; + } + return dynamic_cast(atomic_ptrs[nimStep_index].first); +} + atomic_probit_class *track_atomic_probit(void* tape_mgr_ptr, std::vector* vec_ptr) { return reinterpret_cast(tape_mgr_ptr)->get_atomic_probit(vec_ptr); diff --git a/packages/nimble/inst/CppCode/nimDerivs_atomic_matinverse.cpp b/packages/nimble/inst/CppCode/nimDerivs_atomic_matinverse.cpp index 245141573..84372013e 100644 --- a/packages/nimble/inst/CppCode/nimDerivs_atomic_matinverse.cpp +++ b/packages/nimble/inst/CppCode/nimDerivs_atomic_matinverse.cpp @@ -45,7 +45,7 @@ dFdot = -(dY * Xdot * Y) -(Y * dXdot * Y) - (Y * Xdot * dY), dFdot = (Y * dX * Y * Xdot * Y) - (Y * dXdot * Y) + (Y * dXdot * Y * dX * Y) = + + = + <-Y^T * Ydot_adjoint * Y^T , dXdot > + <(Y * Xdot * Y)^T * Ydot_adjoint * Y^T , dX> - = + <-Y^T * Ydot_adjoint * Y^T , dXdot > + = + <-Y^T * Ydot_adjoint * Y^T , dXdot > = + */ diff --git a/packages/nimble/inst/CppCode/nimDerivs_atomic_zround.cpp b/packages/nimble/inst/CppCode/nimDerivs_atomic_zround.cpp index e830c1af4..e80047c28 100644 --- a/packages/nimble/inst/CppCode/nimDerivs_atomic_zround.cpp +++ b/packages/nimble/inst/CppCode/nimDerivs_atomic_zround.cpp @@ -250,3 +250,33 @@ CppAD::AD nimDerivs_nimRound(const CppAD::AD x) { return out[0]; } +/**************************/ + +double nimDerivs_nimStep_class::operator()(double x) {return nimStep(x);} +double nimDerivs_nimStep_class::operator()(int x) {return nimStep(x);} +CppAD::AD nimDerivs_nimStep_class::operator()(const CppAD::AD &x) {return nimDerivs_nimStep(x);} + +atomic_nimStep_class::atomic_nimStep_class(const std::string& name) : + atomic_discrete_class(name) +{ } + +CppAD::AD nimDerivs_nimStep(const CppAD::AD x) { + atomic_nimStep_class* atomic_nimStep; + bool recording = CppAD::AD::get_tape_handle_nimble() != nullptr; + if(!recording) { + atomic_nimStep = new atomic_nimStep_class("atomic_nimStep"); + } else { + atomic_nimStep = track_atomic_nimStep(CppAD::AD::get_tape_handle_nimble()->nimble_CppAD_tape_mgr_ptr(), + CppAD::local::atomic_index_info_vec_manager_nimble::manage() ); + } + CppAD::vector< CppAD::AD > in(1); + in[0] = x; + CppAD::vector< CppAD::AD > out(1); + (*atomic_nimStep)(in, out); + if(!recording) { + delete atomic_nimStep; + } + return out[0]; +} + + diff --git a/packages/nimble/inst/CppCode/nimDists.cpp b/packages/nimble/inst/CppCode/nimDists.cpp index e976e73bf..46320ec94 100644 --- a/packages/nimble/inst/CppCode/nimDists.cpp +++ b/packages/nimble/inst/CppCode/nimDists.cpp @@ -315,6 +315,92 @@ void nimArr_rmnorm_chol(NimArr<1, double> &ans, NimArr<1, double> &mean, NimArr< // } } +// Drafted by GitHub copilot, modified by NIMBLE team. +// Compute a vector of matrix inverse elements and log determinant of covariance matrix. +// The log determinant is computed based on the Cholesky decomposition. +// The output is in the format needed for dmnorm_inv_ld function. +NimArr<1, double> PDinverse_logdet(NimArr<2, double> &mat) { + int n = mat.dimSize(0); + NimArr<1, double> out; + out.setSize(n * n + 1); + + NimArr<2, double> matCopy; + double* matPtr = nimArrCopyIfNeeded<2, double>(mat, matCopy).getPtr(); + + // Use the internal function for the actual computation + PDinverse_logdet_internal(matPtr, out.getPtr(), n); + + return out; +} + +// drafted by GitHub copilot (begin) +double nimArr_dmnorm_inv_ld(NimArr<1, double> &x, NimArr<1, double> &mean, + NimArr<2, double> &mat, + NimArr<1, double> &inv_ld, int prec_param, + int give_log, int overwrite_inputs) { + double *xptr, *meanptr, *matptr, *inv_ld_ptr; + NimArr<1, double> xCopy, meanCopy, inv_ldCopy; + NimArr<2, double> matCopy; + xptr = nimArrCopyIfNeeded<1, double>(x, xCopy).getPtr(); + int n = x.size(); + meanptr = nimArrCopyIfNeeded<1, double>(mean, meanCopy).getPtr(); + if(mean.size() != n) { + _nimble_global_output<<"Error in nimArr_dmnorm_inv_ld: mean and x are different sizes.\n"; + nimble_print_to_R(_nimble_global_output); + } + inv_ld_ptr = nimArrCopyIfNeeded<1, double>(inv_ld, inv_ldCopy).getPtr(); + if(inv_ld.size() != n*n + 1) { + _nimble_global_output<<"Error in nimArr_dmnorm_inv_ld: inv_ld size ("<(mat, matCopy).getPtr(); + if((mat.dim()[0] != n) || (mat.dim()[1] != n)) { + _nimble_global_output<<"Error in nimArr_dmnorm_inv_ld: mat does not match size of x.\n"; + nimble_print_to_R(_nimble_global_output); + } + + double ans; + ans = dmnorm_inv_ld(xptr, meanptr, matptr, inv_ld_ptr, n, prec_param, give_log, overwrite_inputs); + return(ans); +} + +void nimArr_rmnorm_inv_ld(NimArr<1, double> &ans, NimArr<1, double> &mean, + NimArr<2, double> &mat, NimArr<1, double> &inv_ld, int prec_param) { + NimArr<1, double> ansCopy, meanCopy, inv_ldCopy; + NimArr<2, double> matCopy; + double *ansPtr, *meanPtr, *matPtr, *inv_ldPtr; + + int n = mean.size(); + if(!ans.isMap()) { + ans.setSize(n); + } else { + if(ans.size() != n) { + _nimble_global_output<<"Error in nimArr_rmnorm_inv_ld: answer size ("<< ans.size() <<") does not match mean size ("<(ans, ansCopy).getPtr(); + meanPtr = nimArrCopyIfNeeded<1, double>(mean, meanCopy).getPtr(); + matPtr = nimArrCopyIfNeeded<2, double>(mat, matCopy).getPtr(); + inv_ldPtr = nimArrCopyIfNeeded<1, double>(inv_ld, inv_ldCopy).getPtr(); + + rmnorm_inv_ld(ansPtr, meanPtr, matPtr, inv_ldPtr, n, prec_param); + + if(ansPtr != ans.getPtr()) {ans = ansCopy;} + // if(ans.isMap()) { + // ans = ansCopy; + // } +} +// drafted by GitHub copilot (end) + // Begin multivariate t double nimArr_dmvt_chol(NimArr<1, double> &x, NimArr<1, double> &mu, NimArr<2, double> &chol, double df, double prec_param, int give_log, int overwrite_inputs) { diff --git a/packages/nimble/inst/CppCode/nimbleCppAD.cpp b/packages/nimble/inst/CppCode/nimbleCppAD.cpp index 918de9869..c79fa211a 100644 --- a/packages/nimble/inst/CppCode/nimbleCppAD.cpp +++ b/packages/nimble/inst/CppCode/nimbleCppAD.cpp @@ -838,6 +838,13 @@ NimArr<1, double> make_vector_if_necessary(NimArr<1, int> a){ return(intArray); } +NimArr<1, CppAD::AD > make_vector_if_necessary(CppAD::AD x) { + NimArr<1, CppAD::AD > NimArr_x; + NimArr_x.setSize(1, 0, 0); + NimArr_x[0] = x; + return NimArr_x; +} + void setValues_AD_AD_taping(NimArr<1, CppAD::AD > &v, ManyVariablesMapAccessor &MVA_AD, ManyVariablesMapAccessor &MVA_orig, @@ -1052,6 +1059,8 @@ void nimble_CppAD_tape_mgr::reset() { ftrunc_index = 0; nimRound_exists = false; nimRound_index = 0; + nimStep_exists = false; + nimStep_index = 0; log_pow_int_exists = false; log_pow_int_index = 0; zb_over_a_exists = false; @@ -1110,6 +1119,8 @@ nimble_CppAD_tape_mgr::nimble_CppAD_tape_mgr() : ftrunc_exists(false), nimRound_index(0), nimRound_exists(false), + nimStep_index(0), + nimStep_exists(false), log_pow_int_index(0), log_pow_int_exists(false), zb_over_a_index(0), diff --git a/packages/nimble/inst/include/cppad/core/ad.hpp b/packages/nimble/inst/include/cppad/core/ad.hpp index 6c68d340e..e54146777 100644 --- a/packages/nimble/inst/include/cppad/core/ad.hpp +++ b/packages/nimble/inst/include/cppad/core/ad.hpp @@ -288,7 +288,17 @@ private : taddr_ = taddr; ad_type_ = variable_enum; } - // --------------------------------------------------------------- + // ----------------------------------------------------------------- + // make_dynamic added by Perry de Valpine + // Make this parameter a new dynamic + void make_dynamic(tape_id_t id, addr_t taddr) + { CPPAD_ASSERT_UNKNOWN( Parameter(*this) ); // currently a par + CPPAD_ASSERT_UNKNOWN( taddr > 0 ); // sure valid taddr + + tape_id_ = id; + taddr_ = taddr; + ad_type_ = dynamic_enum; + }// --------------------------------------------------------------- // tape linking functions // // not static diff --git a/packages/nimble/inst/include/cppad/core/add.hpp b/packages/nimble/inst/include/cppad/core/add.hpp index cb28ceca2..21411d270 100644 --- a/packages/nimble/inst/include/cppad/core/add.hpp +++ b/packages/nimble/inst/include/cppad/core/add.hpp @@ -106,7 +106,17 @@ AD operator + (const AD &left , const AD &right) } } else if( dyn_left | dyn_right ) - { addr_t arg0 = left.taddr_; + { if( (! dyn_left) && IdenticalZero(left.value_) ) + { + result.make_dynamic(right.tape_id_, right.taddr_); + } + else if( (! dyn_right) && IdenticalZero(right.value_) ) + { + result.make_dynamic(left.tape_id_, left.taddr_); + } + else + { + addr_t arg0 = left.taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left.value_); @@ -119,6 +129,7 @@ AD operator + (const AD &left , const AD &right) ); result.tape_id_ = tape_id; result.ad_type_ = dynamic_enum; + } } return result; } diff --git a/packages/nimble/inst/include/cppad/core/add_eq.hpp b/packages/nimble/inst/include/cppad/core/add_eq.hpp index 02c516667..923231dab 100644 --- a/packages/nimble/inst/include/cppad/core/add_eq.hpp +++ b/packages/nimble/inst/include/cppad/core/add_eq.hpp @@ -104,7 +104,16 @@ AD& AD::operator += (const AD &right) } } else if( dyn_left | dyn_right ) - { addr_t arg0 = taddr_; + { if( (! dyn_right) && IdenticalZero(right.value_) ) + { // this is left += 0, so do nothing + } + else if( (! dyn_left) && IdenticalZero(left)) + { // this is 0 += right + make_dynamic(right.tape_id_, right.taddr_); + } + else + { + addr_t arg0 = taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left); @@ -117,6 +126,7 @@ AD& AD::operator += (const AD &right) ); tape_id_ = tape_id; ad_type_ = dynamic_enum; + } } return *this; } diff --git a/packages/nimble/inst/include/cppad/core/azmul.hpp b/packages/nimble/inst/include/cppad/core/azmul.hpp index d091b24e2..a9ac1379f 100644 --- a/packages/nimble/inst/include/cppad/core/azmul.hpp +++ b/packages/nimble/inst/include/cppad/core/azmul.hpp @@ -190,7 +190,23 @@ azmul(const AD& x, const AD& y) } } else if( dyn_x | dyn_y ) - { addr_t arg0 = x.taddr_; + { if( (! dyn_x) && IdenticalZero(x.value_) ) + { // result = 0 * dynamic + } + else if ( ( ! dyn_y ) && IdenticalZero(y.value_) ) + { // result = dynamic * 0 + } + else if( ( ! dyn_x ) && IdenticalOne(x.value_) ) + { // result = 1 * dynamic + result.make_dynamic(y.tape_id_, y.taddr_); + } + else if( ( ! dyn_y ) && IdenticalOne(y.value_) ) + { // result = dynamic * 1 + result.make_dynamic(x.tape_id_, x.taddr_ ); + } + else + { + addr_t arg0 = x.taddr_; addr_t arg1 = y.taddr_; if( ! dyn_x ) arg0 = tape->Rec_.put_con_par(x.value_); @@ -203,6 +219,7 @@ azmul(const AD& x, const AD& y) ); result.tape_id_ = tape_id; result.ad_type_ = dynamic_enum; + } } return result; } diff --git a/packages/nimble/inst/include/cppad/core/div.hpp b/packages/nimble/inst/include/cppad/core/div.hpp index db714ee46..4b114eb2e 100644 --- a/packages/nimble/inst/include/cppad/core/div.hpp +++ b/packages/nimble/inst/include/cppad/core/div.hpp @@ -104,7 +104,15 @@ AD operator / (const AD &left , const AD &right) } } else if( dyn_left | dyn_right ) - { addr_t arg0 = left.taddr_; + { if ( (!dyn_left) && IdenticalZero(left.value_)) + { // 0 / dynamic + result.value_ = Base(0.0); + } else if( (!dyn_right) && IdenticalOne(right.value_)) + { // dynamic / 1 + result.make_dynamic(left.tape_id_, left.taddr_); + } else + { + addr_t arg0 = left.taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left.value_); @@ -117,6 +125,7 @@ AD operator / (const AD &left , const AD &right) ); result.tape_id_ = tape_id; result.ad_type_ = dynamic_enum; + } } return result; } diff --git a/packages/nimble/inst/include/cppad/core/div_eq.hpp b/packages/nimble/inst/include/cppad/core/div_eq.hpp index 365b90e9d..5a4f985f6 100644 --- a/packages/nimble/inst/include/cppad/core/div_eq.hpp +++ b/packages/nimble/inst/include/cppad/core/div_eq.hpp @@ -105,7 +105,12 @@ AD& AD::operator /= (const AD &right) } } else if( dyn_left | dyn_right ) - { addr_t arg0 = taddr_; + { if( (!dyn_right) && IdenticalOne(right.value_)) + { // dynamic /= 1, so do nothing + } else if( (!dyn_left) && IdenticalZero(left)) + { // 0 /= dynamic, so do nothing + } else + { addr_t arg0 = taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left); @@ -118,6 +123,7 @@ AD& AD::operator /= (const AD &right) ); tape_id_ = tape_id; ad_type_ = dynamic_enum; + } } return *this; } diff --git a/packages/nimble/inst/include/cppad/core/mul.hpp b/packages/nimble/inst/include/cppad/core/mul.hpp index 14365be31..4712cc5f1 100644 --- a/packages/nimble/inst/include/cppad/core/mul.hpp +++ b/packages/nimble/inst/include/cppad/core/mul.hpp @@ -111,7 +111,21 @@ AD operator * (const AD &left , const AD &right) } } else if( dyn_left | dyn_right ) - { addr_t arg0 = left.taddr_; + { if ( (!dyn_left) && IdenticalZero(left.value_)) + { // 0 * dynamic + result.value_ = Base(0.0); + } else if( (!dyn_right) && IdenticalZero(right.value_)) + { // dynamic * 0 + result.value_ = Base(0.0); + } else if( (!dyn_left) && IdenticalOne(left.value_)) + { // 1 * dynamic + result.make_dynamic(right.tape_id_, right.taddr_); + } else if( (!dyn_right) && IdenticalOne(right.value_)) + { // dynamic * 1 + result.make_dynamic(left.tape_id_, left.taddr_); + } else + { + addr_t arg0 = left.taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left.value_); @@ -124,6 +138,7 @@ AD operator * (const AD &left , const AD &right) ); result.tape_id_ = tape_id; result.ad_type_ = dynamic_enum; + } } return result; } diff --git a/packages/nimble/inst/include/cppad/core/mul_eq.hpp b/packages/nimble/inst/include/cppad/core/mul_eq.hpp index 6b75e97f3..cd3030b05 100644 --- a/packages/nimble/inst/include/cppad/core/mul_eq.hpp +++ b/packages/nimble/inst/include/cppad/core/mul_eq.hpp @@ -114,7 +114,19 @@ AD& AD::operator *= (const AD &right) } } else if( dyn_left | dyn_right ) - { addr_t arg0 = taddr_; + { if( (!dyn_right) && IdenticalOne(right.value_)) + { // dynamic *= 1, so do nothing + } else if( (!dyn_right) && IdenticalZero(right.value_)) + { // dynamic *= 0, so remove from tape + tape_id_ = 0; + } else if( (!dyn_left) && IdenticalOne(left)) + { // 1 *= dynamic + make_dynamic(right.tape_id_, right.taddr_); + } else if( (!dyn_left) && IdenticalZero(left)) + { // 0 *= dynamic, so do nothing + } else + { + addr_t arg0 = taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left); @@ -127,6 +139,7 @@ AD& AD::operator *= (const AD &right) ); tape_id_ = tape_id; ad_type_ = dynamic_enum; + } } return *this; } diff --git a/packages/nimble/inst/include/cppad/core/sub.hpp b/packages/nimble/inst/include/cppad/core/sub.hpp index 1aaeeac6f..c566d21f5 100644 --- a/packages/nimble/inst/include/cppad/core/sub.hpp +++ b/packages/nimble/inst/include/cppad/core/sub.hpp @@ -100,7 +100,13 @@ AD operator - (const AD &left , const AD &right) result.ad_type_ = variable_enum; } else if( dyn_left | dyn_right ) - { addr_t arg0 = left.taddr_; + { if( (! dyn_right) && IdenticalZero(right.value_) ) + { // this is dynamic - 0 + result.make_dynamic(left.tape_id_, left.taddr_); + } + else + { + addr_t arg0 = left.taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left.value_); @@ -113,6 +119,7 @@ AD operator - (const AD &left , const AD &right) ); result.tape_id_ = tape_id; result.ad_type_ = dynamic_enum; + } } return result; } diff --git a/packages/nimble/inst/include/cppad/core/sub_eq.hpp b/packages/nimble/inst/include/cppad/core/sub_eq.hpp index 2859b9474..16b0fc128 100644 --- a/packages/nimble/inst/include/cppad/core/sub_eq.hpp +++ b/packages/nimble/inst/include/cppad/core/sub_eq.hpp @@ -100,7 +100,12 @@ AD& AD::operator -= (const AD &right) ad_type_ = variable_enum; } else if( dyn_left | dyn_right ) - { addr_t arg0 = taddr_; + { if( (! dyn_right) && IdenticalZero(right.value_) ) + { // this is left -= 0, so do nothing + } + else + { + addr_t arg0 = taddr_; addr_t arg1 = right.taddr_; if( ! dyn_left ) arg0 = tape->Rec_.put_con_par(left); @@ -113,6 +118,7 @@ AD& AD::operator -= (const AD &right) ); tape_id_ = tape_id; ad_type_ = dynamic_enum; + } } return *this; } diff --git a/packages/nimble/inst/include/nimble/Utils.h b/packages/nimble/inst/include/nimble/Utils.h index 2152be672..b908ef667 100644 --- a/packages/nimble/inst/include/nimble/Utils.h +++ b/packages/nimble/inst/include/nimble/Utils.h @@ -213,9 +213,31 @@ inline T nimDerivs_cloglog(T x){ inline int nimEquals(double x1, double x2) {return(x1 == x2 ? 1 : 0);} //template -inline T nimDerivs_nimEquals(T x1, T x2){ - return(CondExpEq(x1, x2, T(1), T(0))); -} +// inline T nimDerivs_nimEquals(T x1, T x2){ +// T cond = nimDerivs_nimStep(x1-x2) * nimDerivs_nimStep(x2-x1); +// return cond; +// // return(CondExpEq(x1, x2, T(1), T(0))); +// } + +// inline T nimDerivs_CondExpGe(T x1, T x2, T y1, T y2) { +// T cond = nimDerivs_nimStep(x1-x2); +// return cond*y1 + (T(1) - cond)*y2; +// } + +// inline T nimDerivs_CondExpGt(T x1, T x2, T y1, T y2) { +// T cond = 1-nimDerivs_nimStep(x2-x1); +// return cond*y1 + (T(1) - cond)*y2; +// } + +// inline T nimDerivs_CondExpLe(T x1, T x2, T y1, T y2) { +// T cond = nimDerivs_nimStep(x2-x1); +// return cond*y1 + (T(1) - cond)*y2; +// } + +// inline T nimDerivs_CondExpLt(T x1, T x2, T y1, T y2) { +// T cond = 1-nimDerivs_nimStep(x1-x2); +// return cond*y1 + (T(1) - cond)*y2; +// } inline double nimbleIfElse(bool condition, double x1, double x2) {return(condition ? x1 : x2);} inline double lfactorial(double x) {return(lgammafn(x + 1));} @@ -230,22 +252,26 @@ inline T nimDerivs_logit(T x) { inline double nimRound(double x) {return(round(x));} inline double pairmax(double x1, double x2) {return(x1 > x2 ? x1 : x2);} //template -inline T nimDerivs_pairmax(T x1, T x2) { - return(CondExpGt(x1, x2, x1, x2)); -} +// inline T nimDerivs_pairmax(T x1, T x2) { +// T cond = nimStep(x1-x2); +// return cond*x1 + (T(1) - cond)*x2; +// // return(CondExpGt(x1, x2, x1, x2)); +// } inline double pairmin(double x1, double x2) {return(x1 < x2 ? x1 : x2);} //template -inline T nimDerivs_pairmin(T x1, T x2) { - return(CondExpLt(x1, x2, x1, x2)); -} +// inline T nimDerivs_pairmin(T x1, T x2) { +// T cond = nimStep(x1-x2); +// return (T(1)-cond)*x1 + cond*x2; +// // return(CondExpLt(x1, x2, x1, x2)); +// } //double phi(double x); inline int nimStep(double x) { return(x >= 0 ? 1 : 0);} //template -inline T nimDerivs_nimStep(T x){ - return(CondExpGe(x, T(0), T(1), T(0))); -} +//inline T nimDerivs_nimStep(T x){ +// return(CondExpGe(x, T(0), T(1), T(0))); +//} inline double cube(double x) {return(x*x*x);} //template diff --git a/packages/nimble/inst/include/nimble/dists.h b/packages/nimble/inst/include/nimble/dists.h index d7c5ce0c9..f966d5eb7 100644 --- a/packages/nimble/inst/include/nimble/dists.h +++ b/packages/nimble/inst/include/nimble/dists.h @@ -48,6 +48,9 @@ extern "C" { // NIMBLE C wrappers called from R SEXP C_dmnorm_chol(SEXP, SEXP, SEXP, SEXP, SEXP); SEXP C_rmnorm_chol(SEXP, SEXP, SEXP); + SEXP C_PDinverse_logdet(SEXP); + SEXP C_dmnorm_inv_ld(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); + SEXP C_rmnorm_inv_ld(SEXP, SEXP, SEXP, SEXP); SEXP C_dmvt_chol(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); SEXP C_rmvt_chol(SEXP, SEXP, SEXP, SEXP); SEXP C_dlkj_corr_cholesky(SEXP, SEXP, SEXP, SEXP); @@ -95,6 +98,9 @@ void rdirch(double*, double*, int); double dmnorm_chol(double*, double*, double*, int, double, int, int); void rmnorm_chol(double *, double*, double*, int, double); +void PDinverse_logdet_internal(double *matPtr, double *ans, int n); +double dmnorm_inv_ld(double*, double*, double*, double*, int, int, int, int); +void rmnorm_inv_ld(double*, double*, double*, double*, int, int); double dmvt_chol(double*, double*, double*, double, int, double, int, int); void rmvt_chol(double *, double*, double*, double, int, double); double dwish_chol(double*, double*, double, int, double, int, int); diff --git a/packages/nimble/inst/include/nimble/nimDerivs_atomic_PDinverse_logdet.h b/packages/nimble/inst/include/nimble/nimDerivs_atomic_PDinverse_logdet.h new file mode 100644 index 000000000..81a1b6039 --- /dev/null +++ b/packages/nimble/inst/include/nimble/nimDerivs_atomic_PDinverse_logdet.h @@ -0,0 +1,86 @@ +#ifndef _NIMDERIVS_ATOMIC_PDINVERSE_LOGDET +#define _NIMDERIVS_ATOMIC_PDINVERSE_LOGDET + +#include +#include +#include +#include +#include +#include "nimDerivs_vecmat_utils.h" +#include "nimDerivs_atomic_matmult.h" +#include "nimDerivs_atomic_cache.h" + +void atomic_PDinverse_logdet(const MatrixXd_CppAD &x, + MatrixXd_CppAD &y); + +// MatrixXd_CppAD nimDerivs_PDinverse_logdet(const MatrixXd_CppAD &x); // not sure this is actually needed. +NimArr<1, CppAD::AD > nimDerivs_PDinverse_logdet(const NimArr<2, CppAD::AD > &x); + +class atomic_PDinverse_logdet_class : public CppAD::atomic_three, public nimble_atomic_base { + public: + atomic_PDinverse_logdet_class(const std::string& name); + private: + friend class atomic_cache_class; + friend class atomic_cache_class >; + atomic_cache_class double_cache; + atomic_cache_class > CppADdouble_cache; + + typedef EigenTemplateTypes::typeEigenConstMapStrd EigenConstMap; + typedef EigenTemplateTypes::typeEigenMapStrd EigenMap; + typedef EigenTemplateTypes >::typeEigenConstMapStrd metaEigenConstMap; + typedef EigenTemplateTypes >::typeEigenMapStrd metaEigenMap; + typedef EigenTemplateTypes>::typeMatrixXd metaEigenMatrixXd; + virtual bool for_type( + const CppAD::vector& parameter_x , + const CppAD::vector& type_x , + CppAD::vector& type_y ); + + virtual bool rev_depend( + const CppAD::vector& parameter_x , + const CppAD::vector& type_x , + CppAD::vector& depend_x , + const CppAD::vector& depend_y + ); + + virtual bool forward( + const CppAD::vector& parameter_x , + const CppAD::vector& type_x , + size_t need_y , + size_t order_low , + size_t order_up , + const CppAD::vector& taylor_x , + CppAD::vector& taylor_y ); + + virtual bool forward( + const CppAD::vector >& parameter_x , + const CppAD::vector& type_x , + size_t need_y , + size_t order_low , + size_t order_up , + const CppAD::vector >& taylor_x , + CppAD::vector >& taylor_y ); + + virtual bool reverse( + const CppAD::vector& parameter_x , + const CppAD::vector& type_x , + size_t order_up , + const CppAD::vector& taylor_x , + const CppAD::vector& taylor_y , + CppAD::vector& partial_x , + const CppAD::vector& partial_y ); + virtual bool reverse( + const CppAD::vector >& parameter_x , + const CppAD::vector& type_x , + size_t order_up , + const CppAD::vector >& taylor_x , + const CppAD::vector >& taylor_y , + CppAD::vector >& partial_x , + const CppAD::vector >& partial_y ); +}; + +atomic_PDinverse_logdet_class* new_atomic_PDinverse_logdet(void* tape_mgr, const std::string& name); +void delete_atomic_PDinverse_logdet(void* tape_mgr, atomic_PDinverse_logdet_class *atomic_PDinverse_logdet); + + + +#endif diff --git a/packages/nimble/inst/include/nimble/nimDerivs_atomic_classes.h b/packages/nimble/inst/include/nimble/nimDerivs_atomic_classes.h index 736e46f23..76e0e42bb 100644 --- a/packages/nimble/inst/include/nimble/nimDerivs_atomic_classes.h +++ b/packages/nimble/inst/include/nimble/nimDerivs_atomic_classes.h @@ -14,6 +14,7 @@ #include "nimDerivs_atomic_log_pow_int.h" #include "nimDerivs_atomic_matmult.h" #include "nimDerivs_atomic_matinverse.h" +#include "nimDerivs_atomic_PDinverse_logdet.h" #include "nimDerivs_atomic_backsolve.h" #include "nimDerivs_atomic_forwardsolve.h" #include "nimDerivs_atomic_cholesky.h" diff --git a/packages/nimble/inst/include/nimble/nimDerivs_atomic_zround.h b/packages/nimble/inst/include/nimble/nimDerivs_atomic_zround.h index 974b7bd68..6c9741e8b 100644 --- a/packages/nimble/inst/include/nimble/nimDerivs_atomic_zround.h +++ b/packages/nimble/inst/include/nimble/nimDerivs_atomic_zround.h @@ -144,5 +144,72 @@ class atomic_nimRound_class : public atomic_discrete_class nimDerivs_nimStep(const CppAD::AD x); + +atomic_nimStep_class *track_atomic_nimStep(void* tape_mgr_ptr, + std::vector* vec_ptr); + +class nimDerivs_nimStep_class { + public: + double operator()(double x); + double operator()(int x); + CppAD::AD operator()(const CppAD::AD &x); +}; + +class atomic_nimStep_class : public atomic_discrete_class { +public: + atomic_nimStep_class(const std::string& name); +}; + +#define TTT_ CppAD::AD + +inline TTT_ nimDerivs_pairmax(TTT_ x1, TTT_ x2) { + TTT_ cond = nimDerivs_nimStep(x1-x2); + return CppAD::azmul(cond, x1) + CppAD::azmul(TTT_(1) - cond, x2); +// return(CondExpGt(x1, x2, x1, x2)); +} + +inline TTT_ nimDerivs_pairmin(TTT_ x1, TTT_ x2) { + TTT_ cond = nimDerivs_nimStep(x1-x2); + return CppAD::azmul((TTT_(1)-cond), x1) + CppAD::azmul(cond, x2); +// return(CondExpLt(x1, x2, x1, x2)); +} + +inline TTT_ nimDerivs_nimEquals(TTT_ x1, TTT_ x2){ + TTT_ cond = nimDerivs_nimStep(x1-x2) * nimDerivs_nimStep(x2-x1); + return cond; +// return(CondExpEq(x1, x2, T(1), T(0))); +} + +inline TTT_ nimDerivs_CondExpEq(TTT_ x1, TTT_ x2, TTT_ y1, TTT_ y2) { + TTT_ cond = nimDerivs_nimEquals(x1, x2); + return CppAD::azmul(cond, y1) + CppAD::azmul((TTT_(1) - cond), y2); +// return(CondExpEq(x1, x2, T(1), T(0))); +} + +inline TTT_ nimDerivs_CondExpGe(TTT_ x1, TTT_ x2, TTT_ y1, TTT_ y2) { + TTT_ cond = nimDerivs_nimStep(x1-x2); + return CppAD::azmul(cond, y1) + CppAD::azmul((TTT_(1) - cond), y2); +} + +inline TTT_ nimDerivs_CondExpGt(TTT_ x1, TTT_ x2, TTT_ y1, TTT_ y2) { + TTT_ cond = 1-nimDerivs_nimStep(x2-x1); + return CppAD::azmul(cond, y1) + CppAD::azmul((TTT_(1) - cond), y2); +} + +inline TTT_ nimDerivs_CondExpLe(TTT_ x1, TTT_ x2, TTT_ y1, TTT_ y2) { + TTT_ cond = nimDerivs_nimStep(x2-x1); + return CppAD::azmul(cond, y1) + CppAD::azmul((TTT_(1) - cond), y2); +} + +inline TTT_ nimDerivs_CondExpLt(TTT_ x1, TTT_ x2, TTT_ y1, TTT_ y2) { + TTT_ cond = 1-nimDerivs_nimStep(x1-x2); + return CppAD::azmul(cond, y1) + CppAD::azmul((TTT_(1) - cond), y2); +} + +#undef TTT_ + #endif diff --git a/packages/nimble/inst/include/nimble/nimDerivs_dists.h b/packages/nimble/inst/include/nimble/nimDerivs_dists.h index e4b1e1538..7772db88c 100644 --- a/packages/nimble/inst/include/nimble/nimDerivs_dists.h +++ b/packages/nimble/inst/include/nimble/nimDerivs_dists.h @@ -19,7 +19,7 @@ /* because there was a bug using CPPAD_DISCRETE_FUNCTION from CppAD with base2ad. */ /* That bug was fixed, but after we worked around it anyway. */ inline double discrete_round(const double &x) -{ +{ double out_x = round(x); return(out_x); } @@ -37,11 +37,18 @@ CPPAD_DISCRETE_FUNCTION(double, discrete_round) /* dnorm: normal distribution */ /* This case is modified from TMB code so that give_log can be different when the tape is used from when it is recorded. */ + +template +Type log_or_exp(const Type &x, const Type &give_log) { + return CppAD::azmul(give_log, x) + CppAD::azmul((Type(1) - give_log), exp(x)); +} + template Type nimDerivs_dnorm(Type x, Type mean, Type sd, Type give_log) { Type res = -log(Type(sqrt(2*M_PI))*sd)-Type(.5)*pow((x-mean)/sd,2); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + //res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); return(res); } @@ -53,10 +60,82 @@ Type nimDerivs_dnorm_logFixed(Type x, Type mean, Type sd, int give_log) return(res); } +/* dmnorm: Multivariate normal distribution */ +/* Version that uses inv_ld, a vector of length n*n+1, with the precision elements followed by log determinant of the covariance.*/ +template +Type nimDerivs_nimArr_dmnorm_inv_ld(NimArr<1, Type> &x, + NimArr<1, Type> &mean, + NimArr<2, Type> &mat, + NimArr<1, Type> &inv_ld, + Type prec_param, // This gets baked in. It is essentially a compile time arg. + Type give_log, + Type overwrite_inputs) { + typedef Eigen::Matrix MatrixXt; + + int n = x.dimSize(0); + int i; + Type dens = Type(-n * M_LN_SQRT_2PI); + double pmhalf = CppAD::Value(prec_param) == 1 ? 0.5 : -0.5; + dens += Type(pmhalf) * inv_ld[n*n]; // log(det(cov)) term + /* Note that inv_ld[n*n] will be log(det(cov)) even if prec_param=TRUE */ + + MatrixXt xCopy(n, 1); + for(i = 0; i < n; i++) + xCopy(i, 0) = x[i] - mean[i]; + + NimArr<1, Type> prec_possible_copy; + NimArr<2, Type> mat_possible_copy; + Eigen::Map mapPrec(0,0,0); + if(CppAD::Value(prec_param) == 1) { // mat is the precision. + new (&mapPrec) Eigen::Map(nimArrCopyIfNeeded<2, Type>(mat, mat_possible_copy).getPtr(), n, n); + } else { // mat is the covariance so inv_ld has the precision flattened. + new (&mapPrec) Eigen::Map(nimArrCopyIfNeeded<1, Type>(inv_ld, prec_possible_copy).getPtr(), n, n); + } + + dens -= Type(0.5) * (xCopy.transpose() * mapPrec.template selfadjointView() * xCopy).sum(); // quadratic form term. sum() makes it a C++ scalar. + dens = log_or_exp(dens, give_log); +// dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); + return(dens); +} + +template +Type nimDerivs_nimArr_dmnorm_inv_ld_logFixed(NimArr<1, Type> &x, + NimArr<1, Type> &mean, + NimArr<2, Type> &mat, + NimArr<1, Type> &inv_ld, + Type prec_param, // This gets baked in. It is essentially a compile time arg. + int give_log, + Type overwrite_inputs) { + typedef Eigen::Matrix MatrixXt; + + int n = x.dimSize(0); + int i; + Type dens = Type(-n * M_LN_SQRT_2PI); + double pmhalf = CppAD::Value(prec_param) == 1 ? 0.5 : -0.5; + dens += Type(pmhalf) * inv_ld[n*n]; // log(det(cov)) term + + MatrixXt xCopy(n, 1); + for(i = 0; i < n; i++) + xCopy(i, 0) = x[i] - mean[i]; + + NimArr<1, Type> prec_possible_copy; + NimArr<2, Type> mat_possible_copy; + Eigen::Map mapPrec(0,0,0); + if(CppAD::Value(prec_param) == 1) { // mat is the precision. + new (&mapPrec) Eigen::Map(nimArrCopyIfNeeded<2, Type>(mat, mat_possible_copy).getPtr(), n, n); + } else { // mat is the covariance so inv_ld has the precision flattened. + new (&mapPrec) Eigen::Map(nimArrCopyIfNeeded<1, Type>(inv_ld, prec_possible_copy).getPtr(), n, n); + } + dens -= Type(0.5) * (xCopy.transpose() * mapPrec.template selfadjointView() * xCopy).sum(); // quadratic form term. sum() makes it a C++ scalar. + if(!give_log){ + dens = exp(dens); + } + return(dens); +} /* dmnorm: Multivariate normal distribution */ template -Type nimDerivs_nimArr_dmnorm_chol(NimArr<1, Type> &x, NimArr<1, Type> &mean, NimArr<2, Type> &chol, Type prec_param, Type give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dmnorm_chol(NimArr<1, Type> &x, NimArr<1, Type> &mean, NimArr<2, Type> &chol, Type prec_param, Type give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(prec_param)) std::cout<<"Warning: In dmnorm, prec_param with value = "< &x, NimArr<1, Type> &mean, Nim int i; Type dens = Type(-n * M_LN_SQRT_2PI); Type sumDens = Type(0.); - for(i = 0; i < n*n; i += n + 1) + for(i = 0; i < n*n; i += n + 1) sumDens += log(chol[i]); - dens += CppAD::CondExpEq(prec_param, Type(1), sumDens, -sumDens); + dens += (Type(2)*prec_param - Type(1)) * sumDens; +// dens += CppAD::CondExpEq(prec_param, Type(1), sumDens, -sumDens); MatrixXt xCopy(n, 1); for(i = 0; i < n; i++) @@ -97,7 +177,8 @@ Type nimDerivs_nimArr_dmnorm_chol(NimArr<1, Type> &x, NimArr<1, Type> &mean, Nim // presumably because of column-major order interacting well with the solve. xCopy = xCopy.array()*xCopy.array(); dens += -Type(0.5)*(xCopy.sum()); - dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); + dens = log_or_exp(dens, give_log); +// dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); return(dens); } @@ -107,15 +188,16 @@ Type nimDerivs_nimArr_dmnorm_chol_logFixed(NimArr<1, Type> &x, NimArr<1, Type> & if(!CppAD::Constant(prec_param)) std::cout<<"Warning: In dmnorm, prec_param with value = "< &x, NimArr<1, Type> & /* dmvt: Multivariate t distribution */ template -Type nimDerivs_nimArr_dmvt_chol(NimArr<1, Type> &x, NimArr<1, Type> &mu, NimArr<2, Type> &chol, Type df, Type prec_param, Type give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dmvt_chol(NimArr<1, Type> &x, NimArr<1, Type> &mu, NimArr<2, Type> &chol, Type df, Type prec_param, Type give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(prec_param)) std::cout<<"Warning: In dmvt, prec_param with value = "< &x, NimArr<1, Type> &mu, NimArr< /* mapChol.template triangularView()*xCopy, */ /* mapChol.template triangularView().transpose().solve(xCopy)); */ xCopy = xCopy.array()*xCopy.array(); - + dens += -Type(0.5)*(df + Type(n)) * log(Type(1.) + xCopy.sum() / df); - - dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); + dens = log_or_exp(dens, give_log); +// dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); return(dens); } template -Type nimDerivs_nimArr_dmvt_chol_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &mu, NimArr<2, Type> &chol, Type df, Type prec_param, int give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dmvt_chol_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &mu, NimArr<2, Type> &chol, Type df, Type prec_param, int give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(prec_param)) @@ -199,12 +282,12 @@ Type nimDerivs_nimArr_dmvt_chol_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &mu int n = x.dimSize(0); int i; - Type dens = nimDerivs_lgammafn((df + Type(n)) / Type(2.)) - nimDerivs_lgammafn(df / Type(2.)) - Type(n) * Type(M_LN_SQRT_PI) - Type(n) * log(df) / Type(2.); + Type dens = nimDerivs_lgammafn((df + Type(n)) / Type(2.)) - nimDerivs_lgammafn(df / Type(2.)) - Type(n) * Type(M_LN_SQRT_PI) - Type(n) * log(df) / Type(2.); Type sumDens = Type(0.); - for(i = 0; i < n*n; i += n + 1) + for(i = 0; i < n*n; i += n + 1) sumDens += log(chol[i]); - - dens += CppAD::CondExpEq(prec_param, Type(1), sumDens, -sumDens); + dens += (Type(2)*prec_param - Type(1)) * sumDens; + // dens += CppAD::CondExpEq(prec_param, Type(1), sumDens, -sumDens); MatrixXt xCopy(n, 1); for(i = 0; i < n; i++) @@ -223,9 +306,9 @@ Type nimDerivs_nimArr_dmvt_chol_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &mu /* mapChol.template triangularView()*xCopy, */ /* mapChol.template triangularView().transpose().solve(xCopy)); */ xCopy = xCopy.array()*xCopy.array(); - + dens += -Type(0.5)*(df + Type(n)) * log(Type(1.) + xCopy.sum() / df); - + if(!give_log){ dens = exp(dens); } @@ -259,13 +342,13 @@ Type nimDerivs_nimArr_dlkj_corr_cholesky(NimArr<2, Type> &x, Type eta, Type p, T */ Type dens = (mapX.diagonal().array().log() * (Type(p) - pseq + Type(2.0)*eta - Type(2.0))).sum(); - - dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); + dens = log_or_exp(dens, give_log); + // dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); return(dens); } /* Type sumDens = Type(0.); - for(i = 0; i < n*n; i += n + 1) + for(i = 0; i < n*n; i += n + 1) sumDens += log(chol[i]);*/ @@ -292,10 +375,10 @@ Type nimDerivs_nimArr_dlkj_corr_cholesky_logFixed(NimArr<2, Type> &x, Type eta, } return(dens); } - + /* dwish: Wishart distribution */ template -Type nimDerivs_nimArr_dwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, Type give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, Type give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(scale_param)) @@ -310,16 +393,17 @@ Type nimDerivs_nimArr_dwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type // Eigen::Map mapX(x.getPtr(), n, n); int p = x.dim()[0]; Type dens = (df * mapChol.diagonal().array().log()).sum(); - dens = CppAD::CondExpEq(scale_param, Type(1), -dens, dens); + dens = (Type(1) - Type(2) * scale_param) * dens; + // dens = CppAD::CondExpEq(scale_param, Type(1), -dens, dens); dens += -(df*p/Type(2.) * Type(M_LN2) + p*(p-Type(1.))*Type(M_LN_SQRT_PI/2.)); for(int i = 0; i < p; i++) dens -= nimDerivs_lgammafn((df - Type(i)) / Type(2.)); - + /* Lower may be more efficient: https://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html#details */ MatrixXt Lx = mapX.template selfadjointView().llt().matrixL(); dens += (df - p - Type(1.)) * Lx.diagonal().array().log().sum(); - + if(CppAD::Value(scale_param) == 1) { dens -= Type(0.5) * (mapChol.template triangularView().transpose().solve(Lx)).squaredNorm(); } else { @@ -329,19 +413,19 @@ Type nimDerivs_nimArr_dwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type /* dens -= Type(0.5) * CppAD::CondExpEq(scale_param, Type(1), */ /* mapChol.template triangularView().transpose().solve(Lx).squaredNorm(), */ /* (mapChol.template triangularView()*Lx).squaredNorm()); */ - - dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); + dens = log_or_exp(dens, give_log); + // dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); return(dens); } template -Type nimDerivs_nimArr_dwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, int give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, int give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(scale_param)) std::cout<<"Warning: In dwish, scale_param with value = "< chol_possible_copy; Eigen::Map mapChol(nimArrCopyIfNeeded<2, Type>(chol, chol_possible_copy).getPtr(), n, n); @@ -351,12 +435,13 @@ Type nimDerivs_nimArr_dwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> &c // Eigen::Map mapX(x.getPtr(), n, n); int p = x.dim()[0]; Type dens = (df * mapChol.diagonal().array().log()).sum(); - dens = CppAD::CondExpEq(scale_param, Type(1), -dens, dens); + dens = (Type(1) - Type(2) * scale_param) * dens; + // dens = CppAD::CondExpEq(scale_param, Type(1), -dens, dens); dens += -(df*p/Type(2.) * Type(M_LN2) + p*(p-Type(1.))*Type(M_LN_SQRT_PI/2.)); for(int i = 0; i < p; i++) dens -= nimDerivs_lgammafn((df - Type(i)) / Type(2.)); - + /* Lower may be more efficient: https://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html#details */ MatrixXt Lx = mapX.template selfadjointView().llt().matrixL(); dens += (df - p - Type(1.)) * Lx.diagonal().array().log().sum(); @@ -379,7 +464,7 @@ Type nimDerivs_nimArr_dwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> &c /* dinvwish: Inverse Wishart distribution */ template -Type nimDerivs_nimArr_dinvwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, Type give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dinvwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, Type give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(scale_param)) @@ -394,12 +479,13 @@ Type nimDerivs_nimArr_dinvwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, T int p = x.dim()[0]; Type dens = (df * mapChol.diagonal().array().log()).sum(); - dens = CppAD::CondExpEq(scale_param, Type(1), dens, -dens); + dens = (Type(2) * scale_param - Type(1)) * dens; + // dens = CppAD::CondExpEq(scale_param, Type(1), dens, -dens); dens += -(df*p/Type(2.) * Type(M_LN2) + p*(p-Type(1.))*Type(M_LN_SQRT_PI/2.)); for(int i = 0; i < p; i++) dens -= nimDerivs_lgammafn((df - Type(i)) / Type(2.)); - + /* Lower may be more efficient: https://eigen.tuxfamily.org/dox/classEigen_1_1LLT.html#details */ MatrixXt Lx = mapX.template selfadjointView().llt().matrixL(); dens -= (df + p + Type(1.)) * Lx.diagonal().array().log().sum(); @@ -416,13 +502,13 @@ Type nimDerivs_nimArr_dinvwish_chol(NimArr<2, Type> &x, NimArr<2, Type> &chol, T /* dens -= Type(0.5) * CppAD::CondExpEq(scale_param, Type(1), */ /* (Lx.template triangularView().solve(mapChol.transpose())).squaredNorm(), */ /* (Lx.template triangularView().solve(mapChol.template triangularView().solve(MatrixXd::Identity(n, n)))).squaredNorm()); */ - - dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); + dens = log_or_exp(dens, give_log); + //dens = CppAD::CondExpEq(give_log, Type(1), dens, exp(dens)); return(dens); } template -Type nimDerivs_nimArr_dinvwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, int give_log, Type overwrite_inputs) { +Type nimDerivs_nimArr_dinvwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> &chol, Type df, Type scale_param, int give_log, Type overwrite_inputs) { typedef Eigen::Matrix MatrixXt; if(!CppAD::Constant(scale_param)) @@ -437,7 +523,8 @@ Type nimDerivs_nimArr_dinvwish_chol_logFixed(NimArr<2, Type> &x, NimArr<2, Type> int p = x.dim()[0]; Type dens = (df * mapChol.diagonal().array().log()).sum(); - dens = CppAD::CondExpEq(scale_param, Type(1), dens, -dens); + dens = (Type(2) * scale_param - Type(1)) * dens; + // dens = CppAD::CondExpEq(scale_param, Type(1), dens, -dens); dens += -(df*p/Type(2.) * Type(M_LN2) + p*(p-Type(1.))*Type(M_LN_SQRT_PI/2.)); for(int i = 0; i < p; i++) @@ -469,7 +556,7 @@ template Type nimDerivs_nimArr_dmulti(const NimArr<1, Type> &x, const Type &size, const NimArr<1, Type> &prob, - const Type &give_log) + const Type &give_log) { Type sumProb = 0.0; Type sumX = 0.0; @@ -484,12 +571,14 @@ Type nimDerivs_nimArr_dmulti(const NimArr<1, Type> &x, for(int i = 0; i < K; i++) { logProb += nimDerivs_log_pow_int(prob[i]/sumProb, x[i]) - nimDerivs_lgammafn(x[i] + Type(1.)); } - logProb = CppAD::CondExpEq(sumX, size, + logProb = nimDerivs_CondExpEq(sumX, size, //CppAD::CondExpEq(sumX, size, logProb, -Type(std::numeric_limits::infinity())); - Type ans = CppAD::CondExpEq(give_log, Type(1), - logProb, - exp(logProb)); + + Type ans = log_or_exp(logProb, give_log); + //Type ans = CppAD::CondExpEq(give_log, Type(1), + // logProb, + // exp(logProb)); return ans; } @@ -497,7 +586,7 @@ template Type nimDerivs_nimArr_dmulti_logFixed(const NimArr<1, Type> &x, const Type &size, const NimArr<1, Type> &prob, - int give_log) + int give_log) { Type sumProb(0.0); Type sumX(0.0); @@ -512,7 +601,7 @@ Type nimDerivs_nimArr_dmulti_logFixed(const NimArr<1, Type> &x, for(int i = 0; i < K; i++) { logProb += nimDerivs_log_pow_int(prob[i]/sumProb, x[i]) - nimDerivs_lgammafn(x[i] + Type(1.)); } - logProb = CppAD::CondExpEq(sumX, size, + logProb = nimDerivs_CondExpEq(sumX, size, //CppAD::CondExpEq(sumX, size, logProb, -Type(std::numeric_limits::infinity())); if(give_log) { @@ -532,9 +621,10 @@ Type nimDerivs_nimArr_dcat(const Type &x, sumProb += prob[ j ]; } Type ansProb = stoch_ind_get(prob, x - 1) / sumProb; - Type ans = CppAD::CondExpEq(give_log, Type(1), - log(ansProb), - ansProb); + Type ans = log_or_exp(ansProb, give_log); + //Type ans = CppAD::CondExpEq(give_log, Type(1), + // log(ansProb), + // ansProb); return(ans); } @@ -556,7 +646,7 @@ Type nimDerivs_nimArr_dcat_logFixed(const Type &x, } // Some of the following functions seem to be missing Type() casting of constants. --CP 2020-04-22 - + /* dt: Student's t distribution */ /* This case is modified from TMB code so that give_log can be different when the tape is used from when it is recorded. @@ -565,7 +655,8 @@ template Type nimDerivs_dt(Type x, Type df, Type give_log) { Type res = nimDerivs_lgammafn((df+1)/2) - Type(1)/2*log(df*M_PI) -nimDerivs_lgammafn(df/2) - (df+1)/2*log(1+x*x/df); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -584,7 +675,8 @@ template Type nimDerivs_dt_nonstandard(Type x, Type df, Type mu, Type sigma, Type give_log) { Type res = nimDerivs_dt_logFixed( (x - mu)/sigma, df, 1) - log(sigma); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -604,7 +696,8 @@ template Type nimDerivs_dgamma(Type y, Type shape, Type scale, Type give_log) { Type res = -nimDerivs_lgammafn(shape)+(shape-Type(1.0))*log(y)-y/scale-shape*log(scale); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -625,7 +718,8 @@ Type nimDerivs_dinvgamma(Type x, Type shape, Type rate, Type give_log) { Type xinv = Type(1.0)/x; Type res = nimDerivs_dgamma_logFixed(xinv, shape, rate, 1) - 2*log(x); - res = CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -644,7 +738,8 @@ template Type nimDerivs_dsqrtinvgamma(Type x, Type shape, Type rate, Type give_log) { Type res = nimDerivs_dinvgamma_logFixed(x*x, shape, rate, 1) + log(2*x); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -671,8 +766,9 @@ inline Type nimDerivs_dnbinom(const Type &x, const Type &size, const Type &prob, // but x is integer and p may be == 1, so nimDerivs_log_pow_int(Type(1.)-p, x) // is needed instead of x * log(1.-p) Type res = nimDerivs_lgammafn(x+n)-nimDerivs_lgammafn(n)-nimDerivs_lgammafn(x+Type(1))+ - n*log(p) + nimDerivs_log_pow_int(Type(1.)-p, x); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + n*log(p) + nimDerivs_log_pow_int(Type(1.)-p, x); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -683,7 +779,7 @@ inline Type nimDerivs_dnbinom_logFixed(const Type &x, const Type &size, const Ty Type n=size; Type p=prob; Type res = nimDerivs_lgammafn(x+n)-nimDerivs_lgammafn(n)-nimDerivs_lgammafn(x+Type(1))+ - n*log(p) + nimDerivs_log_pow_int(Type(1.)-p, x); + n*log(p) + nimDerivs_log_pow_int(Type(1.)-p, x); if(!give_log){ res = exp(res); } @@ -696,7 +792,8 @@ template inline Type nimDerivs_dpois(const Type &x, const Type &lambda, Type give_log) { Type res = -lambda + nimDerivs_log_pow_int(lambda, x) - nimDerivs_lgammafn(x+Type(1)); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -712,26 +809,27 @@ inline Type nimDerivs_dpois_logFixed(const Type &x, const Type &lambda, int give /* dexp: exponential distribution */ /* Parameterization available in nimble. Code modified from TMB*/ -template +template Type nimDerivs_dexp_nimble(Type x, Type rate, Type give_log) { Type res = log(rate)-rate*x; - res = CppAD::CondExpEq(give_log, Type(1), - res, - exp(res)); - /* Type res = CppAD::CondExpEq(give_log, Type(0), */ - /* rate*exp(-rate*x), */ - /* log(rate)-rate*x); */ + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), + // res, + // exp(res)); + /* Type res = CppAD::CondExpEq(give_log, Type(0), + /* rate*exp(-rate*x), + /* log(rate)-rate*x); */ return(res); } -template +template Type nimDerivs_dexp_nimble_logFixed(Type x, Type rate, int give_log) { Type res; if(!give_log){ res = rate*exp(-rate*x); - } + } else{ res = log(rate)-rate*x; } @@ -745,8 +843,9 @@ template Type nimDerivs_dexp(Type x, Type scale, Type give_log) { Type res = -log(scale)-x/scale; - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); - // Type res = CppAD::CondExpEq(give_log, Type(0), exp(-x/scale)/scale, -log(scale)-x/scale); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + // Type res = CppAD::CondExpEq(give_log, Type(0), exp(-x/scale)/scale, -log(scale)-x/scale); return(res); } @@ -767,7 +866,8 @@ Type nimDerivs_dexp_logFixed(Type x, Type scale, int give_log) template Type nimDerivs_dflat(Type x, Type give_log) { - Type res = CppAD::CondExpEq(give_log, Type(1), Type(0), Type(1)); + Type res = Type(1) - give_log; //* Type(0) + (Type(1) - give_log) * Type(1); + //Type res = CppAD::CondExpEq(give_log, Type(1), Type(0), Type(1)); return(res); } template @@ -779,8 +879,10 @@ Type nimDerivs_dflat_logFixed(Type x, int give_log) template Type nimDerivs_dhalfflat(Type x, Type give_log) { - Type res = CppAD::CondExpGe(x, Type(0), Type(0), -CppAD::numeric_limits::max()); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + Type res = nimDerivs_CondExpGe(x, Type(0), Type(0), -CppAD::numeric_limits::max()); + // Type res = CppAD::CondExpGe(x, Type(0), Type(0), -CppAD::numeric_limits::max()); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } template @@ -788,34 +890,41 @@ Type nimDerivs_dhalfflat_logFixed(Type x, int give_log) { Type res; if(give_log) { - res = CppAD::CondExpGe(x, Type(0), Type(0), -CppAD::numeric_limits::max()); + res = nimDerivs_CondExpGe(x, Type(0), Type(0), -CppAD::numeric_limits::max()); + // res = CppAD::CondExpGe(x, Type(0), Type(0), -CppAD::numeric_limits::max()); } else { - res = CppAD::CondExpGe(x, Type(0), Type(1), Type(0)); + res = nimDerivs_nimStep(x); + // res = CppAD::CondExpGe(x, Type(0), Type(1), Type(0)); } return(res); } /* dunif: uniform distribution */ /* TMB does not have these cases. */ -template +template Type nimDerivs_dunif(Type x, Type a, Type b, Type give_log) -{ +{ Type res = Type(1.0)/(b-a); - res = CppAD::CondExpEq(give_log, Type(1), log(res), res); - Type ans = CppAD::CondExpLe(x, b, res, -CppAD::numeric_limits::max()); - ans = CppAD::CondExpGe(x, a, ans, -CppAD::numeric_limits::max()); + res = CppAD::azmul(give_log, log(res)) + CppAD::azmul(Type(1)-give_log, res); +// res = CppAD::CondExpEq(give_log, Type(1), log(res), res); + Type ans = nimDerivs_CondExpLe(x, b, res, -CppAD::numeric_limits::max()); +// Type ans = CppAD::CondExpLe(x, b, res, -CppAD::numeric_limits::max()); + ans = nimDerivs_CondExpGe(x, a, ans, -CppAD::numeric_limits::max()); +// ans = CppAD::CondExpGe(x, a, ans, -CppAD::numeric_limits::max()); return(ans); } -template +template Type nimDerivs_dunif_logFixed(Type x, Type a, Type b, int give_log) -{ +{ Type res = Type(1.)/(b-a); if(give_log){ res = log(res); } - Type ans = CppAD::CondExpLe(x, b, res, -CppAD::numeric_limits::max()); - ans = CppAD::CondExpGe(x, a, ans, -CppAD::numeric_limits::max()); + Type ans = nimDerivs_CondExpLe(x, b, res, -CppAD::numeric_limits::max()); + // Type ans = CppAD::CondExpLe(x, b, res, -CppAD::numeric_limits::max()); + ans = nimDerivs_CondExpGe(x, a, ans, -CppAD::numeric_limits::max()); + // ans = CppAD::CondExpGe(x, a, ans, -CppAD::numeric_limits::max()); return(ans); } @@ -823,16 +932,17 @@ Type nimDerivs_dunif_logFixed(Type x, Type a, Type b, int give_log) /* This case is modified from TMB code so that give_log can be different when the tape is used from when it is recorded. */ /* Also the conditional check on x > 0 has been removed. */ -template +template Type nimDerivs_dweibull(Type x, Type shape, Type scale, Type give_log) { Type res = shape/scale * pow(x/scale,shape-1) * exp(-pow(x/scale,shape)); - res = CppAD::CondExpEq(give_log, Type(0), res, log(res)); + res = CppAD::azmul(give_log, log(res)) + CppAD::azmul(Type(1)-give_log, res); + // res = CppAD::CondExpEq(give_log, Type(0), res, log(res)); return(res); } /* This case is modified from TMB to remove the conditional check on x > 0. */ -template +template Type nimDerivs_dweibull_logFixed(Type x, Type shape, Type scale, int give_log) { Type res = shape/scale * pow(x/scale,shape-1) * exp(-pow(x/scale,shape)); @@ -853,7 +963,8 @@ Type nimDerivs_dbinom(Type k, Type size, Type prob, Type give_log) { Type res = nimDerivs_lgammafn(size+1)-nimDerivs_lgammafn(k+1)-nimDerivs_lgammafn(size-k+1)+ CppAD::azmul(k, log(prob))+ CppAD::azmul(size-k, log(1-prob)); - res = CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -872,15 +983,16 @@ Type nimDerivs_dbinom_logFixed(Type k, Type size, Type prob, int give_log) /* TMB does not appear to have this distribution. */ template Type nimDerivs_dchisq(Type x, Type df, Type give_log) -{ +{ Type res = (df/Type(2.0) - Type(1.0))*log(x) - (x/Type(2.0)) - (df/Type(2.0))*log(Type(2.0)) - nimDerivs_lgammafn(df/Type(2.0)); - res = CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CondExpEq(give_log, Type(1), res, exp(res)); return(res); } template Type nimDerivs_dchisq_logFixed(Type x, Type df, int give_log) -{ +{ Type res = (df/Type(2.0) - Type(1.0))*log(x) - (x/Type(2.0)) - (df/Type(2.0))*log(Type(2.0)) - nimDerivs_lgammafn(df/Type(2.0)); if(!give_log){ res = exp(res); @@ -893,9 +1005,10 @@ Type nimDerivs_dchisq_logFixed(Type x, Type df, int give_log) template Type nimDerivs_dbeta(Type x, Type shape1, Type shape2, Type give_log) { - Type res = nimDerivs_lgammafn(shape1+shape2) - nimDerivs_lgammafn(shape1) - nimDerivs_lgammafn(shape2) + + Type res = nimDerivs_lgammafn(shape1+shape2) - nimDerivs_lgammafn(shape1) - nimDerivs_lgammafn(shape2) + (shape1-1)*log(x) + (shape2-1)*log(1-x); - res = CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CondExpEq(give_log, Type(1), res, exp(res)); /* The following code, when used with tape optimize()ation, gives a crash when the wrong bracnch of CondExpEq is triggered */ /* Type res = CondExpEq(give_log, Type(0), */ @@ -911,13 +1024,13 @@ Type nimDerivs_dbeta_logFixed(Type x, Type shape1, Type shape2, int give_log) { Type res; if(give_log){ - res = nimDerivs_lgammafn(shape1+shape2) - nimDerivs_lgammafn(shape1) - nimDerivs_lgammafn(shape2) + + res = nimDerivs_lgammafn(shape1+shape2) - nimDerivs_lgammafn(shape1) - nimDerivs_lgammafn(shape2) + (shape1-1)*log(x) + (shape2-1)*log(1-x); } else{ - res = exp(nimDerivs_lgammafn(shape1+shape2) - nimDerivs_lgammafn(shape1) - + res = exp(nimDerivs_lgammafn(shape1+shape2) - nimDerivs_lgammafn(shape1) - nimDerivs_lgammafn(shape2)) * pow(x,shape1-1) * pow(1-x,shape2-1); - } + } return(res); } @@ -928,7 +1041,7 @@ template Type nimDerivs_dlogis(Type x, Type location, Type scale, Type give_log) { Type res = -(x-location)/scale - log(scale) - 2*log(1+exp(-(x-location)/scale)); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); return(res); } @@ -945,15 +1058,15 @@ Type nimDerivs_dlogis_logFixed(Type x, Type location, Type scale, int give_log) template Type nimDerivs_dlnorm(Type x, Type mean, Type sd, Type give_log) { - Type res = -log(x) - log(sd) - Type(0.5)*log(Type(2.0*M_PI)) - Type(.5)*pow((log(x)-mean)/sd, 2); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + Type res = -log(x) - log(sd) - Type(0.5)*log(Type(2.0*M_PI)) - Type(.5)*pow((log(x)-mean)/sd, 2); + res = log_or_exp(res, give_log); return(res); } template Type nimDerivs_dlnorm_logFixed(Type x, Type mean, Type sd, int give_log) { - Type res = -log(x) - log(sd) - Type(0.5)*log(Type(2.0*M_PI)) - Type(.5)*pow((log(x)-mean)/sd, 2); + Type res = -log(x) - log(sd) - Type(0.5)*log(Type(2.0*M_PI)) - Type(.5)*pow((log(x)-mean)/sd, 2); if(!give_log){ res = exp(res); } @@ -967,7 +1080,8 @@ Type nimDerivs_ddexp(Type x, Type location, Type scale, Type give_log) Type res = log(Type(0.5)) - log(scale) - abs(x-location)/scale; // Type res = (Type(1.0/2.0) / scale) * exp(-abs(x - location) / scale); // res = CppAD::CondExpEq(give_log, Type(1), log(res), res); - res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); + res = log_or_exp(res, give_log); + // res = CppAD::CondExpEq(give_log, Type(1), res, exp(res)); return(res); } @@ -989,26 +1103,30 @@ Type nimDerivs_nimArr_ddirch(NimArr<1, Type> &x, NimArr<1, Type> &alpha, Type gi { Type K = alpha.size(); Type n = x.size(); - Type logres = CppAD::CondExpLt(K, Type(1), CppAD::numeric_limits::quiet_NaN(), Type(0)); - logres = CppAD::CondExpEq(K, n, logres, CppAD::numeric_limits::quiet_NaN()); - + Type logres = nimDerivs_CondExpLt(K, Type(1), CppAD::numeric_limits::quiet_NaN(), Type(0)); + //Type logres = CppAD::CondExpLt(K, Type(1), CppAD::numeric_limits::quiet_NaN(), Type(0)); + logres = nimDerivs_CondExpEq(K, n, logres, CppAD::numeric_limits::quiet_NaN()); + // logres = CppAD::CondExpEq(K, n, logres, CppAD::numeric_limits::quiet_NaN()); + Type sumAlpha = Type(0.0); Type sumX = Type(0.0); Type dens = Type(0.0); - + for(int i = 0; i < K; i++) { - logres += CppAD::CondExpGt(alpha[i], Type(0), Type(0), CppAD::numeric_limits::quiet_NaN()); - logres += CppAD::CondExpGe(x[i], Type(0), Type(0), -Type(std::numeric_limits::infinity())); - logres += CppAD::CondExpLe(x[i], Type(1), Type(0), -Type(std::numeric_limits::infinity())); + logres += nimDerivs_CondExpGt(alpha[i], Type(0), Type(0), CppAD::numeric_limits::quiet_NaN()); + // logres += CppAD::CondExpGt(alpha[i], Type(0), Type(0), CppAD::numeric_limits::quiet_NaN()); + logres += nimDerivs_CondExpGe(x[i], Type(0), Type(0), -Type(std::numeric_limits::infinity())); + logres += nimDerivs_CondExpLe(x[i], Type(1), Type(0), -Type(std::numeric_limits::infinity())); logres += (alpha[i]-Type(1)) * log(x[i]) - nimDerivs_lgammafn(alpha[i]) ; sumAlpha += alpha[i]; sumX += x[i]; } - logres += CppAD::CondExpLe(sumX, Type(1.0 + 10.0*std::numeric_limits::epsilon()), - CppAD::CondExpGe(sumX, Type(1.0 - 10.0*std::numeric_limits::epsilon()), Type(0), -Type(std::numeric_limits::infinity())), + logres += nimDerivs_CondExpLe(sumX, Type(1.0 + 10.0*std::numeric_limits::epsilon()), + nimDerivs_CondExpGe(sumX, Type(1.0 - 10.0*std::numeric_limits::epsilon()), Type(0), -Type(std::numeric_limits::infinity())), -Type(std::numeric_limits::infinity())); logres += nimDerivs_lgammafn(sumAlpha); - logres = CppAD::CondExpEq(give_log, Type(1), logres, exp(logres)); + logres = log_or_exp(logres, give_log); + // logres = nimDerivs_CondExpEq(give_log, Type(1), logres, exp(logres)); return(logres); } @@ -1017,23 +1135,23 @@ Type nimDerivs_nimArr_ddirch_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &alpha { Type K = alpha.size(); Type n = x.size(); - Type logres = CppAD::CondExpLt(K, Type(1), CppAD::numeric_limits::quiet_NaN(), Type(0)); - logres = CppAD::CondExpEq(K, n, logres, CppAD::numeric_limits::quiet_NaN()); - + Type logres = nimDerivs_CondExpLt(K, Type(1), CppAD::numeric_limits::quiet_NaN(), Type(0)); + logres = nimDerivs_CondExpEq(K, n, logres, CppAD::numeric_limits::quiet_NaN()); + Type sumAlpha = Type(0.0); Type sumX = Type(0.0); Type dens = Type(0.0); - + for(int i = 0; i < K; i++) { - logres += CppAD::CondExpGt(alpha[i], Type(0), Type(0), CppAD::numeric_limits::quiet_NaN()); - logres += CppAD::CondExpGe(x[i], Type(0), Type(0), -Type(std::numeric_limits::infinity())); - logres += CppAD::CondExpLe(x[i], Type(1), Type(0), -Type(std::numeric_limits::infinity())); + logres += nimDerivs_CondExpGt(alpha[i], Type(0), Type(0), CppAD::numeric_limits::quiet_NaN()); + logres += nimDerivs_CondExpGe(x[i], Type(0), Type(0), -Type(std::numeric_limits::infinity())); + logres += nimDerivs_CondExpLe(x[i], Type(1), Type(0), -Type(std::numeric_limits::infinity())); logres += (alpha[i]-Type(1)) * log(x[i]) - nimDerivs_lgammafn(alpha[i]) ; sumAlpha += alpha[i]; sumX += x[i]; } - logres += CppAD::CondExpLe(sumX, Type(1.0 + 10.0*std::numeric_limits::epsilon()), - CppAD::CondExpGe(sumX, Type(1.0 - 10.0*std::numeric_limits::epsilon()), Type(0), -Type(std::numeric_limits::infinity())), + logres += nimDerivs_CondExpLe(sumX, Type(1.0 + 10.0*std::numeric_limits::epsilon()), + nimDerivs_CondExpGe(sumX, Type(1.0 - 10.0*std::numeric_limits::epsilon()), Type(0), -Type(std::numeric_limits::infinity())), -Type(std::numeric_limits::infinity())); logres += nimDerivs_lgammafn(sumAlpha); if(!give_log){ @@ -1042,12 +1160,12 @@ Type nimDerivs_nimArr_ddirch_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &alpha return(logres); } -//Type nimDerivs_nimArr_dcar_normal(NimArr<1, Type> &x, double* adj, NimArr<1, Type> &weights, double* num, Type tau, int c, int zero_mean, Type give_log) +//Type nimDerivs_nimArr_dcar_normal(NimArr<1, Type> &x, double* adj, NimArr<1, Type> &weights, double* num, Type tau, int c, int zero_mean, Type give_log) // Should adj,weights,num,c,zero_mean be `Type`? And should they be `const`? // do we need `const Type &adj` because used as an index? //formerly did (int) adj[count] template -Type nimDerivs_nimArr_dcar_normal(NimArr<1, Type> &x, NimArr<1, Type> &adj, NimArr<1, Type> &weights, NimArr<1, Type> &num, Type tau, Type c, Type zero_mean, Type give_log) +Type nimDerivs_nimArr_dcar_normal(NimArr<1, Type> &x, NimArr<1, Type> &adj, NimArr<1, Type> &weights, NimArr<1, Type> &num, Type tau, Type c, Type zero_mean, Type give_log) { // This method implements the following density calculation: // p(x1, ..., xn, tau) = (tau/2/pi)^((N-c)/2) * exp(-tau/2 * sum_{i != j) w_ij (xi-xj)^2 ), @@ -1059,7 +1177,7 @@ Type nimDerivs_nimArr_dcar_normal(NimArr<1, Type> &x, NimArr<1, Type> &adj, NimA //PRINTF("c is equal to %d\n", c); int N = x.size(); int L = adj.size(); - + Type lp = 0; int count = 0; Type xi, xj; @@ -1077,13 +1195,13 @@ Type nimDerivs_nimArr_dcar_normal(NimArr<1, Type> &x, NimArr<1, Type> &adj, NimA lp *= Type(1/2.0); // accounts for double-summing over all (xi,xj) pairs lp *= Type(-1/2.0) * tau; lp += (Type(N)-c)/Type(2.0) * (log(tau) - Type(M_LN_2PI)); - - lp = CppAD::CondExpEq(give_log, Type(1), lp, exp(lp)); + lp = log_or_exp(lp, give_log); + // lp = CppAD::CondExpEq(give_log, Type(1), lp, exp(lp)); return(lp); } template -Type nimDerivs_nimArr_dcar_normal_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &adj, NimArr<1, Type> &weights, NimArr<1, Type> &num, Type tau, Type c, Type zero_mean, int give_log) +Type nimDerivs_nimArr_dcar_normal_logFixed(NimArr<1, Type> &x, NimArr<1, Type> &adj, NimArr<1, Type> &weights, NimArr<1, Type> &num, Type tau, Type c, Type zero_mean, int give_log) { // This method implements the following density calculation: // p(x1, ..., xn, tau) = (tau/2/pi)^((N-c)/2) * exp(-tau/2 * sum_{i != j) w_ij (xi-xj)^2 ), @@ -1157,7 +1275,8 @@ Type nimDerivs_nimArr_dcar_proper(NimArr<1, Type> &x, NimArr<1, Type> &mu, NimAr lp += (log(Type(1) - gamma*evs[i]) - log(M[i])) / Type(2.0); } lp += Type(N / 2.0) * (log(tau) - Type(M_LN_2PI)); - lp = CppAD::CondExpEq(give_log, Type(1), lp, exp(lp)); + lp = log_or_exp(lp, give_log); + // lp = CppAD::CondExpEq(give_log, Type(1), lp, exp(lp)); return(lp); } @@ -1219,9 +1338,9 @@ Type nimDerivs_nimArr_dcar_proper_logFixed(NimArr<1, Type> &x, NimArr<1, Type> & #define Type CppAD::AD inline Type nimDerivs_pow(Type x, Type y) { - /* We experimented with CppAD conditionals but still had cases that + /* We experimented with CppAD conditionals but still had cases that crashed or had constants baked into tapes when they shouldn't be - As a result we implemented nimDerivs_pow_int for cases where + As a result we implemented nimDerivs_pow_int for cases where y will only have integer values and never need its derivatives taken. */ @@ -1231,9 +1350,9 @@ inline Type nimDerivs_pow(Type x, Type y) { /* CppAD::pow(x, CppAD::Integer(y)), */ /* Type(CppAD::numeric_limits::quiet_NaN()))); */ - /* There is still a need to identify hard-coded cases with integer y and + /* There is still a need to identify hard-coded cases with integer y and divert them to nimDerivs_pow_int. Currently nimble-generated code - does not cast the second argument to Type() for pow, so one of the + does not cast the second argument to Type() for pow, so one of the over-loaded cases below (for y being int or double) will be called.*/ Type outVal = CppAD::pow(x, y); return(outVal); diff --git a/packages/nimble/inst/include/nimble/nimDists.h b/packages/nimble/inst/include/nimble/nimDists.h index 1e098108e..b0c0fc7ef 100644 --- a/packages/nimble/inst/include/nimble/nimDists.h +++ b/packages/nimble/inst/include/nimble/nimDists.h @@ -51,6 +51,13 @@ double nimArr_rcat(NimArr<1, double> &prob); double nimArr_dmnorm_chol(NimArr<1, double> &x, NimArr<1, double> &mean, NimArr<2, double> &chol, double prec_param, int give_log, int overwrite_inputs); void nimArr_rmnorm_chol(NimArr<1, double> &ans, NimArr<1, double> &mean, NimArr<2, double> &chol, double prec_param); +NimArr<1, double> PDinverse_logdet(NimArr<2, double> &mat); + +double nimArr_dmnorm_inv_ld(NimArr<1, double> &x, NimArr<1, double> &mean, + NimArr<2, double> &mat, NimArr<1, double> &inv_ld, int prec_param, int give_log, int overwrite_inputs); +void nimArr_rmnorm_inv_ld(NimArr<1, double> &ans, NimArr<1, double> &mean, + NimArr<2, double> &mat, NimArr<1, double> &inv_ld, int prec_param); + double nimArr_dmvt_chol(NimArr<1, double> &x, NimArr<1, double> &mean, NimArr<2, double> &chol, double df, double prec_param, int give_log, int overwrite_inputs); void nimArr_rmvt_chol(NimArr<1, double> &ans, NimArr<1, double> &mean, NimArr<2, double> &chol, double df, double prec_param); diff --git a/packages/nimble/inst/include/nimble/nimOptim.h b/packages/nimble/inst/include/nimble/nimOptim.h index 9ce540f30..260da407d 100644 --- a/packages/nimble/inst/include/nimble/nimOptim.h +++ b/packages/nimble/inst/include/nimble/nimOptim.h @@ -3,17 +3,17 @@ * Copyright (C) 2014-2017 Perry de Valpine, Christopher Paciorek, * Daniel Turek, Clifford Anderson-Bergman, Nick Michaud, Fritz Obermeyer, * Duncan Temple Lang. - * + * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. - * + * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. - * + * * You should have received a copy of the GNU General Public License * along with this program; if not, a copy is available at * https://www.R-project.org/Licenses/ @@ -246,7 +246,7 @@ class NimOptimProblem_Fun_Grad_Hess : public NimOptimProblem_Fun_Grad { double* parscale = NimOptimProblem::working_parscale.getPtr(); for(int i = 0; i < n; ++i) { for(int j = 0; j < n; ++j) { - NimOptimProblem::ans_hessian_[i,j] *= parscale[i] * parscale[j]; + NimOptimProblem::ans_hessian_(i,j) *= parscale[i] * parscale[j]; } } } diff --git a/packages/nimble/inst/include/nimble/nimbleCppAD.h b/packages/nimble/inst/include/nimble/nimbleCppAD.h index b8418db43..cbe572d7e 100644 --- a/packages/nimble/inst/include/nimble/nimbleCppAD.h +++ b/packages/nimble/inst/include/nimble/nimbleCppAD.h @@ -130,6 +130,7 @@ class atomic_forwardsolve_class; class atomic_cholesky_class; class atomic_matmult_class; class atomic_matinverse_class; +class atomic_PDinverse_logdet_class; class atomic_zround_class; class nimDerivs_floor_class; template class atomic_discrete_class; @@ -137,6 +138,7 @@ class atomic_floor_class; class atomic_ceil_class; class atomic_ftrunc_class; class atomic_nimRound_class; +class atomic_nimStep_class; class atomic_log_pow_int_class; class atomic_zb_over_a_class; class atomic_probit_class; @@ -297,6 +299,12 @@ class nimble_CppAD_tape_mgr { void delete_atomic_nimRound(atomic_nimRound_class *atomic_nimRound); atomic_nimRound_class *get_atomic_nimRound(std::vector* vec_ptr); + int nimStep_index; + bool nimStep_exists; + atomic_nimStep_class* new_atomic_nimStep(const std::string& name); + void delete_atomic_nimStep(atomic_nimStep_class *atomic_nimStep); + atomic_nimStep_class *get_atomic_nimStep(std::vector* vec_ptr); + int log_pow_int_index; bool log_pow_int_exists; atomic_log_pow_int_class* new_atomic_log_pow_int(const std::string& name); @@ -344,6 +352,9 @@ class nimble_CppAD_tape_mgr { atomic_matinverse_class* new_atomic_matinverse(const std::string& name); void delete_atomic_matinverse(atomic_matinverse_class *atomic_matinverse); + atomic_PDinverse_logdet_class* new_atomic_PDinverse_logdet(const std::string& name); + void delete_atomic_PDinverse_logdet(atomic_PDinverse_logdet_class *atomic_PDinverse_logdet); + std::vector > dummyOutputs; void add_dummyOutput(CppAD::AD &dummy); void sum_dummyOutputs_to_dependentVars(std::vector > &depVars); @@ -564,6 +575,7 @@ NimArr<1, double> make_vector_if_necessary(int); NimArr<1, double> make_vector_if_necessary(double); NimArr<1, double> make_vector_if_necessary(NimArr<1, double>); NimArr<1, double> make_vector_if_necessary(NimArr<1, int>); +NimArr<1, CppAD::AD > make_vector_if_necessary(CppAD::AD x); template NimArr<1, double> make_vector_if_necessary(const T &x) { // This case catches eigen ops diff --git a/packages/nimble/inst/include/nimble/nimbleCppADbaseClass.cpp b/packages/nimble/inst/include/nimble/nimbleCppADbaseClass.cpp index 1bae2c2e9..403944fad 100644 --- a/packages/nimble/inst/include/nimble/nimbleCppADbaseClass.cpp +++ b/packages/nimble/inst/include/nimble/nimbleCppADbaseClass.cpp @@ -17,8 +17,7 @@ void setOrdersFound(const NimArr<1, double> &derivOrders, int orderSize = derivOrders.size(); double const* array_derivOrders = derivOrders.getConstPtr(); maxOrder = 0; - maxOrder = - *std::max_element(array_derivOrders, array_derivOrders + orderSize); + maxOrder = *std::max_element(array_derivOrders, array_derivOrders + orderSize); std::fill(ordersFound, ordersFound + 3, false); for (int i = 0; i < orderSize; i++) { if ((array_derivOrders[i] > 2) | (array_derivOrders[i] < 0)) { @@ -41,186 +40,539 @@ bool check_inf_nan_gdi(CppAD::AD v) { return false; } +template +inline double gd_getValue(const BASE& x) { + return x; +} + +template<> +inline double gd_getValue(const CppAD::AD& x) { + return CppAD::Value(x); +} + template void getDerivs_internal(vector &independentVars, - TAPETYPE *ADtape, - const NimArr<1, double> &derivOrders, - const NimArr<1, double> &wrtVector, - nimSmartPtr &ansList) { - // std::cout<<"entering getDerivs_internal"< &derivOrders, + const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, BASE> &inDir, + const NimArr<1, BASE> &outDir, + nimSmartPtr &ansList) { + // std::cout<<"entering getDerivs_internal"<Range(); + + int maxOrder; + bool ordersFound[3]; + setOrdersFound(derivOrders, ordersFound, maxOrder); + + // std::cout<<"orders: "< 0 && outInds[0] != -1; + bool inDir_provided = inDir.size() > 0 && !ISNA(gd_getValue(inDir[0])); + bool outDir_provided = outDir.size() > 0 && !ISNA(gd_getValue(outDir[0])); + bool wrt_provided = wrtVector.size() > 0 && wrtVector[0] != -1; + //std::cout<<"outInds_provided: "<jacobian.setSize(q, wrt_n, false, false); + if(outDir_provided){ + if(outDir.size() != m) { + std::cout << "outDir must have the same length as output." << std::endl; + } } - if (ordersFound[2]) { - ansList->hessian.setSize(wrt_n, wrt_n, q, false, false); + + std::size_t length_wrt = wrtVector.size(); // dim of wrt vars + if(length_wrt == 2){ + if(wrtVector[1] == -1){ // 2nd element -1 is a filler to ensure there is a vector out of compilation + // The other input vectors were added at a later stage and don't need this kluge. + length_wrt = 1; + } } - vector cppad_derivOut; - std::vector w(q, 0); - for (size_t dy_ind = 0; dy_ind < q; dy_ind++) { - w[dy_ind] = 1; - if (maxOrder == 1) { - if(!infIndicators[dy_ind]){ -#ifdef _TIME_AD_GENERAL - derivs_run_tape_timer_start(); -#endif - cppad_derivOut = ADtape->Reverse(1, w); -#ifdef _TIME_AD_GENERAL - derivs_run_tape_timer_stop(); -#endif - } + std::size_t length_outInds = outInds.size(); // dim of outInds + + // defaults are for all of wrt, outInds, inDir and outDir provided EMPTY + size_t res_dimx_o1(n), res_dimx1_o2(n), res_dimx2_o2(n); + size_t res_dimy_o0(m), res_dimy_o1(m), res_dimy_o2(m); + // N.B. The "All" flags may not be checked when using inDir or outDir, but they are set for logical coherence anyway. + bool wrtAllx_o1(true), wrtAllx1_o2(true), wrtAllx2_o2(true); + bool outAlly_o0(true), outAlly_o1(true), outAlly_o2(true); + if(!wrt_provided) { + if(inDir_provided) { + wrtAllx_o1 = false; // not checked anyway + wrtAllx1_o2 = false; // not checked anyway + res_dimx_o1 = 1; + res_dimx1_o2 = 1; + } + } else { // wrt_provided + wrtAllx_o1 = false; + wrtAllx1_o2 = false; + wrtAllx2_o2 = false; + res_dimx2_o2 = length_wrt; + if(inDir_provided) { // inDir_provided + res_dimx_o1 = 1; + res_dimx1_o2 = 1; } else { - // std::cout<<"wrtAll = "< x1(n, 0); // vector specifying first derivatives. - // first specify coeffs for first dim - // of s across all directions r, then - // second dim, ... - x1[dx1_ind] = 1; -#ifdef _TIME_AD_GENERAL - derivs_run_tape_timer_start(); -#endif - // std::cout<<"Forward 1 x1: "; - // for(int ijk = 0; ijk < x1.size(); ++ijk) - // std::cout< forwardOut; - // forwardOut = ADtape->Forward(1, x1); - ADtape->Forward(1, x1); - // std::cout<<"forwardOut 1 result (dx1_ind = "<< dx1_ind << ", forwardOut.size() = "<< forwardOut.size() <<"): "; - // for(int ijk = 0; ijk < forwardOut.size(); ++ijk) - // std::cout<Reverse(2, w); - // std::cout<<"reverse 2 result: "; - // for(int ijk = 0; ijk < cppad_derivOut.size(); ++ijk) - // std::cout<hessian[wrt_n * wrt_n * dy_ind + wrt_n * vec_ind + vec_ind2] = - cppad_derivOut[dx2_ind * 2 + 1]; - } - else{ - ansList->hessian[wrt_n * wrt_n * dy_ind + wrt_n * vec_ind + vec_ind2] = - CppAD::numeric_limits::quiet_NaN(); - } - } - } + res_dimx_o1 = length_wrt; + res_dimx1_o2 = length_wrt; + } + } + if(outInds_provided) { + outAlly_o0 = false; + outAlly_o1 = false; + outAlly_o2 = false; + res_dimy_o0 = length_outInds; + res_dimy_o1 = length_outInds; + res_dimy_o2 = length_outInds; + } + if(outDir_provided) { + if(maxOrder == 2) { + outAlly_o2 = false; + res_dimy_o2 = 1; + } + if(maxOrder == 1) { + outAlly_o1 = false; + res_dimy_o1 = 1; } + } + + vector value_ans; + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_start(); + #endif + + int strategy_zero = 1; // 1 means do it, 0 means skip. + int strategy_one = 0; // 0=skip. 1=forward regular. 2=reverse subgraph_rev_jac. 3=reverse normal + int strategy_two= 0; // 0=skip. 1=reverse regular. + if(maxOrder > 0) { if (ordersFound[1]) { - BASE *LHS = ansList->jacobian.getPtr() + dy_ind; - if(!infIndicators[dy_ind]){ - if(wrtAll) { - for (size_t vec_ind3 = 0; vec_ind3 < wrt_n; ++vec_ind3, LHS += q) { - *LHS = cppad_derivOut[vec_ind3 * maxOrder]; - } - } else { - double const *wrtVector_p = wrtVector.getConstPtr(); - double const *wrtVector_p_end = wrtVector_p + wrt_n; - for(; wrtVector_p != wrtVector_p_end; LHS += q ) { - *LHS = cppad_derivOut[(static_cast(*wrtVector_p++) - 1) * maxOrder]; - } - } - } else { - for (size_t vec_ind = 0; vec_ind < wrt_n; vec_ind++) { - *LHS = CppAD::numeric_limits::quiet_NaN(); - LHS += q; - } - } - - // for (size_t vec_ind = 0; vec_ind < wrt_n; vec_ind++) { - // if(!infIndicators[dy_ind]){ - // int dx1_ind = wrtVector[vec_ind] - 1; - // ansList->jacobian[vec_ind * q + dy_ind] = - // cppad_derivOut[dx1_ind * maxOrder + 0]; - // } - // else{ - // ansList->jacobian[vec_ind * q + dy_ind] = - // CppAD::numeric_limits::quiet_NaN(); - // } - // } - + // At least Jacobian is needed. + if(!ordersFound[2]) { + // Hessian is not needed. + if((m <= n || outDir_provided) && !inDir_provided) { // possibly we should compare to wrt_n, but we have use cases where comparing to n works better + // It's hard to really know, because reverse offers the subgraph_rev_jac option. + // Use reverse mode b/c we have fewer (or equal) number of outputs as inputs, or because outDir is provided. + if(m > 1 && !outDir_provided) { + // If there are multiple outputs, use subgraph_rev_jac + strategy_one = 2; // use subgraph_rev_jac if jacobian is the final order + if(!ordersFound[0]) { + // If value is not needed, we can skip Forward 0 because subgraph_rev_jac does it itself. + strategy_zero = 0; + } + } else { + // m == 1, or outDir_provided, so we can use regular reverse mode. + strategy_one = 3;// simple reverse + } + } else { + // Use forward mode b/c we have more outputs than inputs, or inDir was provided + strategy_one = 1; // use forward + } + } else { + // Use Forward mode because we need to continue to Hessian. + // In this case, strategy_one is not actually checked below, but we set it here for consistency. + strategy_one = 1; + } + } // done with ordersFound[1]==true + if(ordersFound[2]) // Hessian is needed. + strategy_two = 1; // not really used below, but set up for consistency. + } + // std::cout<<"n: "<Forward(0, independentVars); + // std::cout<<"value_ans.size() = "<value.setSize(res_dimy_o0, false, false); + if(outAlly_o0) { + std::copy(value_ans.begin(), value_ans.end(), ansList->value.getPtr()); + } else { + BASE *LHS = ansList->value.getPtr(); + for(size_t iii=0;iii 0){ + // std::cout<<"entering maxOrder> 0\n"; + vector infIndicators(m, false); // default values will be false + for(size_t inf_ind = 0; inf_ind < m; inf_ind++){ + if(check_inf_nan_gdi(value_ans[inf_ind])) { + infIndicators[inf_ind] = true; + } + } + //std::cout<<"about to set jacobian and hessian sizes\n"; + if (ordersFound[1]) { + ansList->jacobian.setSize(res_dimy_o1, res_dimx_o1, false, true); // true for fill_zero + } + //std::cout<<"done setting jacobian size\n"; + if (ordersFound[2]) { + ansList->hessian.setSize(res_dimx1_o2, res_dimx2_o2, res_dimy_o2, false, true); + } + //std::cout<<"done setting hessian size\n"; + vector cppad_derivOut; + std::vector w(m, 0); + + // begin replacement + if (maxOrder == 1) { + if(strategy_one == 2) { // subgraph_rev_jac + // Reverse mode + // std::cout<<"subgraph_rev_jac version of jacobian\n"; + // std::cout<<"select_domain: "; + // for (size_t i = 0; i < select_domain.size(); i++) { + // std::cout< select_domain(n, false); + for (size_t i_wrt = 0; i_wrt < res_dimx_o1; i_wrt++) { + int x_ind = wrtAllx_o1 ? i_wrt : wrtVector[i_wrt] - 1; + select_domain[x_ind] = true; + } + std::vector select_range(m, false); + for (size_t i_out = 0; i_out < res_dimy_o1; i_out++) { + int y_ind = outAlly_o1 ? i_out : outInds[i_out] - 1; + select_range[y_ind] = true; + } + + std::vector col_2_wrtVecm1; + if(!wrtAllx_o1) { + col_2_wrtVecm1.resize(n, -1); + for(size_t i_wrt = 0; i_wrt < length_wrt; i_wrt++) { + col_2_wrtVecm1[wrtVector[i_wrt] - 1] = i_wrt; + } + } + std::vector row_2_outIndm1; + if(!outAlly_o1) { + row_2_outIndm1.resize(m, -1); + for(size_t i_out = 0; i_out < length_outInds; i_out++) { + row_2_outIndm1[outInds[i_out] - 1] = i_out; + } + } + // std::cout<<"done setting select_domain and select_range\n"; + CppAD::sparse_rcv , std::vector > matrix_out; + ADtape->subgraph_jac_rev(select_domain, select_range, + independentVars, matrix_out); + // std::cout<<"done with subgraph_jac_rev\n"; + size_t nnz = matrix_out.nnz(); + std::vector row_major = matrix_out.row_major(); + size_t this_row, this_col, this_ind; + BASE *LHS = ansList->jacobian.getPtr(); + // size_t max_row(0), max_col(0); + // bool give_output = false; + for(size_t k = 0; k < nnz; k++) + { + this_ind = row_major[k]; + this_row = matrix_out.row()[this_ind]; + this_col = matrix_out.col()[this_ind]; + int x1_ind = wrtAllx_o1 ? this_col : col_2_wrtVecm1[this_col]; + if(x1_ind < 0) { + std::cout<<"Error: x1_ind is negative. This should not happen."< max_row) { + // max_row = this_row; + // give_output = true; + // } + // if(this_col > max_col) { + // max_col = this_col; + // give_output = true; + // } + // if(give_output) { + // give_output = false; + // // std::cout<<"this_ind: "<Reverse(1, w); + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_stop(); + #endif + } else { + cppad_derivOut.resize(maxOrder*n, 0.0); + } + } else // end use_out + { + y_ind = 0; // outDir is used, so we only need one output. + for(size_t iy = 0; iy < m; iy++) w[iy] = outDir[iy]; + cppad_derivOut = ADtape->Reverse(1, w); + } // end use outDir. + if (ordersFound[1]) { // will always be true if maxOrder == 1 + BASE *LHS = ansList->jacobian.getPtr() + i_out; + if(!infIndicators[y_ind]){ + for (size_t i_wrt = 0; i_wrt < res_dimx_o1; ++i_wrt, LHS += res_dimy_o1) { + int x_ind = wrtAllx_o1 ? i_wrt : wrtVector[i_wrt] - 1; + *LHS = cppad_derivOut[y_ind + x_ind * res_dimy_o1]; + } + } else { + for (size_t i_wrt = 0; i_wrt < res_dimx_o1; i_wrt++) { + *LHS = CppAD::numeric_limits::quiet_NaN(); + LHS += outAlly_o1; + } + } + } + if(!outDir_provided) w[y_ind] = 0; + } // end out_ind loop + // end strategy_one == 3 + } else if(strategy_one == 1) + { // m > n + // Forward mode + //std::cout<<"new version of jacobian\n"; + std::vector x1(n, 0); + for (size_t i_wrt = 0; i_wrt < res_dimx_o1; i_wrt++) { + int x1_ind(0); + if(!inDir_provided) { + x1_ind = wrtAllx_o1 ? i_wrt : wrtVector[i_wrt] - 1; + x1[x1_ind] = 1; + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_start(); + #endif + cppad_derivOut = ADtape->Forward(1, x1); + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_stop(); + #endif + } else { // end use_wrt + x1_ind = 0; // inDir is used, so we only need one input. + for(size_t ix = 0; ix < n; ix++) x1[ix] = inDir[ix]; + cppad_derivOut = ADtape->Forward(1, x1); + } // end use inDir. + if (ordersFound[1]) { // will always be true if maxOrder == 1 + BASE *LHS = ansList->jacobian.getPtr() + res_dimy_o1*i_wrt; + + for(size_t i_out = 0; i_out < res_dimy_o1; ++i_out, ++LHS) { + int y_ind = outAlly_o1 ? i_out : outInds[i_out] - 1; + *LHS = infIndicators[y_ind] ? + CppAD::numeric_limits::quiet_NaN() : + cppad_derivOut[y_ind]; + } + } + if(!inDir_provided) x1[x1_ind] = 0; + }// end vec_ind loop + // end strategy_one == 1 + } else { + std::cout<<"Error: in getDerivs: strategy_one is not set correctly for maxOrder == 1."< 1\n"; + // maxOrder > 1: outer loop over vec_ind, inner loop over dy_ind + // strategy_one and strategy_two are actually ignored here for now, b/c we wouldn't be here if not needed. + std::vector cppad_derivOut_F1; + std::vector x1(n, 0); + for (size_t i1_wrt = 0; i1_wrt < res_dimx1_o2; i1_wrt++) { + int x1_ind(0); + if(!inDir_provided) { + x1_ind = wrtAllx1_o2 ? i1_wrt : wrtVector[i1_wrt] - 1; + x1[x1_ind] = 1; + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_start(); + #endif + cppad_derivOut_F1 = ADtape->Forward(1, x1); + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_stop(); + #endif + } else { // end use_wrt + x1_ind = 0; // inDir is used, so we only need one input. + for(size_t ix = 0; ix < n; ix++) x1[ix] = inDir[ix]; + cppad_derivOut_F1 = ADtape->Forward(1, x1); + } // end use inDir. + // std::cout<<"contents of cppad_derivOut_F1: "; + // for(size_t ijk = 0; ijk < cppad_derivOut_F1.size(); ++ijk) + // std::cout<Reverse(2, w); + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_stop(); + #endif + } else { + cppad_derivOut.resize(maxOrder * n, 0.0); + } + } else { // outDir_provided is true + y_ind = 0; // outDir is used, so we only need one output. + for(size_t iy = 0; iy < m; iy++) w[iy] = outDir[iy]; + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_start(); + #endif + cppad_derivOut = ADtape->Reverse(2, w); + #ifdef _TIME_AD_GENERAL + derivs_run_tape_timer_stop(); + #endif + } // end outDir_provided + // record Hessian outputs + if(!infIndicators[y_ind]){ + for (size_t i2_wrt = 0; i2_wrt < res_dimx2_o2; i2_wrt++) { + int x2_ind(0); + x2_ind = wrtAllx2_o2 ? i2_wrt : wrtVector[i2_wrt] - 1; + + ansList->hessian[res_dimx1_o2 * res_dimx2_o2 * i_out + res_dimx2_o2 * i1_wrt + i2_wrt] = + cppad_derivOut[x2_ind * 2 + 1]; + } + } else { + for (size_t i2_wrt = 0; i2_wrt < res_dimx2_o2; i2_wrt++) { + ansList->hessian[res_dimx1_o2 * res_dimx2_o2 * i_out + res_dimx2_o2 * i1_wrt + i2_wrt] = + CppAD::numeric_limits::quiet_NaN(); + } + } + // end use_out2 + if(!outDir_provided) w[y_ind] = 0; // outDir is used, so we only need one output. + } // end out_ind loop + + if (ordersFound[1]) { + BASE *LHS = ansList->jacobian.getPtr() + res_dimy_o1 * i1_wrt; + + for(size_t i_out = 0; i_out < res_dimy_o1; ++i_out, ++LHS) { + int y_ind = outAlly_o1 ? i_out : outInds[i_out] - 1; + *LHS = infIndicators[y_ind] ? + CppAD::numeric_limits::quiet_NaN() : + cppad_derivOut_F1[y_ind]; + } + } + if(!inDir_provided) x1[x1_ind] = 0; + } + } + } // end else + + // end replacement + #ifdef _TIME_AD_GENERAL + derivs_getDerivs_timer_stop(); + #endif } -#ifdef _TIME_AD_GENERAL - derivs_getDerivs_timer_stop(); -#endif -}; void nimbleFunctionCppADbase::getDerivs_meta(nimbleCppADinfoClass &ADinfo, const NimArr<1, double> &derivOrders, const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, CppAD::AD > &inDir, + const NimArr<1, CppAD::AD > &outDir, const nimbleCppADrecordingInfoClass &nimRecInfo, nimSmartPtr &ansList) { // std::cout<<"Entering getDerivs_meta"< &derivOrders, const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, double> &inDir, + const NimArr<1, double> &outDir, nimSmartPtr &ansList) { // std::cout<<"Entering getDerivs"< outInds, inDir, outDir; // being length 0 will make them not used. + NimArr<1, CppAD::AD > meta_inDir, meta_outDir; if(ADinfo.ADtape_empty() || reset) { // Delete previous tape if it exists. if(!ADinfo.ADtape_empty()) @@ -515,6 +878,9 @@ void nimbleFunctionCppADbase::getDerivs_calculate_internal(nimbleCppADinfoClass &innerTape, derivOrders_meta, wrtVector, + outInds, + meta_inDir, + meta_outDir, ansList_meta); for(int iii = 0; iii < length_wrt; ++iii) dependentVars[iii] = ansList_meta->jacobian[iii]; @@ -596,6 +962,9 @@ void nimbleFunctionCppADbase::getDerivs_calculate_internal(nimbleCppADinfoClass ADinfo.ADtape(), derivOrders_nested, wrtVector, // NOTE: This will not behave fully correctly in non-default use without further thought. + outInds, + inDir, + outDir, ansList_nested); // std::cout<<"getDerivs_calculate_internal E"< &derivOrders, const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, double> &inDir, + const NimArr<1, double> &outDir, nimSmartPtr &ansList); nimSmartPtr getDerivs_wrapper(nimbleCppADinfoClass &ADinfo, const NimArr<1, double> &derivOrders, - const NimArr<1, double> &wrtVector){ + const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, double> &inDir, + const NimArr<1, double> &outDir){ nimSmartPtr ansList = new NIMBLE_ADCLASS; - getDerivs(ADinfo, derivOrders, wrtVector, ansList); + getDerivs(ADinfo, derivOrders, wrtVector, outInds, inDir, outDir, ansList); return(ansList); } void getDerivs_meta(nimbleCppADinfoClass &ADinfo, const NimArr<1, double> &derivOrders, const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, CppAD::AD > &inDir, + const NimArr<1, CppAD::AD > &outDir, const nimbleCppADrecordingInfoClass &nimRecInfo, nimSmartPtr &ansList); nimSmartPtr getDerivs_wrapper_meta(nimbleCppADinfoClass &ADinfo, const NimArr<1, double> &derivOrders, const NimArr<1, double> &wrtVector, + const NimArr<1, double> &outInds, + const NimArr<1, CppAD::AD > &inDir, + const NimArr<1, CppAD::AD > &outDir, const nimbleCppADrecordingInfoClass &nimRecInfo){ nimSmartPtr ansList = new NIMBLE_ADCLASS_META; - getDerivs_meta(ADinfo, derivOrders, wrtVector, nimRecInfo, ansList); + getDerivs_meta(ADinfo, derivOrders, wrtVector, outInds, inDir, outDir, nimRecInfo, ansList); return(ansList); } diff --git a/packages/nimble/src/nimble.cpp b/packages/nimble/src/nimble.cpp index 22d1a4bff..d27e1c065 100644 --- a/packages/nimble/src/nimble.cpp +++ b/packages/nimble/src/nimble.cpp @@ -39,6 +39,9 @@ R_CallMethodDef CallEntries[] = { FUN(C_rmvt_chol, 4), FUN(C_dlkj_corr_cholesky, 4), FUN(C_rlkj_corr_cholesky, 2), + FUN(C_PDinverse_logdet, 2), + FUN(C_dmnorm_inv_ld, 4), + FUN(C_rmnorm_inv_ld, 2), FUN(C_dexp_nimble, 3), FUN(C_rexp_nimble, 2), FUN(C_pexp_nimble, 4), diff --git a/packages/nimble/tests/testthat/AD_test_utils.R b/packages/nimble/tests/testthat/AD_test_utils.R index adce75f6f..61e55045a 100644 --- a/packages/nimble/tests/testthat/AD_test_utils.R +++ b/packages/nimble/tests/testthat/AD_test_utils.R @@ -253,7 +253,6 @@ test_AD2_oneCall <- function(Robj, Cobj, wrt_all <- 1:len_record else wrt_all <- wrt - do_one_set("run", name = "run", doR = TRUE, doC = TRUE, tapingLevel = 0, fixedOrder = 0) do_one_set("derivsRun", order = order, metaLevel = 0, name = "derivsRun", @@ -301,7 +300,7 @@ test_AD2_oneCall <- function(Robj, Cobj, abs_threshold = RRabsThresh, info = paste0("(RR order ", o,")")) pass <- pass && localPass - if(verbose && !pass) { + if(verbose && !localPass) { cat(paste('Some R-to-R derivatives do not match for order', o, '.\n')) } } @@ -310,7 +309,7 @@ test_AD2_oneCall <- function(Robj, Cobj, abs_threshold = RCabsThresh, info = paste0("(RC order ", o,")")) pass <- pass && localPass - if(verbose && !pass) { + if(verbose && !localPass) { cat(paste('Some C-to-R derivatives do not match for order', o, '.\n')) print(RansSet[[1]]) print(CansSet) @@ -321,7 +320,7 @@ test_AD2_oneCall <- function(Robj, Cobj, abs_threshold = CCabsThresh, info = paste0("(CC order ", o, ")")) pass <- pass && localPass - if(verbose && !pass) { + if(verbose && !localPass) { cat(paste('Some C-to-C derivatives do not match for order', o, '.\n')) print(CansSet[[1]]) print(CansSet[-1]) @@ -1788,7 +1787,7 @@ calcNodesForDerivs <- nimbleFunction( test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = NULL, calcNodes = NULL, wrt = NULL, newUpdateNodes = NULL, newConstantNodes = NULL, relTol = c(1e-15, 1e-8, 1e-3, 1e-3, 1e-14), absTolThreshold = 0, useFasterRderivs = FALSE, useParamTransform = FALSE, - checkDoubleTape = TRUE, checkCompiledValuesIdentical = TRUE, checkDoubleUncHessian = TRUE, + checkDoubleTape = TRUE, checkCompiledValuesIdentical = TRUE, check01vs012jacIdentical = TRUE, checkDoubleUncHessian = TRUE, doAllUncHessian = TRUE, seed = 1, verbose = FALSE, debug = FALSE){ if(!is.null(seed)) set.seed(seed) @@ -1830,6 +1829,7 @@ test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = N useFasterRderivs = useFasterRderivs, useParamTransform = useParamTransform, checkDoubleTape = checkDoubleTape, checkCompiledValuesIdentical = checkCompiledValuesIdentical, + check01vs012jacIdentical = check01vs012jacIdentical, checkDoubleUncHessian = checkDoubleUncHessian, doAllUncHessian = doAllUncHessian, verbose = verbose, debug = debug)) ## max. lik. use case @@ -1866,6 +1866,7 @@ test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = N useFasterRderivs = useFasterRderivs, useParamTransform = useParamTransform, checkDoubleTape = checkDoubleTape, checkCompiledValuesIdentical = checkCompiledValuesIdentical, + check01vs012jacIdentical = check01vs012jacIdentical, checkDoubleUncHessian = checkDoubleUncHessian, doAllUncHessian = doAllUncHessian, verbose = verbose, debug = debug)) @@ -1910,6 +1911,7 @@ test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = N useFasterRderivs = useFasterRderivs, useParamTransform = useParamTransform, checkDoubleTape = checkDoubleTape, checkCompiledValuesIdentical = checkCompiledValuesIdentical, + check01vs012jacIdentical = check01vs012jacIdentical, checkDoubleUncHessian = checkDoubleUncHessian, doAllUncHessian = doAllUncHessian, verbose = verbose, debug = debug)) @@ -1957,6 +1959,7 @@ test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = N useFasterRderivs = useFasterRderivs, useParamTransform = useParamTransform, checkDoubleTape = checkDoubleTape, checkCompiledValuesIdentical = checkCompiledValuesIdentical, + check01vs012jacIdentical = check01vs012jacIdentical, checkDoubleUncHessian = checkDoubleUncHessian, doAllUncHessian = doAllUncHessian, verbose = verbose, debug = debug)) @@ -1993,6 +1996,7 @@ test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = N useFasterRderivs = useFasterRderivs, useParamTransform = useParamTransform, checkDoubleTape = checkDoubleTape, checkCompiledValuesIdentical = checkCompiledValuesIdentical, + check01vs012jacIdentical = check01vs012jacIdentical, checkDoubleUncHessian = checkDoubleUncHessian, doAllUncHessian = doAllUncHessian, verbose = verbose, debug = debug)) } else { @@ -2023,6 +2027,7 @@ test_ADModelCalculate <- function(model, name = 'unknown', x = 'given', xNew = N useFasterRderivs = useFasterRderivs, useParamTransform = useParamTransform, checkDoubleTape = checkDoubleTape, checkCompiledValuesIdentical = checkCompiledValuesIdentical, + check01vs012jacIdentical = check01vs012jacIdentical, checkDoubleUncHessian = checkDoubleUncHessian, doAllUncHessian = doAllUncHessian, verbose = verbose, debug = debug)) } @@ -2036,7 +2041,7 @@ test_ADModelCalculate_internal <- function(model, name = 'unknown', xOrig = NULL newUpdateNodes = NULL, newConstantNodes = NULL, relTol = c(1e-15, 1e-8, 1e-3, 1e-3, 1e-14), absTolThreshold = 0, useFasterRderivs = FALSE, useParamTransform = FALSE, checkDoubleTape = TRUE, - checkCompiledValuesIdentical = TRUE, checkDoubleUncHessian = TRUE, + checkCompiledValuesIdentical = TRUE, check01vs012jacIdentical = TRUE, checkDoubleUncHessian = TRUE, doAllUncHessian = TRUE, verbose = FALSE, debug = FALSE){ @@ -2045,7 +2050,6 @@ test_ADModelCalculate_internal <- function(model, name = 'unknown', xOrig = NULL on.exit(local_edition(saved_edition)) test_that(paste0("Derivatives of calculate for model ", name), { - if(exists('paciorek') && paciorek == 0) browser() if(is.null(calcNodes)) calcNodes <- model$getNodeNames() @@ -2144,7 +2148,6 @@ test_ADModelCalculate_internal <- function(model, name = 'unknown', xOrig = NULL for(case in 1:2) { for(idx in seq_along(xList)) { - if(exists('paciorek') && paciorek == idx) browser() if(verbose) { if(case == 1) { if(idx == 1) { @@ -2544,8 +2547,11 @@ test_ADModelCalculate_internal <- function(model, name = 'unknown', xOrig = NULL } ## explicit comparison of first derivs; - ## both of these are reverse mode because 2nd order reverse also invokes first order reverse - expect_identical(cOutput01$jacobian, cOutput012$jacobian) + ## Originally, both of these are reverse mode because 2nd order reverse also invokes first order reverse. + ## As of version 1.4.0, these are not identical for at least test-ADdmnorm test. + if(check01vs012jacIdentical) { + expect_identical(cOutput01$jacobian, cOutput012$jacobian) + } else expect_equal(cOutput01$jacobian, cOutput012$jacobian) expect_identical(cOutput12$jacobian, cOutput012$jacobian) diff --git a/packages/nimble/tests/testthat/test-ADdirsAndInds.R b/packages/nimble/tests/testthat/test-ADdirsAndInds.R new file mode 100644 index 000000000..21503f69f --- /dev/null +++ b/packages/nimble/tests/testthat/test-ADdirsAndInds.R @@ -0,0 +1,720 @@ +# Tests of controlling nimDerivs behavior with +# combinations of wrt, outInds, inDir, and/or outDir. +source(system.file(file.path('tests', 'testthat', 'AD_test_utils.R'), package = 'nimble')) +EDopt <- nimbleOptions("enableDerivs") +BMDopt <- nimbleOptions("buildModelDerivs") +nimbleOptions(enableDerivs = TRUE) +nimbleOptions(buildModelDerivs = TRUE) + +test_that('Derivatives with double-taping work', + { + ADfun1 <- nimbleFunction( + setup = function(){}, + run = function(){}, + methods = list( + jacobianRun = function(x = double(1)) { + out <- derivs(zo(x), wrt = 1:2, order = 1) + ans <- nimNumeric(length = 6, value = out$jacobian) + returnType(double(1)) + return(ans) + }, + zo = function(x = double(1)) { + ans <- c(x[1]^3, x[2]^3, (x[1]^2*x[2]^2)) + return(ans) + returnType(double(1)) + }, + metaDerivsRun = function(x = double(1), + order = double(1)) { + out <- derivs(jacobianRun(x), wrt = 1:2, + order = order) + returnType(ADNimbleList()) + return(out) + } + ), buildDerivs = c('jacobianRun', 'zo') + ) + + ADfunInst <- ADfun1() + xRec <- c(2.2, 3.3) + x <- c(1.6, 2.8) + Rderivs <- ADfunInst$jacobianRun(x) + Rderivs <- ADfunInst$metaDerivsRun(x, order = 1) + temporarilyAssignInGlobalEnv(ADfunInst) + cADfunInst <- compileNimble(ADfunInst) + cADfunInst$metaDerivsRun(xRec, 1) + cderivs <- cADfunInst$metaDerivsRun(x, 1) + expect_equal(cderivs$value, Rderivs$value) + expect_equal(cderivs$jacobian, Rderivs$jacobian, tolerance = 1e-5) + expect_equal(cderivs$hessian, Rderivs$hessian) + } +) + + +test_that('Derivatives with double-taping, indices, and directions work', + { + ADfun1 <- nimbleFunction( + setup = function(){}, + run = function(x = double(1)) { + ans <- c(x[1]^3, x[2]^3, (x[1]^2*x[2]^2)) + return(ans) + returnType(double(1)) + }, + methods = list( + derivsRun = function(x = double(1), + wrt = double(1), + outInds = double(1), + inDir = double(1), + outDir = double(1), + order = double(1)) { + ans <- derivs(run(x), wrt = wrt, outInds = outInds, + inDir = inDir, outDir = outDir, order = order) + return(ans) + returnType(ADNimbleList()) + }, + jacobianRun = function(x = double(1)) { + out <- derivs(run(x), wrt = 1:2, order = 1) + ans <- nimNumeric(length = 6, value = out$jacobian) + returnType(double(1)) + return(ans) + }, + + metaDerivsRun = function(x = double(1), + order = double(1)) { + out <- derivs(jacobianRun(x), wrt = 1:2, + order = order) + returnType(ADNimbleList()) + return(out) + } + ), buildDerivs = c('jacobianRun', 'run') + ) + + ADfunInst <- ADfun1() + xRec <- c(2.2, 3.3) + x <- c(1.6, 2.8) + Rderivs <- ADfunInst$jacobianRun(x) + Rderivs <- ADfunInst$metaDerivsRun(x, order = 1) + temporarilyAssignInGlobalEnv(ADfunInst) + cADfunInst <- compileNimble(ADfunInst) + + ## record + cADfunInst$derivsRun(x = xRec, + wrt = 1:2, + outInds = numeric(), + inDir = numeric(), + outDir = numeric(), + order = 1:2) + + ## Future args: + wrt <- list(1:2, 1, 2, c(2, 1), c(2, 1, 2)) + outInds <- list(1:3, 1, 2, 3, c(3, 1), c(3, 1, 3)) + inDirs <- list(c(-0.3, -0.7)) + outDirs <- list(c(-2, -1, -3)) + orders <- list(0, 1, 2, c(0, 1), c(0, 2), c(1, 2), 0:2) + + correct <- cADfunInst$derivsRun(x = x, + wrt = numeric(), + outInds = numeric(), + inDir = numeric(), + outDir = numeric(), + order = 0:2) + + # Here is a comprehensive testing function + + check_over_orders <- function(x, orders, wrt=numeric(), + outInds=numeric(), inDir=numeric(), + outDir=numeric()) { + n <- length(x) + for(ord in orders) { + this_x <- x + rnorm(length(x), 0, 0.1) + check <- cADfunInst$derivsRun(x = this_x, + wrt = wrt, + outInds = outInds, + inDir = inDir, + outDir = outDir, + order = ord) + correct <- cADfunInst$derivsRun(x = this_x, + wrt = numeric(), + outInds = numeric(), + inDir = numeric(), + outDir = numeric(), + order = 0:2) + m <- nrow(correct$jacobian) + if(0 %in% ord) { + correct_value <- correct$value + if(length(outInds)) + correct_value <- correct_value[outInds,drop=FALSE] + expect_equal(check$value, correct_value) + } + if(1 %in% ord) { + if(length(inDir) > 0) { + correct_jac <- correct$jacobian %*% matrix(inDir, ncol=1) + if(length(outInds) > 0) + correct_jac <- correct_jac[outInds,, drop=FALSE] + } else { + if(length(outDir) > 0 && !(2 %in% ord)) + correct_jac <- matrix(outDir, nrow=1) %*% correct$jacobian + else { + correct_jac <- correct$jacobian + if(length(outInds) > 0) + correct_jac <- correct_jac[outInds,, drop=FALSE] + } + if(length(wrt) > 0) + correct_jac <- correct_jac[,wrt,drop=FALSE] + } + expect_equal(check$jacobian, + correct_jac) + } + if(2 %in% ord) { + if(length(inDir) && length(outDir)) { + correct_hess <- inDir %*% + apply(correct$hessian, 1, \(H) H %*% outDir) |> + array(dim=c(1, n, 1)) + if(length(wrt)) + correct_hess <- correct_hess[,wrt,, drop=FALSE] + } else if(length(inDir) && !length(outDir)) { + correct_hess <- correct$hessian |> + apply(1, \(x) inDir %*% x) |> t() |> array(dim=c(1, n, m)) + if(length(wrt)) + correct_hess <- correct_hess[,wrt,, drop=FALSE] + if(length(outInds)) + correct_hess <- correct_hess[,, outInds, drop=FALSE] + } else if(!length(inDir) && length(outDir)) { + correct_hess <- apply(correct$hessian, 1, \(H) H %*% outDir) |> + array(dim=c(n,n,1)) + if(length(wrt)) + correct_hess <- correct_hess[wrt, wrt, , drop=FALSE] + } else { + correct_hess <- correct$hessian + if(length(wrt)) + correct_hess <- correct_hess[wrt, wrt, , drop=FALSE] + if(length(outInds)) + correct_hess <- correct_hess[,, outInds, drop=FALSE] + } + expect_equal(check$hessian, + correct_hess) + } + } + } + + ## Next we use the above with all combinations of inputs + ## This is rather verbose. If we included a case "numeric()" + ## in each of the input lists above, then we could simply use + ## one big quad-nested loop. + ## However, in building up, it was useful to work in + ## systematic combinations as done below. + ## We could consolidate this in the future. + ## But they run fast, as they use a single compiled object. + + ## wrt alone + for(this_wrt in wrt) { + check_over_orders(x=x,orders,wrt=this_wrt) + } + ## outInds alone + for(this_outInds in outInds) { + check_over_orders(x = x, orders, outInds = this_outInds) + } + ## wrt and outInds together -- the two kinds of indexing + for(this_wrt in wrt) { + for(this_outInds in outInds) { + check_over_orders(x=x, orders=orders, wrt=this_wrt, outInds=this_outInds) + } + } + ## inDir alone + for(this_inDir in inDirs) { + check_over_orders(x=x, orders=orders, inDir = this_inDir) + } + ## inDir and wrt -- the two kinds of "x" controls + for(this_wrt in wrt) { + for(this_inDir in inDirs) { + check_over_orders(x=x, orders=orders, wrt = this_wrt, inDir = this_inDir) + } + } + ## outDir alone + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, outDir = this_outDir) + } + ## outDir and outInds -- the two kinds of "y" controls + for(this_outInds in outInds) { + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, outInds = this_outInds, outDir = this_outDir) + } + } + ## inDir and outDir -- the two kinds of directions + for(this_inDir in inDirs) { + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, outDir = this_outDir, inDir=this_inDir) + } + } + + ## wrt, inDir, and outDir + for(this_wrt in wrt) { + for(this_inDir in inDirs) { + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, wrt=this_wrt, + outDir = this_outDir, inDir=this_inDir) + } + } + } + + ## wrt and outDir + for(this_wrt in wrt) { + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, wrt=this_wrt, outDir = this_outDir) + } + } + + ## outInds, inDir, and outDir + for(this_outInds in outInds) { + for(this_inDir in inDirs) { + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, + outInds = this_outInds, outDir = this_outDir, inDir=this_inDir) + } + } + } + + ## outInds and inDir + for(this_outInds in outInds) { + for(this_inDir in inDirs) { + check_over_orders(x=x, orders=orders, + outInds = this_outInds, inDir=this_inDir) + } + } + + ## all 4 controls + for(this_wrt in wrt) { + for(this_outInds in outInds) { + for(this_inDir in inDirs) { + for(this_outDir in outDirs) { + check_over_orders(x=x, orders=orders, wrt=this_wrt, + outInds = this_outInds, outDir = this_outDir, inDir=this_inDir) + } + } + } + } + + + ## ## What follows is some initial testing for manually inspecting cases. + ## ## It may be useful if anything breaks to be able to run through these. + ## foo <- function(t, x_, inDir) cADfunInst$run(x_ + t*inDir) + ## ###################### + ## ## play with order 0 # + ## ###################### + ## ## both wrt ## + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + ## ## wrt flipped, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 1, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 2, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 3, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + + ## ## all wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## ## single wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## # permuted wrt and permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 0) + + ## ## inDir, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = numeric(), + ## outInds = numeric(), + ## inDir = c(-.3, -.7), + ## outDir = numeric(), + ## order = 0) + + ## ## outDir, all in + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = c(-2, -1, -3), + ## order = 0) + + + ## ##################### + ## ## play with order 1 + ## ##################### + ## ## both wrt + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + ## ## wrt flipped, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 1, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 2, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 3, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + + ## ## all wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## ## single wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## # permuted wrt and permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1) + + ## ## inDir, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = numeric(), + ## outInds = numeric(), + ## inDir = c(-.3, -.7), + ## outDir = numeric(), + ## order = 1) + + ## ## outDir, all in + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = c(-2, -1, -3), + ## order = 1) + + + ## ##################### + ## ## play with order 2 + ## ##################### + ## ## both wrt + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + ## ## wrt flipped, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 1, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 2, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 3, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + + ## ## all wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## single wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## # permuted wrt and permuted out + ## cADfunInst$Derivsrun(x = x, + ## wrt = c(2, 1), + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 2) + + ## ## inDir, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = numeric(), + ## outInds = numeric(), + ## inDir = c(-.3, -.7), + ## outDir = numeric(), + ## order = 2) + + ## ## outDir, all in + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = c(-2, -1, -3), + ## order = 2) + + ## # apply(correct$hessian, 1, \(H) H %*% c(-2, -1, -3)) + + ## ## inDir and outDir + ## cADfunInst$derivsRun(x = x, + ## wrt = numeric(), + ## outInds = numeric(), + ## inDir = c(-.3, -.7), + ## outDir = c(-2, -1, -3), + ## order = 2) + + ## #c(-.3, -.7) %*% apply(correct$hessian, 1, \(H) H %*% c(-2, -1, -3)) + + + ## ##################### + ## ## play with orders 1 & 2 + ## ##################### + ## ## both wrt + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + ## ## wrt flipped, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## single wrt, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 1, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 2, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## all wrt, single out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = 3, + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + + ## ## all wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## single wrt, permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = 2, + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## # permuted wrt and permuted out + ## cADfunInst$derivsRun(x = x, + ## wrt = c(2, 1), + ## outInds = c(3, 1), + ## inDir = numeric(), + ## outDir = numeric(), + ## order = 1:2) + + ## ## inDir, all out + ## cADfunInst$derivsRun(x = x, + ## wrt = numeric(), + ## outInds = numeric(), + ## inDir = c(-.3, -.7), + ## outDir = numeric(), + ## order = 1:2) + + ## ## outDir, all in + ## cADfunInst$derivsRun(x = x, + ## wrt = 1:2, + ## outInds = numeric(), + ## inDir = numeric(), + ## outDir = c(-2, -1, -3), + ## order = 1:2) + + ## # apply(correct$hessian, 1, \(H) H %*% c(-2, -1, -3)) + + ## ## inDir and outDir + ## cADfunInst$derivsRun(x = x, + ## wrt = numeric(), + ## outInds = numeric(), + ## inDir = c(-.3, -.7), + ## outDir = c(-2, -1, -3), + ## order = 1:2) + + ## #c(-.3, -.7) %*% apply(correct$hessian, 1, \(H) H %*% c(-2, -1, -3)) + } +) diff --git a/packages/nimble/tests/testthat/test-ADdists.R b/packages/nimble/tests/testthat/test-ADdists.R index 2474bfbbd..6843da113 100644 --- a/packages/nimble/tests/testthat/test-ADdists.R +++ b/packages/nimble/tests/testthat/test-ADdists.R @@ -74,7 +74,6 @@ dirch_test_log <- make_AD_test2( wrt = c('x', 'alpha'), inputs = list(record = list(x = c(log(.2/.5), log(.3/.5)), alpha = c(2, 4, 4), log = 1), test = list(x = c(log(.4/.45), log(.15/.45)), alpha = c(3, 2, 5), log = 0))) - dirch_test_out <- test_AD2(dirch_test_log) ## Dirichlet without log argument @@ -145,7 +144,7 @@ wish_test_log <- make_AD_test2( for(j in 1:4) { for(i in 1:j) { x2D[i,j] <- x2D[j,i] <- x[indOrig] - chol2D[i,j] <- chol[indOrig] + chol2D[i,j] <- chol[indOrig] indOrig <- indOrig + 1 } } @@ -169,7 +168,7 @@ wish_test_out <- test_AD2(wish_test_log) ## inverse Wishart -set.seed(123) +set.seed(123) wRec <- rinvwish_chol(1, cholRec, df = 7) wTest <- rinvwish_chol(1, cholTest, df = 7) @@ -192,7 +191,7 @@ invwish_test_log <- make_AD_test2( for(j in 1:4) { for(i in 1:j) { x2D[i,j] <- x2D[j,i] <- x[indOrig] - chol2D[i,j] <- chol[indOrig] + chol2D[i,j] <- chol[indOrig] indOrig <- indOrig + 1 } } @@ -225,10 +224,10 @@ lkj_test_log <- make_AD_test2( name = "dlkj_corr_cholesky manual", opParam = list(name = "dlkj_corr_cholesky manual"), expr = quote({ - # correlation matrix inverse transform + # correlation matrix inverse transform z <- nimMatrix(nrow = 5, ncol = 5) u <- nimMatrix(nrow = 5, ncol = 5) - + j <- 1L i <- 1L cnt <- 1L diff --git a/packages/nimble/tests/testthat/test-ADdmnorm.R b/packages/nimble/tests/testthat/test-ADdmnorm.R new file mode 100644 index 000000000..b5e33a73b --- /dev/null +++ b/packages/nimble/tests/testthat/test-ADdmnorm.R @@ -0,0 +1,1012 @@ +source(system.file(file.path('tests', 'testthat', 'AD_test_utils.R'), package = 'nimble')) +BMDopt <- nimbleOptions("buildModelDerivs") +nimbleOptions(enableDerivs = TRUE) +nimbleOptions(buildModelDerivs = TRUE) + +nimbleOptions(useADcholAtomic = TRUE) +nimbleOptions(useADsolveAtomic = TRUE) +nimbleOptions(useADmatMultAtomic = TRUE) +nimbleOptions(useADmatInverseAtomic = TRUE) + +RwarnLevel <- options('warn')$warn +options(warn = 1) + +relTol <- eval(formals(test_ADModelCalculate)$relTol) +relTol[3] <- 1e-6 +relTol[4] <- 1e-4 + +# Standalone test +test_that("basics of R calls to PDinverse_logdet and dmnormAD work", { + set.seed(1) + cov <- crossprod(matrix(rnorm(9), nrow=3)) + PDild <- PDinverse_logdet(cov) + prec <- solve(cov) + expect_equal(PDild[1:9][upper.tri(prec, diag=TRUE)], prec[upper.tri(prec, diag=TRUE)]) + expect_equal(PDild[10], determinant(cov, TRUE)$modulus[1]) + mu <- c(.2, .3, .1) + x <- c(.1, .2, .3) + chol_cov <- chol(cov) + expect_equal(dmnorm_inv_ld(x, mu, mat = cov, inv_ld = PDild, prec_param=FALSE), + dmnorm_chol(x, mu, chol_cov, prec_param=FALSE)) + expect_equal(dmnorm_inv_ld(x, mu, mat = cov, inv_ld = PDild, prec_param=FALSE, log=TRUE), + dmnorm_chol(x, mu, chol_cov, prec_param=FALSE, log=TRUE)) + + PDild2 <- PDinverse_logdet(prec) + expect_equal(PDild2[1:9][upper.tri(cov, diag=TRUE)], cov[upper.tri(cov, diag=TRUE)]) + expect_equal(PDild2[10], determinant(prec, TRUE)$modulus[1]) + chol_prec <- chol(prec) + expect_equal(dmnorm_inv_ld(x, mu, mat = prec, inv_ld = PDild2, prec_param=TRUE), + dmnorm_chol(x, mu, chol_prec, prec_param=TRUE)) + expect_equal(dmnorm_inv_ld(x, mu, mat = prec, inv_ld = PDild2, prec_param=TRUE, log=TRUE), + dmnorm_chol(x, mu, chol_prec, prec_param=TRUE, log=TRUE)) +}) + +# Test in a model with dmnormAD + +set.seed(1) +n <- 3 +locs <- cbind(runif(n),runif(n)) +dd <- matrix(nrow = n, ncol = n) +for(i in 1:n) for(j in 1:n) dd[i,j] <- sqrt(sum((locs[i,]-locs[j,])^2)) +code <- nimbleCode({ + sigma ~ dunif(0, 20) + rho ~ dunif(0,3) + cov[1:3,1:3] <- sigma*sigma*exp(-dd[1:3,1:3]/rho) + inv_ld[1:10] <- PDinverse_logdet(cov[1:3, 1:3]) + x[1:3] ~ dmnormAD(mean = mu[1:3], cov = cov[1:3, 1:3]) +}) +inits <- list(sigma = 0.8, rho = 1.2, mu = c(.2, .3, .1)) +constants = list(dd=dd) +data = list(x = c(.1, .2, .3)) +model <- nimbleModel(code, data = data, constants = constants, inits = inits, buildDerivs=TRUE) + +relTolTmp <- relTol +relTolTmp[2] <- 1e-6 +relTolTmp[3] <- 1e-2 +relTolTmp[4] <- 1e-2 +relTolTmp[5] <- 1e-13 + +test_ADModelCalculate(model, useParamTransform = TRUE, relTol = relTolTmp, + newConstantNodes = list(dd = dd*1.1, x=c(.15, .25, .35)), + newUpdateNodes = list(mu = c(.22, .32, .12)), + checkCompiledValuesIdentical = FALSE, check01vs012jacIdentical = FALSE, + checkDoubleUncHessian = TRUE, + useFasterRderivs = TRUE, verbose = FALSE, name = 'dmnormAD with inv_ld') + +## Depending on the order in which the test above appears in the sequence of testing, +## we have seen the following crash: +# ============================================ +# testing HMC/MAP-based scenario +# -------------------------------------------- +# Testing initial wrt values with initial constantNodes +# Using wrt: sigma rho +# corrupted size vs. prev_size while consolidating +# Aborted + +cov <- crossprod(matrix(rnorm(25), 5)) +mu <- c(0.2,0.1,0.3,0.15,0.25) +x <- c(0.1,0.2,0.15,0.3,0.12) + +PDinverse_logdet_test <- make_AD_test2( + op = list( + name = "PDinverse_logdet with prec_param FALSE: positive definite inverse and log determinant test", + opParam = list(name = "PDinverse_logdet"), + expr = quote({ + pdl <- PDinverse_logdet(mat) + out <- pdl + }), + args = list( + mat = quote(double(2)) + ), + outputType = quote(double(1)) + ), + argTypes = c(mat='double(2)'), + wrt = c('mat'), + inputs = list(record = list(mat = cov), + test = list(mat = cov+0.1)) +) +PDinverse_logdet_test_out <- test_AD2(PDinverse_logdet_test) + +dmnormAD_test1 <- make_AD_test2( + op = list( + name = "dmnormAD test with cov input and log fixed", + opParam = list(name = "dmormAD test 1"), + expr = quote({ + pld <- PDinverse_logdet(cov) + logProb <- dmnorm_inv_ld(x, mu, mat=cov, inv_ld=pld, prec_param=0, log=TRUE) + out <- logProb + }), + args = list( + x = quote(double(1)), + mu = quote(double(1)), + cov = quote(double(2)) + ), + outputType = quote(double(0)) + ), + argTypes = c(x='double(1)', mu='double(1)', cov='double(2)'), + wrt = c('x','mu','mat'), + inputs = list(record = list(cov = cov, mu=mu, x=x), + test = list(cov = cov+0.1, mu=mu-0.07, x=x+0.03)) +) +dmnormAD_test1_out <- test_AD2(dmnormAD_test1) + +dmnormAD_test1b <- make_AD_test2( + op = list( + name = "dmnormAD test with cov input and log dynamics", + opParam = list(name = "dmormAD test 1"), + expr = quote({ + pld <- PDinverse_logdet(cov) + logProb <- dmnorm_inv_ld(x, mu, mat=cov, inv_ld=pld, prec_param=0, log=log) + out <- logProb + }), + args = list( + x = quote(double(1)), + mu = quote(double(1)), + cov = quote(double(2)), + log = quote(double(0)) + ), + outputType = quote(double(0)) + ), + argTypes = c(x='double(1)', mu='double(1)', cov='double(2)', log='double(0)'), + wrt = c('x','mu','mat'), + inputs = list(record = list(cov = cov, mu=mu, x=x, log=1), + test = list(cov = cov+0.1, mu=mu-0.07, x=x+0.03, log=0)) +) +dmnormAD_test1b_out <- test_AD2(dmnormAD_test1b) + +prec <- solve(cov) +dmnormAD_test2 <- make_AD_test2( + op = list( + name = "dmnormAD test with prec input and log fixed", + opParam = list(name = "dmormAD test 1"), + expr = quote({ + pld <- PDinverse_logdet(prec) + logProb <- dmnorm_inv_ld(x, mu, mat=prec, inv_ld=pld, prec_param=1, log=TRUE) + out <- logProb + }), + args = list( + x = quote(double(1)), + mu = quote(double(1)), + prec = quote(double(2)) + ), + outputType = quote(double(0)) + ), + argTypes = c(x='double(1)', mu='double(1)', prec='double(2)'), + wrt = c('x','mu','mat'), + inputs = list(record = list(prec = prec, mu=mu, x=x), + test = list(prec = prec+0.1, mu=mu-0.07, x=x+0.03)) +) +dmnormAD_test2_out <- test_AD2(dmnormAD_test2) + +dmnormAD_test2b <- make_AD_test2( + op = list( + name = "dmnormAD test with prec input and log dynamic", + opParam = list(name = "dmormAD test 1"), + expr = quote({ + pld <- PDinverse_logdet(prec) + logProb <- dmnorm_inv_ld(x, mu, mat=prec, inv_ld=pld, prec_param=1, log=log) + out <- logProb + }), + args = list( + x = quote(double(1)), + mu = quote(double(1)), + prec = quote(double(2)), + log = quote(double(0)) + ), + outputType = quote(double(0)) + ), + argTypes = c(x='double(1)', mu='double(1)', prec='double(2)', log='double(0)'), + wrt = c('x','mu','mat'), + inputs = list(record = list(prec = prec, mu=mu, x=x, log=0), + test = list(prec = prec+0.1, mu=mu-0.07, x=x+0.03, log=1)) +) +dmnormAD_test2b_out <- test_AD2(dmnormAD_test2b) + +test_that("non-assignment of dmnormAD if cholesky param used", { + code <- nimbleCode({ + y[1:3] ~ dmnorm(mu[1:3], cholesky = L[1:3,1:3], prec_param = 0) + }) + expect_message(m <- nimbleModel(code, inits=list(mu=rep(0,3),L=diag(3))), + "Detected use of `cholesky` parameterization") +}) + +test_that("lifting of PDinverse_logdet", { + code <- nimbleCode({ + y[1:3] ~ dmnorm(mu[1:3], Q[2, 1:3, 5:7, 1]) + }) + m <- nimbleModel(code, data = list(y=rnorm(3))) + which_PDinverse <- grep("PDinverse_logdet", m$getNodeNames()) + expect_identical(length(m[[m$getNodeNames()[which_PDinverse]]]), 10L) +}) + +# getParam + +test_that('getParam', { + + code = nimbleCode({ + a[1:4] ~ dmnorm(mu[1:4],pr[1:4,1:4]) + }) + pr1 = diag(4) + pr1[1,2]=pr1[2,1]=.3 + pr1[1,3]=pr1[3,1]=-.1 + pr1[2,3]=pr1[3,2] = 0.4 + pr2 <- pr1 + pr1[1,2]=pr1[2,1]=.5 + pr1[1,3]=pr1[3,1]=-.2 + pr1[2,3]=pr1[3,2] = 0.3 + + m = nimbleModel(code, inits =list(mu=rep(1,4), pr = pr1)) + cm = compileNimble(m) + + cm$pr <- pr2 + cm$calculate(cm$getDependencies('pr')) + + expect_identical(pr1, m$getParam('a', 'prec')) + expect_identical(pr2, cm$getParam('a', 'prec')) + expect_equal(solve(pr1), m$getParam('a', 'cov')) + expect_equal(solve(pr2), cm$getParam('a', 'cov')) + + code = nimbleCode({ + a[1:3] ~ dmnorm(mu[1:3],cov=pr[1:3,1:3]) + }) + pr1 = diag(3) + pr1[1,2]=pr1[2,1]=.3 + pr2 <- pr1 + pr1[1,2]=pr1[2,1]=.5 + m = nimbleModel(code, inits =list(mu=rep(1,3), pr = pr1)) + cm = compileNimble(m) + + cm$pr <- pr2 + cm$calculate(cm$getDependencies('pr')) + + expect_identical(pr1, m$getParam('a', 'cov')) + expect_identical(pr2, cm$getParam('a', 'cov')) + expect_equal(solve(pr1), m$getParam('a', 'prec')) + expect_equal(solve(pr2), cm$getParam('a', 'prec')) +}) + +# test-models + +K <- 2 +y <- c(.1, .3) +model <- function() { + y[1:K] ~ dmnorm(mu[1:K], prec[1:K,1:K]); + for(i in 1:K) { + mu0[i] <- 0 + } + R[1,1] <- .01 + R[2,2] <- .01 + R[1,2] <- 0 + R[2,1] <- 0 + Omega[1,1] <- .01 + Omega[2,2] <- .01 + Omega[1,2] <- 0 + Omega[2,1] <- 0 + + mu[1:K] ~ dmnorm(mu0[1:K], Omega[1:K,1:K]) + prec[1:K,1:K] ~ dwish(R[1:K,1:K], 5) +} + +inits <- list(mu = c(0,0), prec = matrix(c(.005,.001,.001,.005), 2)) +data <- list(K = K, y = y) + +testBUGSmodel(example = 'testN', dir = "", + model = model, data = data, inits = inits, + useInits = TRUE) + +# Polya-gamma + +test_that('polyagamma validity checks', { + ## Various valid basic model structures + code <- nimbleCode({ + for(i in 1:n) { + p0[i] <- expit(b0+b1*x[i]) + y0[i] ~ dbern(p0[i]) + + logit(p1[i]) <- a0 + a1*x[i] + y1[i]~dbern(p1[i]) + + logit(p2[i]) <- inprod(x2[i,1:p], b2[1:p]) + y2[i]~dbin(p2[i], size = m) + + logit(p3[i]) <- x2[i,1:p] %*% b3[1:p] + y3[i]~dbin(p3[i], 1) + } + m ~ dpois(5) + + b0~dnorm(0, sd=10) + b1~dnorm(0, sd=10) + a0~dnorm(0, sd=10) + a1~dnorm(0, sd=10) + + b2[1:p] ~ dmnorm(mu[1:p], Q[1:p,1:p]) + for(i in 1:p) + b3[i] ~ dnorm(0, sd = 10) + }) + + n <- 10 + p <- 3 + constants <- list(n=n, p=p, x = runif(n), x2 = matrix(runif(n*p),n)) + ys <- rep(1,n) + data <- list(y0 = ys, y1 = ys, y2 = ys, y3 = ys) + + m <- nimbleModel(code, constants = constants, data = data) + conf <- configureMCMC(m, nodes = 'm') + expect_silent(conf$addSampler(type='polyagamma', target=c('b0','b1'))) + expect_silent(conf$addSampler(type='polyagamma', target=c('a0','a1'))) + expect_silent(conf$addSampler(type='polyagamma', target=c('b2'))) + expect_silent(conf$addSampler(type='polyagamma', target=c('b3'))) + expect_silent(mcmc <- buildMCMC(conf)) + + + ## dnorm + dmnorm + code <- nimbleCode({ + for(i in 1:n) { + p0[i] <- expit(b[1]+u[k[i]]+b[2]*x[i]) + y0[i] ~ dbern(p0[i]) + } + + b[1:2] ~ dmnorm(z[1:2], pr[1:2,1:2]) + + for(i in 1:p) + u[i] ~ dnorm(0,1) + }) + + + n <- 9 + p <- 3 + constants <- list(n=n, p=p, x = runif(n), k = rep(1:3, each = 3), z=rep(0,2), pr=diag(2)) + ys <- rep(1,n) + data <- list(y0 = ys) + + m <- nimbleModel(code, constants = constants, data = data) + conf <- configureMCMC(m, nodes = NULL) + expect_silent(conf$addSampler(type='polyagamma', target=c('b','u'))) + expect_silent(mcmc <- buildMCMC(conf)) + mcmc$run(1) + expect_equal(mcmc$samplerFunctions[[1]]$X, + cbind(rep(1,n), constants$x, c(rep(1,3),rep(0,6)), + c(rep(0,3),rep(1,3),rep(0,3)), + c(rep(0,6),rep(1,3)))) + lens <- c(2,1,1,1); names(lens) <- c('b[1:2]','u[1]','u[2]','u[3]') + expect_equal(mcmc$samplerFunctions[[1]]$nodeLengths, lens) + +}) + +test_that('Wishart-dmnorm conjugacy', { + n <- 3 + R <- diag(rep(1,3)) + mu <- 1:3 + Y <- matrix(rnorm(9), 3) + data <- list(Y = Y, mu = mu, R = R) + code <- nimbleCode( { + for(i in 1:3) { + Y[i, 1:3] ~ dmnorm(mu[1:3], Omega[1:3,1:3]); + } + Omega[1:3,1:3] ~ dwish(R[1:3,1:3], 4); + }) + m <- nimbleModel(code, data = data) + conf <- configureMCMC(m) + expect_identical(conf$getSamplers()[[1]]$name, "conjugate_dwish_dmnormAD_identity") + + code <- nimbleCode( { + for(i in 1:3) { + Y[i, 1:3] ~ dmnorm(mu[1:3], cov = Omega[1:3,1:3]); + } + Omega[1:3,1:3] ~ dwish(R[1:3,1:3], 4); + }) + m <- nimbleModel(code, data = data) + conf <- configureMCMC(m) + expect_identical(conf$getSamplers()[[1]]$name, "RW_wishart") + + code <- nimbleCode( { + for(i in 1:3) { + Y[i, 1:3] ~ dmnorm(mu[1:3], Omega[1:3,1:3]); + } + Omega[1:3,1:3] ~ dinvwish(R[1:3,1:3], 4); + }) + m <- nimbleModel(code, data = data) + conf <- configureMCMC(m) + expect_identical(conf$getSamplers()[[1]]$name, "RW_wishart") + + code <- nimbleCode( { + for(i in 1:3) { + Y[i, 1:3] ~ dmnorm(mu[1:3], cov = Omega[1:3,1:3]); + } + Omega[1:3,1:3] ~ dinvwish(R[1:3,1:3], 4); + }) + m <- nimbleModel(code, data = data) + conf <- configureMCMC(m) + expect_identical(conf$getSamplers()[[1]]$name, "conjugate_dinvwish_dmnormAD_identity") +}) + +# Run various NIMBLE tests that use dmnorm with dmnormAD instead. + +# From test-mcmc.R + +## verbose: set to FALSE +nimbleVerboseSetting <- nimbleOptions('verbose') +nimbleOptions(verbose = FALSE) + +## MCMC progress bar: set to FALSE +nimbleProgressBarSetting <- nimbleOptions('MCMCprogressBar') +nimbleOptions(MCMCprogressBar = FALSE) + +## MCMC orderSamplersPosteriorPredictiveLast - save current setting +nimbleReorderPPsamplersSetting <- getNimbleOption('MCMCorderPosteriorPredictiveSamplersLast') + +## MCMC use usePosteriorPredictiveSampler - save current setting +nimbleUsePosteriorPredictiveSamplerSetting <- getNimbleOption('MCMCusePosteriorPredictiveSampler') + +## MCMC calculation include predictive dependencies - save current setting +nimbleUsePredictiveDependenciesSetting <- nimbleOptions('MCMCusePredictiveDependenciesInCalculations') + +## MCMC warn about unsampled nodes - save current setting +nimbleWarnUnsampledNodesSetting <- nimbleOptions('MCMCwarnUnsampledStochasticNodes') + +test_that('elliptical slice sampler setup', { + set.seed(0) + ESScode <- quote({ + x[1:d] ~ dmnorm(mu_x[1:d], prec = prec_x[1:d, 1:d]) + y[1:d] ~ dmnorm(x[1:d], prec = prec_y[1:d, 1:d]) + }) + d <- 3 + mu_x <- rnorm(d) + temp <- array(rnorm(d^2), c(d,d)) + prec_x <- solve(temp %*% t(temp)) + temp <- array(rnorm(d^2), c(d,d)) + prec_y <- solve(temp %*% t(temp)) + y <- rnorm(d) + ESSconstants <- list(d = d, mu_x = mu_x, prec_x = prec_x, prec_y = prec_y) + ESSdata <- list(y = y) + ESSinits <- list(x = rep(0, d)) + + test_mcmc(model = ESScode, data = c(ESSconstants, ESSdata), inits = ESSinits, + name = 'exact values of elliptical slice sampler', + seed = 0, + exactSample = list(`x[1]` = c(-0.492880566939352, -0.214539223107114, 1.79345037297218, 1.17324496091208, 2.14095077672555, 1.60417482445964, 1.94196916651627, 2.66737323347255, 2.66744178776022, 0.253966883192744), `x[2]` = c(-0.161210109217102, -0.0726534676226932, 0.338308532423757, -0.823652445515156, -0.344130712698579, -0.132642244861469, -0.0253168895009594, 0.0701624130921676, 0.0796842215444978, -0.66369112443311), `x[3]` = c(0.278627475932455, 0.0661336950029345, 0.407055002920732, 1.98761228946318, 1.0839897275519, 1.00262648370199, 0.459841485268785, 2.59229443025387, 1.83769567435409, 1.92954706515119)), + samplers = list(list(type = 'ess', target = 'x'))) + + test_mcmc(model = ESScode, data = c(ESSconstants, ESSdata), inits = ESSinits, + name = 'results to tolerance of elliptical slice sampler', + results = list(mean = list(x = c(1.0216463, -0.4007247, 1.1416904))), + resultsTolerance = list(mean = list(x = c(0.01, 0.01, 0.01))), + numItsC = 100000, + samplers = list(list(type = 'ess', target = 'x')), avoidNestedTest = TRUE) + + }) + +test_that('block sampler on MVN node setup', { + code <- nimbleCode({ + mu[1] <- 10 + mu[2] <- 20 + mu[3] <- 30 + x[1:3] ~ dmnorm(mu[1:3], prec = Q[1:3,1:3]) + }) + + Q = matrix(c(1.0,0.2,-1.0,0.2,4.04,1.6,-1.0,1.6,10.81), nrow=3) + data = list(Q = Q) + inits = list(x = c(10, 20, 30)) + + test_mcmc(model = code, name = 'block sampler on multivariate node', data = data, seed = 0, numItsC = 10000, + results = list(mean = list(x = c(10,20,30)), + var = list(x = diag(solve(Q)))), + resultsTolerance = list(mean = list(x = rep(1,3)), + var = list(x = c(.1, .03, .01))), + samplers = list( + list(type = 'RW_block', target = 'x[1:3]')), avoidNestedTest = TRUE) + # caution: setting targetNodes='x' works but the initial end sampler is not removed because x[1:3] in targetNode in default sampler != 'x' in targetNodes passed in + if(FALSE) { + Rmodel <- nimbleModel(code, constants = list(Q=Q)) + mcmcspec <- MCMCspec(Rmodel, nodes = NULL) + mcmcspec$addSampler(type = 'RW_block', target = 'x', control = list(adaptInterval=500)) + mcmcspec$getMonitors() + Rmcmc <- buildMCMC(mcmcspec) + Cmodel <- compileNimble(Rmodel) + Cmcmc <- compileNimble(Rmcmc, project = Rmodel) + Cmcmc(200000) ## this runs nearly instantaneously on my computer -DT + samples <- as.matrix(nfVar(Cmcmc, 'mvSamples')) + samples <- samples[50001:200000,] + dim(samples) + apply(samples, 2, mean) + solve(Q) + cov(samples) + propCov <- nfVar(Cmcmc, 'samplerFunctions')[[1]]$propCov + scale <- nfVar(Cmcmc, 'samplerFunctions')[[1]]$scale + propCov * scale^2 + + nfVar(Cmcmc, 'samplerFunctions')[[1]]$scaleHistory + nfVar(Cmcmc, 'samplerFunctions')[[1]]$acceptanceRateHistory + nfVar(Cmcmc, 'samplerFunctions')[[1]]$scale + nfVar(Cmcmc, 'samplerFunctions')[[1]]$propCov + ## why is the proposal cov w/ .99 cross-corrs? + ## also MCMC in C takes a surprisingly long time - this might be threaded lin alg behaving badly on small matrices + } +}) + +test_that('second block sampler on multivariate node', { +### DT's model + mu <- c(1,2,3) + corr <- matrix(c(1,.8,0.3,.8,1,0,0.3,0,1), nrow=3) + varr <- c(1,2,3) + Sig <- diag(sqrt(varr)) + Q <- Sig %*% corr %*% Sig + P <- solve(Q) + + code <- nimbleCode({ + x[1:3] ~ dmnorm(mu[1:3], prec = P[1:3,1:3]) + }) + data = list(P = P, mu = mu) + + test_mcmc(model = code, name = 'second block sampler on multivariate node', data = data, seed = 0, numItsC = 100000, + results = list(mean = list(x = mu), + var = list(x = varr)), + resultsTolerance = list(mean = list(x = rep(.1,3)), + var = list(x = c(.1,.1,.1))), + samplers = list( + list(type = 'RW_block', target = 'x[1:3]')), avoidNestedTest = TRUE) + }) + +test_that('MVN conjugate setup', { +### MVN conjugate update + + set.seed(0) + mu0 = 1:3 + Q0 = matrix(c(1, .2, .8, .2, 2, 1, .8, 1, 2), nrow = 3) + Q = solve(matrix(c(3, 1.7, .9, 1.7, 2, .6, .9, .6, 1), nrow = 3)) + a = c(-2, .5, 1) + B = matrix(rnorm(9), 3) + +##### not currently working - see Perry's email of ~ 10/6/14 + ## code <- nimbleCode({ + ## mu[1:3] ~ dmnorm(mu0[1:3], Q0[1:3, 1:3]) + ## y[1:3] ~ dmnorm(asCol(a[1:3]) + B[1:3, 1:3] %*% asCol(mu[1:3]), Q[1:3, 1:3]) + ## }) + + code <- nimbleCode({ + mu[1:3] ~ dmnorm(mu0[1:3], Q0[1:3, 1:3]) + y_mean[1:3] <- asCol(a[1:3]) + B[1:3, 1:3] %*% asCol(mu[1:3]) + y[1:3] ~ dmnorm(y_mean[1:3], Q[1:3, 1:3]) + }) + + + mu <- mu0 + chol(solve(Q0)) %*% rnorm(3) + # make sure y is a vec not a 1-col matrix or get a dimensionality error + y <- c(a + B%*%mu + chol(solve(Q)) %*% rnorm(3)) + data = list(mu0 = mu0, Q0 = Q0, Q = Q, a = a, B = B, y = y) + + muQtrue = t(B) %*% Q%*%B + Q0 + muMeanTrue = c(solve(muQtrue, crossprod(B, Q%*%(y-a)) + Q0%*%mu0)) + + test_mcmc(model = code, name = 'two-level multivariate normal', data = data, seed = 0, numItsC = 10000, + results = list(mean = list(mu = muMeanTrue), + cov = list(mu = solve(muQtrue))), + resultsTolerance = list(mean = list(mu = rep(.02,3)), + cov = list(mu = matrix(.01, 3, 3))), avoidNestedTest = TRUE) + + +### scalar RW updates in place of conjugate mv update + + test_mcmc(model = code, name = 'two-level multivariate normal with scalar updaters', data = data, seed = 0, numItsC = 100000, + results = list(mean = list(mu = muMeanTrue), + cov = list(mu = solve(muQtrue))), + resultsTolerance = list(mean = list(mu = rep(.03,3)), + cov = list(mu = matrix(.03, 3, 3))), + samplers = list(list(type = 'RW', target = 'mu[1]'), + list(type = 'RW', target = 'mu[2]'), + list(type = 'RW', target = 'mu[3]')), + removeAllDefaultSamplers = TRUE, avoidNestedTest = TRUE) + +}) + + +test_that('another MVN conjugate sampler setup', { + set.seed(0) + prior_mean <- rep(0,5) + tmp <- array(rnorm(25), c(5,5)) + tmp <- tmp + t(tmp) + 5*diag(5) + prior_cov <- tmp + a <- array(rnorm(20), c(4,5)) + B <- array(NA, c(4,5,5)) + for(i in c(2,4)) B[i,,] <- array(rnorm(25), c(5,5)) + B[1,,] <- diag(5) + B[3,,] <- diag(5) + M_y <- array(NA, c(4,5,5)) + for(i in 1:4) { + tmp <- array(rnorm(25,i), c(5,5)) + tmp <- tmp + t(tmp) + 5*i*diag(5) + M_y[i,,] <- tmp + } + x <- rep(0, 5) + y <- array(rnorm(20), c(4,5)) + + code <- nimbleCode({ + x[1:5] ~ dmnorm(mean = prior_mean[1:5], cov = prior_cov[1:5,1:5]) + for(i in 1:4) + mu_y[i,1:5] <- asCol(a[i,1:5]) + B[i,1:5,1:5] %*% asCol(x[1:5]) + y[1,1:5] ~ dmnorm(mu_y[1,1:5], prec = M_y[1,1:5,1:5]) + y[2,1:5] ~ dmnorm(mu_y[2,1:5], cov = M_y[2,1:5,1:5]) + y[3,1:5] ~ dmnorm(mu_y[3,1:5], prec = M_y[3,1:5,1:5]) + y[4,1:5] ~ dmnorm(mu_y[4,1:5], cov = M_y[4,1:5,1:5]) + }) + constants <- list(prior_mean=prior_mean, prior_cov=prior_cov, a=a, B=B, M_y=M_y) + data <- list(y=y) + inits <- list(x=x) + Rmodel <- nimbleModel(code, constants, data, inits) + spec <- configureMCMC(Rmodel) + Rmcmc <- buildMCMC(spec) + + Cmodel <- compileNimble(Rmodel) + Cmcmc <- compileNimble(Rmcmc, project = Rmodel) + + set.seed(0) + Rmcmc$run(10) + Rsamples <- as.matrix(Rmcmc$mvSamples) + set.seed(0) + Cmcmc$run(10) + Csamples <- as.matrix(Cmcmc$mvSamples) + + expect_equal(round(as.numeric(Rsamples), 8), + c(0.97473128, 0.50438666, 1.1251132, 0.83830666, 0.74077066, 0.92935482, 0.83758372, 0.98708273, 1.24199937, 0.67348127, -0.54387714, -0.60713969, -0.51392796, -0.3176801, -0.34416529, -0.08530564, -0.47160157, -0.21996584, -0.20504917, -0.77287122, 0.78462584, 0.46103509, 0.43862813, 0.49343096, 0.61020864, 0.55088287, 0.53887202, 0.49863894, 0.62691318, 0.80142839, 0.34941152, 0.06623608, 0.05624477, 0.21369178, 0.26585415, -0.1439989, -0.03133488, 0.3544062, -0.03518959, 0.27415746, 0.40977, 0.8351078, 0.25719293, 0.05663917, 0.30894028, 0.33113315, 0.47647909, 0.26143962, 0.07180759, 0.27255767), + info = 'R sample not correct compared to known result') + + dif <- as.numeric(Rsamples - Csamples) + expect_true(max(abs(dif)) < 1E-15, info = 'R and C equiv') + + y_prec <- array(NA, c(4,5,5)) + y_prec[1,,] <- M_y[1,,] + y_prec[2,,] <- solve(M_y[2,,]) + y_prec[3,,] <- M_y[3,,] + y_prec[4,,] <- solve(M_y[4,,]) + contribution_mean <- array(NA, c(4,5)) + for(i in 1:4) contribution_mean[i,] <- t(B[i,,]) %*% y_prec[i,,] %*% (y[i,] - a[i,]) + contribution_prec <- array(NA, c(4,5,5)) + for(i in 1:4) contribution_prec[i,,] <- t(B[i,,]) %*% y_prec[i,,] %*% B[i,,] + prior_prec <- solve(prior_cov) + post_prec <- prior_prec + apply(contribution_prec, c(2,3), sum) + post_cov <- solve(post_prec) + post_mean <- (post_cov %*% (prior_prec %*% prior_mean + apply(contribution_mean, 2, sum)))[,1] + + Cmcmc$run(100000) + Csamples <- as.matrix(Cmcmc$mvSamples) + + dif_mean <- as.numeric(apply(Csamples, 2, mean)) - post_mean + expect_true(max(abs(dif_mean)) < 0.001, info = 'posterior mean') + + dif_cov <- as.numeric(cov(Csamples) - post_cov) + expect_true(max(abs(dif_cov)) < 0.001, info = 'posterior cov') +}) + + +test_that('conjugate Wishart setup', { + set.seed(0) + + trueCor <- matrix(c(1, .3, .7, .3, 1, -0.2, .7, -0.2, 1), 3) + covs <- c(3, 2, .5) + + trueCov = diag(sqrt(covs)) %*% trueCor %*% diag(sqrt(covs)) + Omega = solve(trueCov) + + n = 20 + R = diag(rep(1,3)) + mu = 1:3 + Y = mu + t(chol(trueCov)) %*% matrix(rnorm(3*n), ncol = n) + M = 3 + data <- list(Y = t(Y), n = n, M = M, mu = mu, R = R) + + code <- nimbleCode( { + for(i in 1:n) { + Y[i, 1:M] ~ dmnorm(mu[1:M], Omega[1:M,1:M]); + } + Omega[1:M,1:M] ~ dwish(R[1:M,1:M], 4); + }) + + newDf = 4 + n + newR = R + tcrossprod(Y- mu) + OmegaTrueMean = newDf * solve(newR) + + wishRV <- array(0, c(M, M, 10000)) + for(i in 1:10000) { + z <- t(chol(solve(newR))) %*% matrix(rnorm(3*newDf), ncol = newDf) + wishRV[ , , i] <- tcrossprod(z) + } + OmegaSimTrueSDs = apply(wishRV, c(1,2), sd) + + test_mcmc(model = code, name = 'conjugate Wishart', data = data, seed = 0, numItsC = 1000, inits = list(Omega = OmegaTrueMean), + results = list(mean = list(Omega = OmegaTrueMean ), + sd = list(Omega = OmegaSimTrueSDs)), + resultsTolerance = list(mean = list(Omega = matrix(.05, M,M)), + sd = list(Omega = matrix(0.06, M, M))), avoidNestedTest = TRUE) + # issue with Chol in R MCMC - probably same issue as in jaw-linear + +}) + +test_that('conjugate Wishart setup with scaling', { + set.seed(0) + + trueCor <- matrix(c(1, .3, .7, .3, 1, -0.2, .7, -0.2, 1), 3) + covs <- c(3, 2, .5) + tau <- 4 + + trueCov = diag(sqrt(covs)) %*% trueCor %*% diag(sqrt(covs)) + Omega = solve(trueCov) / tau + + n = 20 + R = diag(rep(1,3)) + mu = 1:3 + Y = mu + t(chol(trueCov)) %*% matrix(rnorm(3*n), ncol = n) + M = 3 + data <- list(Y = t(Y), n = n, M = M, mu = mu, R = R) + + code <- nimbleCode( { + for(i in 1:n) { + Y[i, 1:M] ~ dmnorm(mu[1:M], tmp[1:M,1:M]) + } + tmp[1:M,1:M] <- tau * Omega[1:M,1:M] + Omega[1:M,1:M] ~ dwish(R[1:M,1:M], 4); + }) + + newDf = 4 + n + newR = R + tcrossprod(Y - mu)*tau + OmegaTrueMean = newDf * solve(newR) + + wishRV <- array(0, c(M, M, 10000)) + for(i in 1:10000) { + z <- t(chol(solve(newR))) %*% matrix(rnorm(3*newDf), ncol = newDf) + wishRV[ , , i] <- tcrossprod(z) + } + OmegaSimTrueSDs = apply(wishRV, c(1,2), sd) + + m <- nimbleModel(code, data = data[1],constants=data[2:5],inits = list(Omega = OmegaTrueMean, tau = tau)) + conf <- configureMCMC(m) + expect_equal(conf$getSamplers()[[1]]$name, "conjugate_dwish_dmnormAD_multiplicativeScalar", + info = "conjugate dmnormAD-dwish with scaling not detected") + + test_mcmc(model = code, name = 'conjugate Wishart, scaled', data = data, seed = 0, numItsC = 1000, inits = list(Omega = OmegaTrueMean, tau = tau), + results = list(mean = list(Omega = OmegaTrueMean ), + sd = list(Omega = OmegaSimTrueSDs)), + resultsTolerance = list(mean = list(Omega = matrix(.05, M,M)), + sd = list(Omega = matrix(0.06, M, M))), avoidNestedTest = TRUE) + # issue with Chol in R MCMC - probably same issue as in jaw-linear + +}) + +test_that('using RW_wishart sampler on non-conjugate Wishart node', { + set.seed(0) + trueCor <- matrix(c(1, .3, .7, .3, 1, -0.2, .7, -0.2, 1), 3) + covs <- c(3, 2, .5) + trueCov = diag(sqrt(covs)) %*% trueCor %*% diag(sqrt(covs)) + Omega = solve(trueCov) + n = 20 + R = diag(rep(1,3)) + mu = 1:3 + Y = mu + t(chol(trueCov)) %*% matrix(rnorm(3*n), ncol = n) + M = 3 + data <- list(Y = t(Y), n = n, M = M, mu = mu, R = R) + code <- nimbleCode( { + for(i in 1:n) { + Y[i, 1:M] ~ dmnorm(mu[1:M], Omega[1:M,1:M]) + } + Omega[1:M,1:M] ~ dwish(R[1:M,1:M], 4) + }) + newDf = 4 + n + newR = R + tcrossprod(Y- mu) + OmegaTrueMean = newDf * solve(newR) + wishRV <- array(0, c(M, M, 10000)) + for(i in 1:10000) { + z <- t(chol(solve(newR))) %*% matrix(rnorm(3*newDf), ncol = newDf) + wishRV[ , , i] <- tcrossprod(z) + } + OmegaSimTrueSDs = apply(wishRV, c(1,2), sd) + allData <- data + constants <- list(n = allData$n, M = allData$M, mu = allData$mu, R = allData$R) + data <- list(Y = allData$Y) + inits <- list(Omega = OmegaTrueMean) + + Rmodel <- nimbleModel(code, constants, data, inits) + Rmodel$calculate() + conf <- configureMCMC(Rmodel, nodes = NULL) + conf$addSampler('Omega', 'RW_wishart') + Rmcmc <- buildMCMC(conf) + compiledList <- compileNimble(list(model=Rmodel, mcmc=Rmcmc)) + Cmodel <- compiledList$model; Cmcmc <- compiledList$mcmc + set.seed(0) + samples <- runMCMC(Cmcmc, 10000) + d1 <- as.numeric(apply(samples, 2, mean)) - as.numeric(OmegaTrueMean) + difference <- sum(round(d1,9) - c(0.024469145, -0.011872571, -0.045035297, -0.011872571, -0.003443918, 0.009363410, -0.045035297, 0.009363410, 0.049971420)) + expect_true(difference == 0) +}) + + +test_that('using RW_wishart sampler on inverse-Wishart distribution', { + code <- nimbleCode( { + for(i in 1:n) { + Y[i, 1:M] ~ dmnorm(mu[1:M], cov = C[1:M,1:M]) + } + C[1:M,1:M] ~ dinvwish(R[1:M,1:M], 4) + }) + + set.seed(0) + trueCor <- matrix(c(1, .3, .7, .3, 1, -0.2, .7, -0.2, 1), 3) + covs <- c(3, 2, .5) + trueCov = diag(sqrt(covs)) %*% trueCor %*% diag(sqrt(covs)) + n = 20 + R = diag(rep(1,3)) + mu = 1:3 + Y = mu + t(chol(trueCov)) %*% matrix(rnorm(3*n), ncol = n) + M = 3 + data <- list(Y = t(Y), n = n, M = M, mu = mu, R = R) + + newDf = 4 + n + newR = R + tcrossprod(Y- mu) + CTrueMean <- newR / newDf + CTrueMeanVec <- as.numeric(CTrueMean) + allData <- data + constants <- list(n = allData$n, M = allData$M, mu = allData$mu, R = allData$R) + data <- list(Y = allData$Y) + inits <- list(C = CTrueMean) + niter <- 50000 + + Rmodel <- nimbleModel(code, constants, data, inits) + Rmodel$calculate() + conf <- configureMCMC(Rmodel) + Rmcmc <- buildMCMC(conf) + compiledList <- compileNimble(list(model=Rmodel, mcmc=Rmcmc)) + Cmodel <- compiledList$model; Cmcmc <- compiledList$mcmc + set.seed(0) + samples <- runMCMC(Cmcmc, niter) + conjMean <- as.numeric(apply(samples, 2, mean)) + conjSD <- as.numeric(apply(samples, 2, sd)) + + Rmodel <- nimbleModel(code, constants, data, inits) + Rmodel$calculate() + conf <- configureMCMC(Rmodel, nodes = NULL) + conf$addSampler('C', 'RW_wishart') + Rmcmc <- buildMCMC(conf) + compiledList <- compileNimble(list(model=Rmodel, mcmc=Rmcmc)) + Cmodel <- compiledList$model; Cmcmc <- compiledList$mcmc + set.seed(0) + samples <- runMCMC(Cmcmc, niter) + RWMean <- as.numeric(apply(samples, 2, mean)) + RWSD <- as.numeric(apply(samples, 2, sd)) + + expect_true(all(round(as.numeric(RWMean - conjMean), 9) == c(-0.001651758, -0.009675571, 0.004894809, -0.009675571, 0.015533882, -0.008256095, 0.004894809, -0.008256095, 0.002119615))) + expect_true(all(round(as.numeric(RWSD - conjSD), 9) == c(0.022803503, -0.010107015, 0.012342044, -0.010107015, 0.006191412, -0.000091101, 0.012342044, -0.000091101, 0.001340032))) +}) + +test_that('detect conjugacy when scaling Wishart, inverse Wishart cases', { + code <- nimbleCode({ + mycov[1:p, 1:p] <- lambda * Sigma[1:p,1:p] / eta + y[1:p] ~ dmnorm(z[1:p], cov = mycov[1:p, 1:p]) + Sigma[1:p,1:p] ~ dinvwish(S[1:p,1:p], nu) + }) + m <- nimbleModel(code, constants = list(p = 3)) + expect_identical(length(m$checkConjugacy('Sigma')), 1L, 'inverse-Wishart case') + + code <- nimbleCode({ + mycov[1:p, 1:p] <- lambda * Sigma[1:p,1:p] / eta + y[1:p] ~ dmnorm(z[1:p], prec = mycov[1:p, 1:p]) + Sigma[1:p,1:p] ~ dwish(S[1:p,1:p], nu) + }) + m <- nimbleModel(code, constants = list(p = 3)) + expect_identical(length(m$checkConjugacy('Sigma')), 1L, 'Wishart case') + + code <- nimbleCode({ + mycov[1:p, 1:p] <- lambda * Sigma[1:p,1:p] / eta + y[1:p] ~ dmnorm(z[1:p], cov = mycov[1:p, 1:p]) + Sigma[1:p,1:p] ~ dwish(S[1:p,1:p], nu) + }) + m <- nimbleModel(code, constants = list(p = 3)) + expect_identical(length(m$checkConjugacy('Sigma')), 0L, 'inverse-Wishart not-conj case') + + code <- nimbleCode({ + mycov[1:p, 1:p] <- lambda * Sigma[1:p,1:p] / eta + y[1:p] ~ dmnorm(z[1:p], prec = mycov[1:p, 1:p]) + Sigma[1:p,1:p] ~ dinvwish(S[1:p,1:p], nu) + }) + m <- nimbleModel(code, constants = list(p = 3)) + expect_identical(length(m$checkConjugacy('Sigma')), 0L, 'Wishart not-conj case') + + code <- nimbleCode({ + mycov[1:p, 1:p] <- lambda[1:p,1:p] * Sigma[1:p,1:p] + y[1:p] ~ dmnorm(z[1:p], cov = mycov[1:p, 1:p]) + Sigma[1:p,1:p] ~ dinvwish(S[1:p,1:p], nu) + }) + m <- nimbleModel(code, constants = list(p = 3)) + expect_identical(length(m$checkConjugacy('Sigma')), 0L, 'Wishart case') + + code <- nimbleCode({ + mycov[1:p, 1:p] <- lambda[1:p,1:p] %*% Sigma[1:p,1:p] + y[1:p] ~ dmnorm(z[1:p], cov = mycov[1:p, 1:p]) + Sigma[1:p,1:p] ~ dinvwish(S[1:p,1:p], nu) + }) + m <- nimbleModel(code, constants = list(p = 3)) + expect_identical(length(m$checkConjugacy('Sigma')), 0L, 'Wishart case') +}) + + +test_that("realized conjugacy links are working", { + + ## dmnorm cases + code <- nimbleCode({ + mu[1:3] ~ dmnorm(z[1:3], pr[1:3,1:3]) + for(i in 1:2) + y1[i, 1:3] ~ dmnorm(mu[1:3], pr[1:3,1:3]) + mn2[1:3] <- b0[1:3] + mu[1:3] + for(i in 1:2) { + y2[i, 1:3] ~ dmnorm(mn2[1:3], pr[1:3,1:3]) + } + mn3[1:3] <- A[1:3,1:3]%*%mu[1:3] + for(i in 1:2) { + y3[i, 1:3] ~ dmnorm(mn3[1:3], pr[1:3,1:3]) + } + mn4[1:3] <- b0[1:3] + A[1:3,1:3]%*%mu[1:3] + for(i in 1:2) { + y4[i, 1:3] ~ dmnorm(mn4[1:3], pr[1:3,1:3]) + } + }) + m <- nimbleModel(code, data = list (y1 = matrix(rnorm(6),2), + y2 = matrix(rnorm(6),2), + y3 = matrix(rnorm(6),2), + y4 = matrix(rnorm(6),2)), + inits = list(b0 = rnorm(3), A=matrix(1:9, 3), pr = diag(3))) + conf <- configureMCMC(m) + mcmc <- buildMCMC(conf) + + expect_identical(conf$getSamplers()[[1]]$name, "conjugate_dmnormAD_dmnormAD_additive_dmnormAD_identity_dmnormAD_linear_dmnormAD_multiplicative") + expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnormAD_identity, 2L) + expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnormAD_additive, 2L) + expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnormAD_multiplicative, 2L) + expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnormAD_linear, 2L) + + expect_identical(c('dep_dmnormAD_identity_coeff', 'dep_dmnormAD_additive_coeff') %in% + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + expect_identical(c('dep_dmnormAD_multiplicative_coeff', 'dep_dmnormAD_linear_coeff') %in% + ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) + expect_identical(c('dep_dmnormAD_identity_offset', 'dep_dmnormAD_multiplicative_offset') %in% + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + expect_identical(c('dep_dmnormAD_additive_offset', 'dep_dmnormAD_linear_offset') %in% + ls(mcmc$samplerFunctions[[1]]), rep(TRUE, 2)) + + expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnormAD_identity_nodeNames, c('y1[1, 1:3]', 'y1[2, 1:3]')) + expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnormAD_additive_nodeNames, c('y2[1, 1:3]', 'y2[2, 1:3]')) + expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnormAD_multiplicative_nodeNames, c('y3[1, 1:3]', 'y3[2, 1:3]')) + expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnormAD_linear_nodeNames, c('y4[1, 1:3]', 'y4[2, 1:3]')) + + ## dwish case + code <- nimbleCode({ + for(i in 1:2) { + y1[i, 1:3] ~ dmnorm(mu[1:3], pr[1:3,1:3]) + } + pr2[1:3,1:3] <- d*pr[1:3,1:3] + for(i in 1:2) { + y2[i, 1:3] ~ dmnorm(mu[1:3], pr2[1:3,1:3]) + } + pr[1:3,1:3] ~ dwish(R[1:3,1:3], 8) + }) + m <- nimbleModel(code, data = list (y1 = matrix(rnorm(6),2), + y2 = matrix(rnorm(6),2)), + inits = list(pr = diag(3), R = diag(3))) + conf <- configureMCMC(m) + mcmc <- buildMCMC(conf) + + expect_identical(conf$getSamplers()[[1]]$name, "conjugate_dwish_dmnormAD_identity_dmnormAD_multiplicativeScalar") + expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnormAD_identity, 2L) + expect_identical(mcmc$samplerFunctions[[1]]$N_dep_dmnormAD_multiplicativeScalar, 2L) + + expect_identical('dep_dmnormAD_identity_coeff' %in% + ls(mcmc$samplerFunctions[[1]]), FALSE) + expect_identical('dep_dmnormAD_multiplicativeScalar_coeff' %in% + ls(mcmc$samplerFunctions[[1]]), TRUE) + expect_identical(c('dep_dmnormAD_identity_offset', 'dep_dmnormAD_multiplicativeScalar_offset') %in% + ls(mcmc$samplerFunctions[[1]]), rep(FALSE, 2)) + + expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnormAD_identity_nodeNames, c('y1[1, 1:3]', 'y1[2, 1:3]')) + expect_identical(mcmc$samplerFunctions[[1]]$dep_dmnormAD_multiplicativeScalar_nodeNames, c('y2[1, 1:3]', 'y2[2, 1:3]')) +}) + + +options(warn = RwarnLevel) +nimbleOptions(verbose = nimbleVerboseSetting) +nimbleOptions(MCMCprogressBar = nimbleProgressBarSetting) +nimbleOptions(buildModelDerivs = BMDopt) diff --git a/packages/nimble/tests/testthat/test-ADmodels.R b/packages/nimble/tests/testthat/test-ADmodels.R index 6732ed888..3632d9754 100644 --- a/packages/nimble/tests/testthat/test-ADmodels.R +++ b/packages/nimble/tests/testthat/test-ADmodels.R @@ -5,9 +5,9 @@ nimbleOptions(enableDerivs = TRUE) nimbleOptions(buildModelDerivs = TRUE) nimbleOptions(useADcholAtomic = TRUE) -nimbleOptions(useADsolveAtomic = TRUE) -nimbleOptions(useADmatMultAtomic = TRUE) -nimbleOptions(useADmatInverseAtomic = TRUE) +nimbleOptions(useADsolveAtomic = TRUE) +nimbleOptions(useADmatMultAtomic = TRUE) +nimbleOptions(useADmatInverseAtomic = TRUE) relTol <- eval(formals(test_ADModelCalculate)$relTol) relTol[3] <- 1e-6 @@ -58,6 +58,8 @@ test_that('AD support is correctly determined for distributions', { c ~ userObj1$dNO() d ~ userObj1$dYES() }) + temporarilyAssignInGlobalEnv(userObj1) + m <- nimbleModel(mc) expect_false(m$modelDef$checkADsupportForDistribution("duserNO")) expect_true(m$modelDef$checkADsupportForDistribution("duserYES")) @@ -102,7 +104,7 @@ test_that('pow and pow_int work', { ## 28 sec. expect_equal(wrapperDerivs$jacobian, cDerivs$jacobian) expect_equal(wrapperDerivs$hessian, cDerivs$hessian, tolerance = 1e-4) - + mc <- nimbleCode({ d ~ dnorm(0, sd = 10) trt <- -1 @@ -153,7 +155,7 @@ test_that('makeModelDerivsInfo works correctly', { inits = list(mu = rnorm(n), sigma2 = 2.1, tau = 1.7, mu0 = 0.8)) ## Various wrt and calcNodes cases; not exhaustive, but hopefully capturing most possibilities - + ## distinct wrt and calcNodes: updateNodes should have 'mu' and determ parents not in calcNodes result <- makeModelDerivsInfo(model = model, wrtNodes = c('sigma2', 'mu0', 'tau'), calcNodes = c('mu','y')) expect_identical(result$updateNodes, c("lifted_sqrt_oPsigma2_cP", "lifted_d1_over_sqrt_oPtau_cP", @@ -266,7 +268,7 @@ test_that('makeModelDerivsInfo works correctly', { result <- makeModelDerivsInfo(model = model, wrtNodes = c('sigma2', 'mu0', 'mu'), calcNodes = c('sigma2', 'mu0', 'mu')) expect_identical(result$updateNodes, c(lftChElems)) expect_identical(result$constantNodes, character(0)) - + }) ## basic model, with lifted nodes @@ -278,7 +280,7 @@ code <- nimbleCode({ tau0 ~ dgamma(1.5, 2.5) tau ~ dgamma(1.2, 2.3) for(j in 1:3) { - for(i in 1:2) + for(i in 1:2) a[j, i] ~ dnorm(mu[j], var = tau) mu[j] ~ dnorm(mu0, var = tau0) } @@ -293,7 +295,7 @@ test_ADModelCalculate(model, relTol = relTolTmp, checkCompiledValuesIdentical = verbose = verbose, name = 'basic model, lifted nodes') ## basic state space model -set.seed(1) +set.seed(1) code <- nimbleCode({ x0 ~ dnorm(0,1) x[1] ~ dnorm(x0, 1) @@ -313,7 +315,7 @@ relTolTmp[2] <- 1e-7 relTolTmp[3] <- 1e-4 relTolTmp[4] <- 1e-2 ## 204 sec. -test_ADModelCalculate(model, relTol = relTolTmp, useFasterRderivs = TRUE, verbose = verbose, name = 'basic state space') +test_ADModelCalculate(model, relTol = relTolTmp, useFasterRderivs = TRUE, verbose = verbose, name = 'basic state space') ## basic tricky indexing set.seed(1) @@ -328,7 +330,7 @@ model <- nimbleModel(code, dimensions = list(x = 2, y = 2, z = 3), inits = list( relTolTmp <- relTol relTolTmp[4] <- 1e-3 ### 185 sec. -test_ADModelCalculate(model, relTol = relTolTmp, verbose = verbose, name = 'basic tricky indexing', +test_ADModelCalculate(model, relTol = relTolTmp, verbose = verbose, name = 'basic tricky indexing', newUpdateNodes = list(covMat = matrix(c(0.7, .25, .25, .7), 2))) @@ -343,7 +345,7 @@ code <- nimbleCode({ sigma ~ dinvgamma(shape = a, rate = b) a ~ dexp(scale = 1.5) b ~ dexp(2.1) - mu0 ~ dnorm(0, .00001) # dflat() # dflat not handled + mu0 ~ dnorm(0, .00001) # dflat() # dflat not handled }) n <- 10 log_mu_init <- rnorm(1) @@ -364,7 +366,7 @@ test_ADModelCalculate(model, relTol = relTolTmp, verbose = verbose, name = 'stoc -## complicated indexing +## complicated indexing set.seed(1) code <- nimbleCode({ x[2:4,3:5] <- S[1:3,1:3] @@ -454,7 +456,7 @@ code <- nimbleCode({ rho ~ dgamma(2, 3) }) -set.seed(1) +set.seed(1) n <- 5 locs <- runif(n) dd <- fields::rdist(locs) @@ -475,8 +477,8 @@ relTolTmp[5] <- 1e-13 ## 350 sec. test_ADModelCalculate(model, useParamTransform = TRUE, useFasterRderivs = TRUE, - newUpdateNodes = list(dist = newDist, pr = newPr), checkCompiledValuesIdentical = FALSE, - relTol = relTolTmp, absTolThreshold = 1e-12, verbose = verbose, + newUpdateNodes = list(dist = newDist, pr = newPr), checkCompiledValuesIdentical = FALSE, + relTol = relTolTmp, absTolThreshold = 1e-12, verbose = verbose, name = 'dnorm with user-defined fxn for covariance with loops') ## other dmnorm parameterizations @@ -542,7 +544,7 @@ rGPdist <- nimbleFunction( assign('dGPdist', dGPdist, envir = .GlobalEnv) assign('rGPdist', rGPdist, envir = .GlobalEnv) -code <- nimbleCode({ +code <- nimbleCode({ y[1:n] ~ dGPdist(dist[1:n, 1:n], rho) rho ~ dgamma(2, 3) }) @@ -636,11 +638,11 @@ test_ADModelCalculate(model, useParamTransform = TRUE, relTol = relTol, verbose set.seed(1) code <- nimbleCode({ a ~ dgamma(1.1, 0.8) - for(i in 1:4) + for(i in 1:4) b[i] ~ dnorm(z[i], var = a) z[1] ~ dnorm(0, 1) z[2] ~ dnorm(0, 1) - z[3:4] ~ dmnorm(mu0[1:2], pr[1:2,1:2]) + z[3:4] ~ dmnorm(mu0[1:2], pr[1:2,1:2]) }) model <- nimbleModel(code, data = list(b = rnorm(4)), inits = list(a = 1.3, z = runif(4), pr = diag(2), mu0 = rep(0, 2))) ## calcNodes excludes det intermediates @@ -720,7 +722,7 @@ relTolTmp[5] <- 1e-9 ## 335 sec. test_ADModelCalculate(model, x = 'prior', useParamTransform = TRUE, newUpdateNodes = list(p = newP), newConstantNodes = list(y = newY), checkDoubleUncHessian = FALSE, relTol = relTolTmp, absTolThreshold = 1e-12, checkCompiledValuesIdentical = FALSE, - useFasterRderivs = TRUE, verbose = verbose, name = 'Dirichlet paramTransform') + useFasterRderivs = TRUE, verbose = verbose, name = 'Dirichlet paramTransform') ## Various likelihood-level non-differentiable constructs/distributions @@ -754,7 +756,7 @@ model <- nimbleModel(code, data = list(y = 1.5), inits = list(mu = rnorm(1))) relTolTmp <- relTol -test_ADModelCalculate(model, relTol = relTolTmp, +test_ADModelCalculate(model, relTol = relTolTmp, verbose = verbose, name = 'dflat') ## dcat as likelihood @@ -765,7 +767,7 @@ code <- nimbleCode({ p[1:3] ~ ddirch(alpha[1:3]) }) -model <- nimbleModel(code, data = list(z = 2), inits = list(alpha = 1:3, p = c(.2,.3,.5))) +model <- nimbleModel(code, data = list(z = 2), inits = list(alpha = 1:3, p = c(.2,.3,.5))) relTolTmp <- relTol relTolTmp <- relTol @@ -793,7 +795,7 @@ code <- nimbleCode({ p[1:3] ~ ddirch(alpha[1:3]) }) -model <- nimbleModel(code, data = list(y = rnorm(5)), inits = list(z = c(1,2,2,1,3), mu = rnorm(3), alpha = 1:3, p = c(.2,.3,.5))) +model <- nimbleModel(code, data = list(y = rnorm(5)), inits = list(z = c(1,2,2,1,3), mu = rnorm(3), alpha = 1:3, p = c(.2,.3,.5))) relTolTmp <- relTol relTolTmp[1] <- 1e-14 @@ -801,7 +803,7 @@ relTolTmp[3] <- 1e-5 relTolTmp[4] <- 1e-3 test_ADModelCalculate(model, wrt = c('mu','p'), relTol = relTolTmp, absTolThreshold = 1e-12, - newUpdateNodes = list(z=c(3,3,2,1,2)), xNew = list(p = c(.35,.25,.4)), + newUpdateNodes = list(z=c(3,3,2,1,2)), xNew = list(p = c(.35,.25,.4)), verbose = verbose, name = 'dcat latent, stochastic indexing', checkCompiledValuesIdentical = FALSE, useFasterRderivs = TRUE, useParamTransform = TRUE) diff --git a/packages/nimble/tests/testthat/test-parameterTransform.R b/packages/nimble/tests/testthat/test-parameterTransform.R index 63e29b3a2..6df961947 100644 --- a/packages/nimble/tests/testthat/test-parameterTransform.R +++ b/packages/nimble/tests/testthat/test-parameterTransform.R @@ -3,7 +3,10 @@ source(system.file(file.path('tests', 'testthat', 'test_utils.R'), package = 'ni RwarnLevel <- options('warn')$warn options(warn = 1) nimbleVerboseSetting <- nimbleOptions('verbose') +buildModelDerivsSetting = nimbleOptions('buildModelDerivs') nimbleOptions(verbose = FALSE) +nimbleOptions(buildModelDerivs = TRUE) + ## @@ -277,5 +280,4 @@ test_that('parameterTransform works for zero nodes', { options(warn = RwarnLevel) nimbleOptions(verbose = nimbleVerboseSetting) - - +nimbleOptions(buildModelDerivs = buildModelDerivsSetting)