R/RcppExports.R

In magi: MAnifold-Constrained Gaussian Process Inference

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

solveMagi <- function(yFull, odeModel, tvecFull, sigmaExogenous, phiExogenous, xInitExogenous, thetaInitExogenous, muExogenous, dotmuExogenous, priorTemperatureLevel = 1, priorTemperatureDeriv = 1, priorTemperatureObs = 1, kernel = "generalMatern", nstepsHmc = 500L, burninRatioHmc = 0.5, niterHmc = 10000L, stepSizeFactorHmc, nEpoch = 10L, bandSize = 20L, useFrequencyBasedPrior = FALSE, useBand = TRUE, useMean = TRUE, useScalerSigma = FALSE, useFixedSigma = FALSE, skipMissingComponentOptimization = FALSE, positiveSystem = FALSE, verbose = FALSE) {
    .Call('_magi_solveMagi', PACKAGE = 'magi', yFull, odeModel, tvecFull, sigmaExogenous, phiExogenous, xInitExogenous, thetaInitExogenous, muExogenous, dotmuExogenous, priorTemperatureLevel, priorTemperatureDeriv, priorTemperatureObs, kernel, nstepsHmc, burninRatioHmc, niterHmc, stepSizeFactorHmc, nEpoch, bandSize, useFrequencyBasedPrior, useBand, useMean, useScalerSigma, useFixedSigma, skipMissingComponentOptimization, positiveSystem, verbose)
}

#' sample from GP ODE for latent x, theta, sigma, and phi
#' the phi can be sampled together. However this incurs huge computational cost because the matrix inverse whenever phi2 changed
#' not used in the final method because of the computational cost.
#' @noRd
xthetaphisigmaSample <- function(xInitial, thetaInitial, phiInitial, sigmaInitial, yobs, xtimes, step, modelName = "FN", nsteps = 1L, traj = FALSE) {
    .Call('_magi_xthetaphisigmaSample', PACKAGE = 'magi', xInitial, thetaInitial, phiInitial, sigmaInitial, yobs, xtimes, step, modelName, nsteps, traj)
}

#' R wrapper for xthetaphisigmallik
#' Calculate the log posterior of x, theta, sigma, and phi
#' @param xlatent the matrix of latent ODE curve
#' @param theta the parameter of ODE function
#' @param phi the hyper-parameter of GP kernel
#' @param sigma the observation noise for each component of y
#' @param yobs matrix of observations
#' @param xtimes the time index of discretizations
#' @param modelName string of model name
#' export removed: function specific to pre-coded models only
#' @noRd
xthetaphisigmallikRcpp <- function(xlatent, theta, phi, sigma, yobs, xtimes, modelName = "FN") {
    .Call('_magi_xthetaphisigmallikRcpp', PACKAGE = 'magi', xlatent, theta, phi, sigma, yobs, xtimes, modelName)
}

#' R wrapper for basic_hmcC for normal distribution
#' @noRd
hmcNormal <- function(initial, step, lb, ub, nsteps = 1L, traj = FALSE) {
    .Call('_magi_hmcNormal', PACKAGE = 'magi', initial, step, lb, ub, nsteps, traj)
}

#' R wrapper for xthetallik
#' @noRd
xthetallik_rescaledC <- function(yobs, covVr, covRr, sigma, initial) {
    .Call('_magi_xthetallik_rescaledC', PACKAGE = 'magi', yobs, covVr, covRr, sigma, initial)
}

#' R wrapper for xthetallikBandApprox
#' @noRd
xthetallikBandApproxC <- function(yobs, covVr, covRr, sigma, initial) {
    .Call('_magi_xthetallikBandApproxC', PACKAGE = 'magi', yobs, covVr, covRr, sigma, initial)
}

#' R wrapper for xthetallik
#' @noRd
xthetallik_withmuC <- function(yobs, covVr, covRr, sigma, initial) {
    .Call('_magi_xthetallik_withmuC', PACKAGE = 'magi', yobs, covVr, covRr, sigma, initial)
}

#' R wrapper for xthetallik
#' @noRd
speedbenchmarkXthetallik <- function(yobs, covVr, covRr, sigmaScalar, initial, nrep = 10000L) {
    .Call('_magi_speedbenchmarkXthetallik', PACKAGE = 'magi', yobs, covVr, covRr, sigmaScalar, initial, nrep)
}

#' test for pass by reference for gpcov
#' @noRd
changeGPcovFromC <- function(covVr) {
    .Call('_magi_changeGPcovFromC', PACKAGE = 'magi', covVr)
}

cov_r2cpp_t1 <- function(cov_r) {
    invisible(.Call('_magi_cov_r2cpp_t1', PACKAGE = 'magi', cov_r))
}

cov_r2cpp_t2 <- function(cov_r) {
    invisible(.Call('_magi_cov_r2cpp_t2', PACKAGE = 'magi', cov_r))
}

cov_r2cpp_t3 <- function(cov_r) {
    invisible(.Call('_magi_cov_r2cpp_t3', PACKAGE = 'magi', cov_r))
}

#' R wrapper for phisigllik
#' @noRd
phisigllikHard2DC <- function(phisig, yobs, dist, kernel = "matern") {
    .Call('_magi_phisigllikHard2DC', PACKAGE = 'magi', phisig, yobs, dist, kernel)
}

#' R wrapper for xthetallik
#' @noRd
xthetallikC <- function(yobs, covVr, covRr, sigmaInput, initial, useBand = FALSE, priorTemperatureInput = 1.0) {
    .Call('_magi_xthetallikC', PACKAGE = 'magi', yobs, covVr, covRr, sigmaInput, initial, useBand, priorTemperatureInput)
}

#' test for output custom object
#' @noRd
maternCovTestOutput <- function(phi, dist, complexity = 0L) {
    .Call('_magi_maternCovTestOutput', PACKAGE = 'magi', phi, dist, complexity)
}

#' test for input custom object
#' @noRd
maternCovTestInput <- function(cov_input) {
    .Call('_magi_maternCovTestInput', PACKAGE = 'magi', cov_input)
}

#' R wrapper for phisigllik
#' phi sigma marginal log likelihood in GP smoothing.
#' used to set hyper-parameters phi
#' @param phisig vector of phi and sigma
#' @param yobs matrix of observations. if set each dimension separately, then supply a matrix of only 1 column
#' @param dist distance matrix on time index
#' @param kernel the type of kernel, support "matern", "rbf", "compact1", "periodicMatern", "generalMatern"
#' @noRd
phisigllikC <- function(phisig, yobs, dist, kernel = "matern") {
    .Call('_magi_phisigllikC', PACKAGE = 'magi', phisig, yobs, dist, kernel)
}

#' sample from GP marginal likelihood for phi and sigma
#' Internal function for debugging purpose
#' @noRd
phisigSample <- function(yobs, dist, initial, step, nsteps = 1L, traj = FALSE, kernel = "matern") {
    .Call('_magi_phisigSample', PACKAGE = 'magi', yobs, dist, initial, step, nsteps, traj, kernel)
}

#' general matern cov calculation Rcpp wrapper
#' @param phi the GP kernel hyper-parameters
#' @param distSigned signed distance of time indices
#' @param complexity the complexity of GP kernel calculation. Refer to documentation of `calCov`
#' @noRd
generalMaternCovRcpp <- function(phi, distSigned, complexity = 3L) {
    .Call('_magi_generalMaternCovRcpp', PACKAGE = 'magi', phi, distSigned, complexity)
}

#' sample from GP ODE for latent x and theta using HMC
#' Internal function for debugging purpose
#' @noRd
xthetaSample <- function(yobs, covList, sigmaInput, initial, step, nsteps = 1L, traj = FALSE, loglikflag = "usual", overallTemperature = 1, priorTemperatureInput = 1.0, modelName = "FN") {
    .Call('_magi_xthetaSample', PACKAGE = 'magi', yobs, covList, sigmaInput, initial, step, nsteps, traj, loglikflag, overallTemperature, priorTemperatureInput, modelName)
}

#' R wrapper for xthetallik
#' Calculate the log posterior of x and theta, without sigma
#' @param yobs matrix of observations
#' @param covAllDimInput list of covariance kernel objects
#' @param sigmaInput vector of sigma for each component of y
#' @param initial vector of x and theta, i.e., vectorise x from matrix, and then concatenate the theta vector
#' @param modelName string of model name
#' @param useBand boolean variable to use band matrix approximation
#' @param priorTemperatureInput the prior temperature for derivative, level, and observation, in that order
#' export removed: function specific to pre-coded models only
#' @noRd
xthetallikRcpp <- function(yobs, covAllDimInput, sigmaInput, initial, modelName = "FN", useBand = FALSE, priorTemperatureInput = 1.0) {
    .Call('_magi_xthetallikRcpp', PACKAGE = 'magi', yobs, covAllDimInput, sigmaInput, initial, modelName, useBand, priorTemperatureInput)
}

#' R wrapper for xthetallik
#' Internal function for debugging purpose
#' @noRd
xthetallikWithmuBandC <- function(yobs, covVr, covRr, sigmaInput, initial, useBand = TRUE, priorTemperatureInput = 1.0) {
    .Call('_magi_xthetallikWithmuBandC', PACKAGE = 'magi', yobs, covVr, covRr, sigmaInput, initial, useBand, priorTemperatureInput)
}

#' R wrapper for basic_hmcC
#' Internal function for debugging purpose
#' @noRd
basic_hmcRcpp <- function(rlpr, initial, step, lb, ub, nsteps = 1L, traj = FALSE) {
    .Call('_magi_basic_hmcRcpp', PACKAGE = 'magi', rlpr, initial, step, lb, ub, nsteps, traj)
}

#' R wrapper for chainSamplerRcpp
#' Internal function for debugging purpose
#' @noRd
chainSamplerRcpp <- function(yobs, covAllDimInput, nstepsInput, loglikflagInput, priorTemperatureInput, sigmaSizeInput, modelInput, niterInput, burninRatioInput, xthetasigmaInit, stepLowInit, positiveSystem, verbose) {
    .Call('_magi_chainSamplerRcpp', PACKAGE = 'magi', yobs, covAllDimInput, nstepsInput, loglikflagInput, priorTemperatureInput, sigmaSizeInput, modelInput, niterInput, burninRatioInput, xthetasigmaInit, stepLowInit, positiveSystem, verbose)
}

#' R wrapper for optimizeThetaInit
#' Internal function for debugging purpose
#' @noRd
optimizeThetaInitRcpp <- function(yobs, odeModel, covAllDimInput, sigmaAllDimensionsInput, priorTemperatureInput, xInitInput, useBandInput) {
    .Call('_magi_optimizeThetaInitRcpp', PACKAGE = 'magi', yobs, odeModel, covAllDimInput, sigmaAllDimensionsInput, priorTemperatureInput, xInitInput, useBandInput)
}

#' R wrapper for xthetasigmallik
#' Calculate the log posterior of x, theta, and sigma
#' @param xlatent the matrix of latent ODE curve
#' @param theta the parameter of ODE function
#' @param sigma the observation noise for each component of y
#' @param yobs matrix of observations
#' @param covAllDimInput list of covariance kernel objects
#' @param priorTemperatureInput the prior temperature for derivative, level, and observation, in that order
#' @param useBand boolean variable indicator to use band matrix approximation
#' @param useMean boolean variable indicator to use mean function in GP
#' @param modelName string of model name
#' export removed: function specific to pre-coded models only
#' @noRd
xthetasigmallikRcpp <- function(xlatent, theta, sigma, yobs, covAllDimInput, priorTemperatureInput = 1.0, useBand = FALSE, useMean = FALSE, modelName = "FN") {
    .Call('_magi_xthetasigmallikRcpp', PACKAGE = 'magi', xlatent, theta, sigma, yobs, covAllDimInput, priorTemperatureInput, useBand, useMean, modelName)
}

#' MAGI posterior density
#' 
#' Computes the MAGI log-posterior value and gradient for an ODE model with the given inputs: the observations \eqn{Y}, the latent system trajectories \eqn{X},
#' the parameters \eqn{\theta}, the noise standard deviations \eqn{\sigma}, and covariance kernels.
#' 
#' @param y data matrix of observations
#' @param xlatent matrix of system trajectory values
#' @param theta vector of parameter values \eqn{\theta}
#' @param sigma vector of observation noise for each system component
#' @param covAllDimInput list of covariance kernel objects for each system component. Covariance calculations may be carried out with \code{\link{calCov}}.
#' @param odeModel list of ODE functions and inputs. See details.
#' @param priorTemperatureInput  vector of tempering factors for the GP prior, derivatives, and observations, in that order. Controls the influence of the GP prior relative to the likelihood.  Recommended values: the total number of observations divided by the total number of discretization points for the GP prior and derivatives, and 1 for the observations.
#' @param useBand logical: should the band matrix approximation be used?  If \code{TRUE}, \code{covAllDimInput} must include CinvBand, mphiBand, and KinvBand as computed by \code{\link{calCov}}.
#' 
#' @return A list with elements \code{value} for the value of the log-posterior density and \code{grad} for its gradient.
#' 
#' @examples
#' # Trajectories from the Fitzhugh-Nagumo equations
#' tvec <- seq(0, 20, 2)
#' Vtrue <- c(-1, 1.91, 1.38, -1.32, -1.5, 1.73, 1.66, 0.89, -1.82, -0.93, 1.89)
#' Rtrue <- c(1, 0.33, -0.62, -0.82, 0.5, 0.94, -0.22, -0.9, -0.08, 0.95, 0.3)
#' 
#' # Noisy observations
#' Vobs <- Vtrue + rnorm(length(tvec), sd = 0.05)
#' Robs <- Rtrue + rnorm(length(tvec), sd = 0.1)
#' 
#' # Prepare distance matrix for covariance kernel calculation
#' foo <- outer(tvec, t(tvec), '-')[, 1, ]
#' r <- abs(foo)
#' r2 <- r^2
#' signr <- -sign(foo)
#'   
#' # Choose some hyperparameter values to illustrate
#' rphi <- c(0.95, 3.27)
#' vphi <- c(1.98, 1.12)
#' phiTest <- cbind(vphi, rphi)
#' 
#' # Covariance computations
#' curCovV <- calCov(phiTest[,1], r, signr, kerneltype = "generalMatern")
#' curCovR <- calCov(phiTest[,2], r, signr, kerneltype = "generalMatern")
#' 
#' # Y and X inputs to MagiPosterior
#' yInput <- data.matrix(cbind(Vobs, Robs))
#' xlatentTest <- data.matrix(cbind(Vtrue, Rtrue))
#' 
#' # Create odeModel list for FN equations
#' fnmodel <- list(
#'   fOde = fnmodelODE,
#'   fOdeDx = fnmodelDx,
#'   fOdeDtheta = fnmodelDtheta,
#'   thetaLowerBound = c(0, 0, 0),
#'   thetaUpperBound = c(Inf, Inf, Inf)
#' )
#' 
#' MagiPosterior(yInput, xlatentTest, theta = c(0.2, 0.2, 3), sigma = c(0.05, 0.1),
#'     list(curCovV, curCovR), fnmodel)
#' 
#'
#' @export
MagiPosterior <- function(y, xlatent, theta, sigma, covAllDimInput, odeModel, priorTemperatureInput = 1.0, useBand = FALSE) {
    .Call('_magi_MagiPosterior', PACKAGE = 'magi', y, xlatent, theta, sigma, covAllDimInput, odeModel, priorTemperatureInput, useBand)
}

#' sample from GP ODE for latent x, theta, and sigma
#' Internal function for debugging purpose
#' @noRd
xthetasigmaSample <- function(yobs, covList, sigmaInit, xthetaInit, step, nsteps = 1L, traj = FALSE, loglikflag = "usual", priorTemperatureInput = 1.0, modelName = "FN") {
    .Call('_magi_xthetasigmaSample', PACKAGE = 'magi', yobs, covList, sigmaInit, xthetaInit, step, nsteps, traj, loglikflag, priorTemperatureInput, modelName)
}

calcFrequencyBasedPrior <- function(x) {
    .Call('_magi_calcFrequencyBasedPrior', PACKAGE = 'magi', x)
}

gpsmooth <- function(yobsInput, distInput, kernelInput, sigmaExogenScalar = -1, useFrequencyBasedPrior = FALSE) {
    .Call('_magi_gpsmooth', PACKAGE = 'magi', yobsInput, distInput, kernelInput, sigmaExogenScalar, useFrequencyBasedPrior)
}

calcMeanCurve <- function(xInput, yInput, xOutput, phiCandidates, sigmaCandidates, kerneltype = "generalMatern", useDeriv = FALSE) {
    .Call('_magi_calcMeanCurve', PACKAGE = 'magi', xInput, yInput, xOutput, phiCandidates, sigmaCandidates, kerneltype, useDeriv)
}

optimizeThetaInit <- function(yobsInput, fOdeModelInput, covAllDimensionsInput, sigmaAllDimensionsInput, priorTemperatureInput, xInitInput, useBandInput) {
    .Call('_magi_optimizeThetaInit', PACKAGE = 'magi', yobsInput, fOdeModelInput, covAllDimensionsInput, sigmaAllDimensionsInput, priorTemperatureInput, xInitInput, useBandInput)
}

optimizePhi <- function(yobsInput, tvecInput, fOdeModelInput, sigmaAllDimensionsInput, priorTemperatureInput, xInitInput, thetaInitInput, phiInitInput, missingComponentDim) {
    .Call('_magi_optimizePhi', PACKAGE = 'magi', yobsInput, tvecInput, fOdeModelInput, sigmaAllDimensionsInput, priorTemperatureInput, xInitInput, thetaInitInput, phiInitInput, missingComponentDim)
}

solveMagiRcpp <- function(yFull, odeModel, tvecFull, sigmaExogenous, phiExogenous, xInitExogenous, thetaInitExogenous, muExogenous, dotmuExogenous, priorTemperatureLevel, priorTemperatureDeriv, priorTemperatureObs, kernel, nstepsHmc, burninRatioHmc, niterHmc, stepSizeFactorHmc, nEpoch, bandSize, useFrequencyBasedPrior, useBand, useMean, useScalerSigma, useFixedSigma, skipMissingComponentOptimization, positiveSystem, verbose) {
    .Call('_magi_solveMagiRcpp', PACKAGE = 'magi', yFull, odeModel, tvecFull, sigmaExogenous, phiExogenous, xInitExogenous, thetaInitExogenous, muExogenous, dotmuExogenous, priorTemperatureLevel, priorTemperatureDeriv, priorTemperatureObs, kernel, nstepsHmc, burninRatioHmc, niterHmc, stepSizeFactorHmc, nEpoch, bandSize, useFrequencyBasedPrior, useBand, useMean, useScalerSigma, useFixedSigma, skipMissingComponentOptimization, positiveSystem, verbose)
}

#' basic_hmcC
#' 
#' BASIC HAMILTONIAN MONTE CARLO UPDATE
#' Shihao Yang, 2017, rewriten from Radford M. Neal, 2012.
#'
#' @param lpr       Function returning the log probability of the position part 
#'                  of the state, plus an arbitrary constant, with gradient
#'                  as an attribute if grad=TRUE is passed.
#' @param initial   The initial position part of the state (a vector).
#' @param nsteps    Number of steps in trajectory used to propose a new state.
#'                  (Default is 1, giving the "Langevin" method.)
#' @param step      Stepsize or stepsizes.  May be scalar or a vector of length 
#'                  equal to the dimensionality of the state.
#' @param traj      TRUE if values of q and p along the trajectory should be 
#'                  returned (default is FALSE).
#' @noRd
NULL

#' log likelihood for latent states and ODE theta conditional on phi sigma
#'   latent states have a mean mu
#' @param phisig      the parameter phi and sigma
#' @param yobs        observed data
NULL

#' log likelihood for latent states and ODE theta conditional on phi sigma
#'
#' the scale is in fact taken out and it is the legacy version of xthetallik
#'
#' @param phisig      the parameter phi and sigma
#' @param yobs        observed data
NULL

#' approximate log likelihood for latent states and ODE theta conditional on phi sigma
#'
#' band matrix approximation
#'
#' @param phisig      the parameter phi and sigma
#' @param yobs        observed data
#' @noRd
#' FIXME xtheta currently passed by value for Fortran code
NULL

#' log likelihood for Gaussian Process marginal likelihood with Matern kernel
#'
#' @param phisig      the parameter phi and sigma
#' @param yobs        observed data
NULL

hmcTest <- function() {
    .Call('_magi_hmcTest', PACKAGE = 'magi')
}

bandTest <- function(filename = "data_band.txt") {
    .Call('_magi_bandTest', PACKAGE = 'magi', filename)
}

#' matern variance covariance matrix with derivatives
#' 
#' @param phi         the parameter of (sigma_c_sq, alpha)
#' @param dist        distance matrix
#' @param complexity  how much derivative information should be calculated
NULL

#' matern variance covariance matrix with derivatives
#' 
#' @param phi         the parameter of (sigma_c_sq, alpha)
#' @param dist        distance matrix
#' @param complexity  how much derivative information should be calculated
NULL

#' periodic matern variance covariance matrix with derivatives
#' 
#' @param phi         the parameter of (sigma_c_sq, alpha)
#' @param dist        distance matrix
#' @param complexity  how much derivative information should be calculated
NULL

#' log likelihood for Gaussian Process marginal likelihood with Matern kernel
#' 
#' @param phisig      the parameter phi and sigma
#' @param yobs        observed data
NULL

#' log likelihood for latent states and ODE theta conditional on phi sigma
#' 
#' @param phisig      the parameter phi and sigma
#' @param yobs        observed data
NULL