# R/diffusion.R In ctmm: Continuous-Time Movement Modeling

#### Defines functions diffusionDiff.OUO.fnDiff.OUF.fn

```# NECESSARY NOT-SO-SPECIAL FUNCTIONS

# the z->1 limit is numerically unstable
Diff.OUF.fn <- function(z,deriv=0)
{
w <- 1-z
if(w > 0.004)
{
if(deriv==0)
{ z <- z^(z/w) }
else if(deriv==1)
{ z <- Diff.OUF.fn(z) * ( w + log(z) )/w^2 }
}
else
{
if(deriv==0)
{ z <- exp(-1) * series(w,c(1,1/2,7/24,3/16,743/5760,215/2304)) }
else if(deriv==1)
{ z <- -Diff.OUF.fn(z) * series(w,1/2:10) }
}
return(z)
}

# NaN at 0
Diff.OUO.fn <- function(x,deriv=0)
{
if(deriv==0)
{ x <- sqrt(1+x^2) * exp(-ifelse(x==0,1,atan(x)/x)) }
else if(deriv==1)
{
if(x>0.1)
{ x <- Diff.OUO.fn(x)/x^2 * ( x - 2*x/(1+x^2) + atan(x) ) }
else
{
n <- 10
coef <- seq(5,by=4,length.out=n)/seq(3,by=2,length.out=n) * (-1)^(1+1:n)
x <- Diff.OUO.fn(x) * x*series(x^2,coef)
}
}
return(x)
}

# calculate the max_lag diffusion rate for a movement model
diffusion <- function(CTMM,level=0.95,finish=TRUE)
{
CTMM <- get.taus(CTMM) # pre-calculate stuff

range <- CTMM\$range
tau <- CTMM\$tau
omega <- CTMM\$omega
f <- CTMM\$f.nu[1]
nu <- CTMM\$f.nu[2]
Omega2 <- CTMM\$Omega2
circle <- abs(CTMM\$circle) # only magnitude matters for this
circle <- FALSE # turning circulation off, because it reduces diffusion rate

sigma <- var.covm(CTMM\$sigma,ave=FALSE) # variance
COV <- CTMM\$COV
if(!is.null(COV))
{
COV <- axes2var(CTMM,MEAN=FALSE)
}
else
# this is redundant with D.grad, but needed for mean.ctmm summary
J <- rep(0,length(CTMM\$features))
names(J) <- CTMM\$features

if(!length(tau) || all(tau==0)) # IID
{
return( c(0,Inf,Inf) )
}
else if(!range) # Brownian motion # IOU # max diffusion rate at infinite lag
{ D <- 1 }
else if(length(tau)==1 || tau[2]==0) # OU
{
tau <- tau[1]
NAMES <- CTMM\$tau.names[1]

if(circle*tau <= 1) # max diffusion rate at zero lag
{
D <- 1/tau
}
else # circulation enhanced max diffusion rate
{
NAMES <- c(NAMES,"circle")
z <- circle*tau
# D <- atan((z^2-1)/(2*z))/circle
D <- (2*atan(z)-pi/2)/circle
}
}
else if(length(tau)==2)
{
NAMES <- CTMM\$tau.names

if(!omega && tau[1]!=tau[2] && !circle) # OUF
{
z <- tau[2]/tau[1]
D <- Diff.OUF.fn(z) / tau[1]
}
else if(omega && tau[1]==tau[2] && !circle) # OUO
{
z <- tau[1]*omega
D <- Diff.OUO.fn(z) / tau[1]
}
else if(!omega && tau[1]==tau[2] && !circle) # OUf
{
NAMES <- NAMES[1]
D <- exp(-1)/tau[1]
}
}
# else if(circle) # OUF, OUO, OUf with circulation can't be solved analytically
# {
#   if(!omega && tau[1]!=tau[2]) # OUF
#   {
#     NAMES <- paste("tau",NAMES)
#     S0.fn <- function(t) { -diff(exp(-t/tau)*tau)/diff(tau) }
#     D0.fn <- function(t) { diff(exp(-t/tau))/diff(tau) }
#     D1.fn <- function(t) { -diff(exp(-t/tau)/tau)/diff(tau) }
#     J <- diag(2)
#     LAG0 <- log(tau[1]/tau[2])/(1/tau[2]-1/tau[1])
#   }
#   else if(!omega && tau[1]==tau[2]) # OUf
#   {
#     S0.fn <- function(t) { -(1+f*t) * exp(-f*t) }
#     D0.fn <- function(t) { f^2*t * exp(-f*t) }
#     D1.fn <- function(t) { f^2*(1-f*t) * exp(-f*t) }
#     J <- CTMM\$J.f.tau
#     LAG0 <- tau[1]
#   }
#   else if(omega) # OUO
#   {
#     S0.fn <- function(t) { -(cos(nu*t)+f/nu*sin(nu*t)) * exp(-f*t) }
#     D0.fn <- function(t) { Omega2/nu * sin(nu*t) * exp(-f*t) }
#     D1.fn <- function(t) { Omega2 * ( cos(nu*t) - f/nu*sin(nu*t) ) * exp(-f*t) }
#     J <- CTMM\$J.nu.tau
#     LAG0 <- atan(nu*t)/nu
#   }
#
#   # negative diffusion rate function of lag
#   nD.fn <- function(t) { -(D0.fn(t)*cos(circle*t) - S0.fn(t)*circle*sin(circle*t)) }
#
#   MAX <- optimizer(LAG0,nD.fn,lower=0)
#   LAG <- MAX\$par
#   D <- -MAX\$value
#
#   # UNFINISHED
#   # UNFINISHED
#   # UNFINISHED
#
#   # zero circulation result
#   R0 <- CTMM
#   R0\$circle <- FALSE
#   R0 <- diffusion(R0,finish=FALSE)
#
#   # very annoying to calculate this better (optimize inside gradient)
# }

J <- sigma * J

if(CTMM\$isotropic[1])
{ J["major"] <- 2*D }
else
{ J["major"] <- J["minor"] <- D }
D <- sigma * D

if(!is.null(COV))
{
VAR <- nant(VAR,Inf)
}
else
{ VAR <- Inf }

# this is for chi^2 CIs
DOF <- 2*D^2/VAR # / length(CTMM\$axes) # /2 is for counting

# return information for CIs but not completed CIs