cobyla call behaves differently between nlopt 2.7.1 and 2.10.0 #195
Description
The reprex is:
library(nloptr)
internal_whittle_garma_obj_short <- function(par, params) {
# Just the first part of the objective function.
# This is the "S" function of Giraitis, Hidalgo & Robinson 2001.
ss <- params$ss
spec <- ss$spec
n_freq <- length(spec)
spec_den_inv <- internal_calc_garma_spec_den_inv(par, params)
## I have checked the source for for "spectrum" which is used to calculate the "spec"
## The I/f calc includes spectrum and the "f" we calc. f includes the factor 1/(2 pi)
## Source code for "spectrum" however divides by n but not by 2*pi. Thus we make allowances for that here.
I_f <- spec * spec_den_inv / (2 * pi * n_freq)
res <- sum(I_f, na.rm = TRUE)
return(res)
}
internal_calc_garma_spec_den_inv <- function(par, params) {
# This is the "k" function of Giraitis, Hidalgo & Robinson 2001.
# true spec den is this times sigma^2/(2*pi)
ss <- params$ss
p <- params$p
q <- params$q
k <- params$k
n_freq <- length(ss$freq)
freq <- ss$freq[1:n_freq]
spec <- ss$spec[1:n_freq]
start <- 1
if (is.null(params$periods)) {
u <- fd <- numeric(0)
for (k1 in seq_len(k)) {
u <- c(u, par[start])
fd <- c(fd, par[start + 1])
start <- start + 2
}
} else {
k <- length(params$periods)
u <- cos(2 * pi / params$periods)
fd <- par[start:(start + k - 1)]
start <- start + k
}
if (p > 0) phi_vec <- (-par[start:(start + p - 1)]) else phi_vec <- 1
if (q > 0) theta_vec <- (-par[(p + start):length(par)]) else theta_vec <- 1
cos_2_pi_f <- cos(2 * pi * freq)
mod_phi <- mod_theta <- 1
if (p > 0) mod_phi <- internal_a_fcn(phi_vec, freq)
if (q > 0) mod_theta <- internal_a_fcn(theta_vec, freq)
## As per Giraitis, Hidalgo & Robinson 2001, we do not include sigma^2/2pi factor in the spectral density.
## The below is the inverse of the spectral density.
spec_den_inv <- mod_phi / mod_theta
for (k1 in seq_len(k)) {
u_k <- u[k1]
fd_k <- fd[k1]
spec_den_inv <- spec_den_inv * (4 * ((cos_2_pi_f - u_k)^2))^fd_k
}
return(spec_den_inv)
}
params <- readRDS("params.rds")
initial_pars <- c(0.24868989, 0.05572044)
lb <- c(1e-15, 0e+00)
ub <- c(1.0, 0.5)
nlopt_control <- list(
maxeval = 10000,
ftol_rel = 1e-8,
xtol_rel = 0
)
x <- nloptr::cobyla(
x0 = initial_pars, fn = internal_whittle_garma_obj_short, lower = lb, upper = ub,
params = params, control = nlopt_control
)
x$parWhen I install nloptr using:
install_github("astamm/nloptr")i.e. current master with nlopt 2.10.0 sources, I get:
> x$par
[1] NaN NaN
When I install nloptr using:
install_github("astamm/nloptr@04b973a")i.e. commit right before updating sources to2.10.0 and therefore using 2.7.1, I get:
> x$par
[1] 0.24868989 0.05572042
Needs investigation to see if this is nlopt or nloptr issue.