Running 'bar.R' [0s]
  Comparing 'bar.Rout' to 'bar.Rout.save' ... OK
  Running 'baz.R' [0s]
  Comparing 'baz.Rout' to 'baz.Rout.save' ... OK
  Running 'foo.R' [0s]
  Comparing 'foo.Rout' to 'foo.Rout.save' ... OK
  Running 'fred.R' [159s]
  Running 'goo.R' [0s]
  Comparing 'goo.Rout' to 'goo.Rout.save' ... OK
  Running 'gred.R' [2s]
  Comparing 'gred.Rout' to 'gred.Rout.save' ... OK
  Running 'hero.R' [0s]
  Comparing 'hero.Rout' to 'hero.Rout.save' ... OK
  Running 'hoo.R' [0s]
  Comparing 'hoo.Rout' to 'hoo.Rout.save' ... OK
  Running 'poo.R' [155s]
  Running 'qux.R' [0s]
  Comparing 'qux.Rout' to 'qux.Rout.save' ... OK
Running the tests in 'tests/fred.R' failed.
Complete output:
  > 
  >  library(trust)
  > 
  >  options(digits = 3)
  > 
  >  ##### four-way contingency table with all two-way interactions
  > 
  >  d <- c(3, 4, 5, 6)
  >  n <- 1000
  > 
  >  ##### model matrix
  >  m <- NULL
  >  for (i in 1:d[1]) {
  +      for (j in 1:d[2]) {
  +          mfoo <- array(0, dim = d)
  +          mfoo[i, j, , ] <- 1
  +          mfoo <- as.vector(mfoo)
  +          m <- cbind(m, mfoo)
  +      }
  +  }
  >  for (i in 1:d[1]) {
  +      for (j in 1:d[3]) {
  +          mfoo <- array(0, dim = d)
  +          mfoo[i, , j, ] <- 1
  +          mfoo <- as.vector(mfoo)
  +          m <- cbind(m, mfoo)
  +      }
  +  }
  >  for (i in 1:d[1]) {
  +      for (j in 1:d[4]) {
  +          mfoo <- array(0, dim = d)
  +          mfoo[i, , , j] <- 1
  +          mfoo <- as.vector(mfoo)
  +          m <- cbind(m, mfoo)
  +      }
  +  }
  >  for (i in 1:d[2]) {
  +      for (j in 1:d[3]) {
  +          mfoo <- array(0, dim = d)
  +          mfoo[ , i, j, ] <- 1
  +          mfoo <- as.vector(mfoo)
  +          m <- cbind(m, mfoo)
  +      }
  +  }
  >  for (i in 1:d[2]) {
  +      for (j in 1:d[4]) {
  +          mfoo <- array(0, dim = d)
  +          mfoo[ , i, , j] <- 1
  +          mfoo <- as.vector(mfoo)
  +          m <- cbind(m, mfoo)
  +      }
  +  }
  >  for (i in 1:d[3]) {
  +      for (j in 1:d[4]) {
  +          mfoo <- array(0, dim = d)
  +          mfoo[ , , i, j] <- 1
  +          mfoo <- as.vector(mfoo)
  +          m <- cbind(m, mfoo)
  +      }
  +  }
  >  dimnames(m) <- NULL
  >  foo <- qr(m)
  >  m <- m[ , foo$pivot[seq(1, foo$rank)]]
  > 
  >  ##### true parameter value
  >  set.seed(42)
  >  theta.true <- 0.25 * rnorm(ncol(m))
  >  theta.true <- round(theta.true, 5)
  > 
  >  ##### simulate data
  >  eta <- as.numeric(m %*% theta.true)
  >  p <- exp(eta)
  >  p <- p / sum(p)
  >  x <- sample(nrow(m), n, replace = TRUE, prob = p)
  >  x <- tabulate(x, nbins = nrow(m))
  > 
  >  ##### save data
  >  iffy <- try(read.table("fred.txt"), silent = TRUE)
  >  if (inherits(iffy, "try-error")) {
  +      data <- data.frame(x = x, m = m)
  +      write.table(data, file = "fred.txt", row.names = FALSE)
  +  }
  >  data <- read.table(file = "fred.txt", header = TRUE)
  >  x <- data$x
  >  data$x <- NULL
  >  m <- as.matrix(data)
  >  dimnames(m) <- NULL
  > 
  >  ##### log likelihood
  >  objfun <- function(theta) {
  +      eta <- as.numeric(m %*% theta)
  +      p <- exp(eta)
  +      f <- sum(x * eta - p)
  +      g <- as.numeric(t(x - p) %*% m)
  +      B <- sweep(- m, 1, p, "*")
  +      B <- t(m) %*% B
  +      list(value = f, gradient = g, hessian = B)
  +  }
  > 
  >  ##### check it
  >  sally <- objfun(theta.true)
  >  epsilon <- 1e-8
  >  mygrad <- double(length(theta.true))
  >  for (i in 1:length(mygrad)) {
  +      theta.eps <- theta.true
  +      theta.eps[i] <- theta.true[i] + epsilon
  +      sally.eps <- objfun(theta.eps)
  +      mygrad[i] <- (sally.eps$value - sally$value) / epsilon
  +  }
  >  all.equal(sally$gradient, mygrad, tolerance = length(mygrad) * epsilon)
  [1] TRUE
  >  myhess <- matrix(NA, length(theta.true), length(theta.true))
  >  for (i in 1:length(mygrad)) {
  +      theta.eps <- theta.true
  +      theta.eps[i] <- theta.true[i] + epsilon
  +      sally.eps <- objfun(theta.eps)
  +      myhess[i, ] <- (sally.eps$gradient - sally$gradient) / epsilon
  +  }
  >  all.equal(sally$hessian, myhess, tolerance = length(mygrad) * epsilon)
  [1] TRUE
  > 
  >  fred <- trust(objfun, theta.true, 1, sqrt(ncol(m)), minimize = FALSE)
  > 
  >  fran <- glm.fit(m, x, family = poisson(), intercept = FALSE)
  > 
  >  all.equal(fran$coefficients, fred$argument)
  [1] TRUE
  > 
  >  fred <- trust(objfun, rep(0, length(theta.true)), 1, sqrt(ncol(m)),
  +      minimize = FALSE, blather = TRUE)
  >  fred$converged
  [1] TRUE
  >  #### CRAN don't like: gives different results on different hardware
  >  #### or different compiler flags
  >  ## ceiling(log10(max(abs(fred$gradient))))
  >  ## length(fred$r)
  >  ## data.frame(type = fred$steptype, rho = fred$rho, change = fred$preddiff,
  >  ##     accept = fred$accept, r = fred$r)
  >  ## (fred$stepnorm / fred$r)[fred$accept & fred$steptype != "Newton"]
  > 
  >  fred <- trust(objfun, rep(-5, length(theta.true)), 1, sqrt(ncol(m)),
  +      minimize = FALSE, blather = TRUE)
  >  fred$converged
  [1] TRUE
  >  #### CRAN don't like: gives different results on different hardware
  >  #### or different compiler flags
  >  ## ceiling(log10(max(abs(fred$gradient))))
  >  ## length(fred$r)
  >  ## data.frame(type = fred$steptype, rho = fred$rho, change = fred$preddiff,
  >  ##     accept = fred$accept, r = fred$r)
  >  ## (fred$stepnorm / fred$r)[fred$accept & fred$steptype != "Newton"]
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.89    0.28    1.15 
Running the tests in 'tests/poo.R' failed.
Complete output:
  > 
  >  objfun <- function(x) {
  +      ##### Rosenbrock's function #####
  +      stopifnot(is.numeric(x))
  +      stopifnot(length(x) == 2)
  +      f <- expression(100 * (x2 - x1^2)^2 + (1 - x1)^2)
  +      g1 <- D(f, "x1")
  +      g2 <- D(f, "x2")
  +      h11 <- D(g1, "x1")
  +      h12 <- D(g1, "x2")
  +      h22 <- D(g2, "x2")
  +      x1 <- x[1]
  +      x2 <- x[2]
  +      f <- eval(f)
  +      g <- c(eval(g1), eval(g2))
  +      B <- rbind(c(eval(h11), eval(h12)), c(eval(h12), eval(h22)))
  +      list(value = f, gradient = g, hessian = B)
  +  }
  > 
  >  library(trust)
  > 
  >  parinit <- c(3, 1)
  > 
  >  out <- trust(objfun, parinit, 1, 1e5, blather = TRUE)
  >  out$converged
  [1] TRUE
  >  length(out$r)
  [1] 21
  > 
  >  parscale <- c(5, 1)
  >  shift <- 4
  >  theta <- parscale * (parinit + shift)
  > 
  >  pobjfun <- function(x) {
  +      out <- objfun(x / parscale - shift)
  +      out$gradient <- out$gradient / parscale
  +      out$hessian <- out$hessian / outer(parscale, parscale)
  +      return(out)
  +  }
  > 
  >  pout <- trust(pobjfun, theta, 1, 1e5, blather = TRUE)
  >  pout$converged
  [1] TRUE
  >  length(pout$r)
  [1] 25
  >  all.equal(out$argument, pout$argument / parscale - shift)
  [1] "Mean relative difference: 1.011409e-07"
  > 
  >  qout <- trust(objfun, parinit, 1, 1e5, parscale = parscale, blather = TRUE)
  >  qout$converged
  [1] TRUE
  >  length(qout$r)
  [1] 25
  > 
  >  all.equal(pout$valpath, qout$valpath)
  [1] TRUE
  >  transpath <- pout$argpath
  >  transpath <- sweep(transpath, 2, parscale, "/")
  >  transpath <- sweep(transpath, 2, shift)
  >  all.equal(transpath, qout$argpath)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.20    0.06    0.25