#### Original in ~/R/Pkgs/DistributionUtils/inst/unitTests/runit.incompleteBesselK.R

### "Translated" from RUnit unit tests to simple base R tests by Martin Maechler
##  checkEquals(....)              ==> stopifnot(all.equal(<....>))
##  checkTrue  (...., msg = <msg>) ==> stopifnot(<msg> = <....>)

### Unit tests of function incompleteBesselK and incompleteBesselKR

## Author: David Scott, Date: 5 Jan 2012

## Original version of the test file which was misplaced

library(Bessel)
## Data from Harris (2008)
HarrisCase1 <- c(2.225310761266469,
                 0.213894166822940,
                 0.054503469799701,
                 0.023253121507708,
                 0.013042750996080,
                 0.008567534990649,
                 0.006208676806601,
                 0.004801085238177,
                 0.003884072049627,
                 0.003246798003149)

HarrisCase2 <- 0.000012249987981

HarrisCase3 <- 0.00000041500106421228

HarrisCase4 <- 0.000528504325244

## Calculations for Harris Case 1
numIBF <- numeric(10)
for (i in 0:9){
    ibf <- incompleteBesselK(0.01, 4, i, nmax = 100)
    numIBF[i + 1] <- ibf
}
maxDiff <- max(abs(numIBF - HarrisCase1))
## stopifnot(paste("Harris Case 1: maxDiff =", maxDiff) = maxDiff <= 10^(-13))
if(!(maxDiff <= 10^(-13))) stop(paste("Harris Case 1: maxDiff =", maxDiff))

## Calculations for Harris Case 2
ibf <- incompleteBesselK(4.95, 5, 2)
maxDiff <- max(abs(ibf - HarrisCase2))
stopifnot("Harris Case 2" = maxDiff <= 10^(-13))

## Calculations for Harris Case 3
ibf <- incompleteBesselK(10, 2, 6)
maxDiff <- max(abs(ibf - HarrisCase3))
stopifnot("Harris Case 3" = maxDiff <= 10^(-11))

## Calculations for Harris Case 4
ibf <- incompleteBesselK(3.1, 2.6, 5)
maxDiff <- max(abs(ibf - HarrisCase4))
## stopifnot(paste("Harris Case 4: maxDiff =", maxDiff) = maxDiff <= 10^(-13))
if(!(maxDiff <= 10^(-13))) stop(paste("Harris Case 4: maxDiff =", maxDiff))

## Newer version of the test file: not necessarily better
## But includes testing of the pure R version

## Values given by Harris (2008)
## MM: should ask for more precision than default sqrt(.Machine$double.eps) = 1.4901e-8
stopifnot(all.equal(incompleteBesselK (0.01, 4, 0), 2.225310761266469))
stopifnot(all.equal(incompleteBesselKR(0.01, 4, 0), 2.225310761266469))

stopifnot(all.equal(incompleteBesselK (0.01, 4, 9), 0.003246798003149))
stopifnot(all.equal(incompleteBesselKR(0.01, 4, 9), 0.003246798003149))

stopifnot(all.equal(incompleteBesselK (4.95, 5, 2), 0.000012249987981))
stopifnot(all.equal(incompleteBesselKR(4.95, 5, 2), 0.000012249987981))

stopifnot(all.equal(incompleteBesselK (10, 2, 6), 0.0000004150045864731308))
stopifnot(all.equal(incompleteBesselKR(10, 2, 6), 0.0000004150045864731308))

stopifnot(all.equal(incompleteBesselK (3.1, 2.6, 5), 0.000528504325244))
stopifnot(all.equal(incompleteBesselKR(3.1, 2.6, 5), 0.000528504325244))

### Check values when x > y using numeric integration:
integrBess <- function(t, x, y, nu)  t^-(nu+1) * exp(-x*t - y/t)

nus <- seq(0, 12, by = 1/4)
x. <- 4
y. <- 0.01

numIBF <- sapply(nus, incompleteBesselK, x = x., y = y.)
integIBF <- sapply(nus, integrate, f = integrBess, lower = 1, upper = Inf, x = x., y = y.)
## --> a  5 x #{nus}  matrix of lists
intIBF <- as.numeric(integIBF["value", ])
errIBF <- as.numeric(integIBF["abs.error", ])
numIBF - intIBF # -1.257989e-11 -1.154437e-11 -9.600298e-12 -7.442554e-12 -5.422429e-12 .....
stopifnot(identical(length(nus), as.integer(table(unlist(integIBF["message",])))),
          max(abs(numIBF - intIBF)) < 1e-10,
          abs(numIBF - intIBF) < errIBF)

plot(nus, numIBF, type = "o", col = 2, cex = 1/2)