testDynamicalModel {magi}

Test dynamic system model specification

Description

Given functions for the ODE and its gradients (with respect to the system components and parameters), verify the correctness of the gradients using numerical differentiation.

Usage

testDynamicalModel(modelODE, modelDx, modelDtheta, modelName, x, theta, tvec)

Arguments

modelODE

function that computes the ODEs, specified with the form f(theta, x, tvec). See examples.

modelDx

function that computes the gradients of the ODEs with respect to the system components. See examples.

modelDtheta

function that computes the gradients of the ODEs with respect to the parameters \theta. See examples.

modelName

string giving a name for the model

x

data matrix of system values, one column for each component, at which to test the gradients

theta

vector of parameter values for \theta, at which to test the gradients

tvec

vector of time points corresponding to the rows of x

Details

Calls test_that to test equality of the analytic and numeric gradients.

Value

A list with elements testDx and testDtheta, each with value TRUE if the corresponding gradient check passed and FALSE if not.

Examples

# ODE system and gradients for Fitzhugh-Nagumo equations: fnmodelODE, fnmodelDx, fnmodelDtheta

# Example of incorrect gradient with respect to parameters theta
fnmodelDthetaWrong <- function(theta, x, tvec) {
  resultDtheta <- array(0, c(nrow(x), length(theta), ncol(x)))
  
  V = x[, 1]
  R = x[, 2]
  
  resultDtheta[, 3, 1] = V - V^3 / 3.0 - R
  
  resultDtheta[, 1, 2] =  1.0 / theta[3] 
  resultDtheta[, 2, 2] = -R / theta[3]
  resultDtheta[, 3, 2] = 1.0 / (theta[3]^2) * (V - theta[1] + theta[2] * R)
  
  resultDtheta
}

# Sample data for testing gradient correctness
data(FNdat)
 
# Correct gradients
testDynamicalModel(fnmodelODE, fnmodelDx, fnmodelDtheta, 
    "FN equations", FNdat[, c("V", "R")], c(.5, .6, 2), FNdat$time)
    
# Incorrect theta gradient (test fails)
testDynamicalModel(fnmodelODE, fnmodelDx, fnmodelDthetaWrong, 
    "FN equations", FNdat[, c("V", "R")], c(.5, .6, 2), FNdat$time)