R help - Oct 2018 - Maximum number of iterations (maxit) does not work in hydroPSO-modFit-optimr in r

The argument maxit used in libraries optimr, hydroPSO, and FME to
control the maximal number of iterations to be performed does not work
at ALL!!

This needs to be fixed...



mysamp <- function(n, m, s, lwr, upr, nnorm) {
  samp <- rnorm(nnorm, m, s)
  samp <- samp[samp >= lwr & samp <= upr]
  if (length(samp) >= n) {
    return(sample(samp, n))
  }
  stop(simpleError("Not enough values to sample from. Try increasing
nnorm."))
}

x <- mysamp(1000,0.8117648,0.1281839,0,1,1000000)
y <- mysamp(1000,0.7530346,0.1865946,0,1,1000000)

n <- length(x)
xgroup <- c('a','b')
ygroup <- c('c','d')

for (i in c(2:n)){
  xgroup[i] <-  c('b')
  ygroup[i] <-  c('d')
}

dataset <- data.frame(x = x, y = y, xgroup = as.factor(xgroup),
                      ygroup = as.factor(ygroup))
dataset$ObsNo <- 1:n
dataset <- dataset[order(dataset$x,decreasing=F),]

xopt <- c(dataset$x[1],0.95)    #Initial critical x,y values


my.function <- function(xopt){
  for (i in c(1:n)){
    dataset$xgroup[i] <- if(dataset$x[i] < xopt[1]) 'a' else
'b'
    dataset$ygroup[i] <- if(dataset$y[i] < xopt[2]) 'c' else
'd'
    dataset$q.i[i] <- with(dataset, ifelse
                           (dataset$xgroup[i]=='a' &
dataset$ygroup[i]=='d', 1, 0))
    dataset$q.ii[i] <- with(dataset, ifelse
                            (dataset$xgroup[i]=='b' &
dataset$ygroup[i]=='d', 1, 0))
    dataset$q.iii[i] <- with(dataset, ifelse
                             (dataset$xgroup[i]=='b' &
dataset$ygroup[i]=='c', 1, 0))
    dataset$q.iv[i] <- with(dataset, ifelse
                            (dataset$xgroup[i]=='a' &
dataset$ygroup[i]=='c', 1, 0))

    dataset$q.err[i] <- sum(dataset$q.i[i] + dataset$q.iii[i])
  }
  min.qerr <- sum(dataset$q.err)
  q.I <- sum(dataset$q.i)
  q.II <- sum(dataset$q.ii)
  q.III <- sum(dataset$q.iii)
  q.IV <- sum(dataset$q.iv)
  q.Iandq.III <- sum(dataset$q.err)
  print(c(q.I,
          q.II,
          q.III,
          q.IV,
          q.Iandq.III))
  return(min.qerr)
}

my.function(xopt)

#-------Algorithm

xmin=c(0,0.60)
xmax=c(0.95,0.95)

res= hydroPSO(par=xopt,my.function,method="spso2011",
              lower=xmin,upper=xmax,control=list(maxit=3))


res= modFit(my.function,p=xopt,method="Pseudo",
            lower=xmin,upper=xmax,control=list(numiter=3))

res= optimr(par=xopt,my.function,method="L-BFGS-B",
            lower=xmin,upper=xmax,control=list(maxit=3,trace=T))

#Only the hjkb form dfoptim that has argument maxfeval and it works.

res=hjkb(par=xopt,my.function,
         lower=xmin,upper=xmax,control=list(maxfeval=3,tol=2^-10,info=T))






Ahmed Attia, Ph.D.
Agronomist & Crop Modeler

As the author of optimx (there is a new version just up), I can say that
maxit only "works" if the underlying solver is set up to use it.

You seem to be mistaken, however, in this case, as the following example
shows.
> library(optimx)
> library(adagio)
> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Linux Mint 18.3

Matrix products: default
BLAS: /usr/lib/openblas-base/libblas.so.3
LAPACK: /usr/lib/libopenblasp-r0.2.18.so

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C
LC_TIME=en_CA.UTF-8
 [4] LC_COLLATE=en_CA.UTF-8     LC_MONETARY=en_CA.UTF-8
LC_MESSAGES=en_CA.UTF-8
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                  LC_ADDRESS=C

[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_CA.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] adagio_0.7.1     optimx_2018-7.10

loaded via a namespace (and not attached):
[1] compiler_3.5.1    tools_3.5.1      
numDeriv_2016.8-1
> test <- optimr(c(-1.2, 1), fn=fnRosenbrock, method="L-BFGS-B",
control=list(maxit=2, trace=2))
Parameter scaling:[1] 1 1
N = 2, M = 5 machine precision = 2.22045e-16
This problem is unconstrained.
final  value 4.120516
stopped after 3 iterations

So it seems to have worked, albeit one iteration later than specified
(the test is not always immediate).

Note optimr() is just a wrapper. There may be improvements that could be
made to use maxit and similar controls. Note that the code is in plain R
and optimr() was written in a way that allows it to be maintained. The
most difficult problem -- which may be related to the OP's complaint in
some cases -- is translating maxit to whatever is used by the solver,
then copying back the relevant information about number of iterations.

John Nash


On 2018-10-03 01:30 PM, Ahmed Attia wrote:
> The argument maxit used in libraries optimr, hydroPSO, and FME to
> control the maximal number of iterations to be performed does not work
> at ALL!!
> 
> This needs to be fixed...
> 
> 
> 
> mysamp <- function(n, m, s, lwr, upr, nnorm) {
>   samp <- rnorm(nnorm, m, s)
>   samp <- samp[samp >= lwr & samp <= upr]
>   if (length(samp) >= n) {
>     return(sample(samp, n))
>   }
>   stop(simpleError("Not enough values to sample from. Try increasing
nnorm."))
> }
> 
> x <- mysamp(1000,0.8117648,0.1281839,0,1,1000000)
> y <- mysamp(1000,0.7530346,0.1865946,0,1,1000000)
> 
> n <- length(x)
> xgroup <- c('a','b')
> ygroup <- c('c','d')
> 
> for (i in c(2:n)){
>   xgroup[i] <-  c('b')
>   ygroup[i] <-  c('d')
> }
> 
> dataset <- data.frame(x = x, y = y, xgroup = as.factor(xgroup),
>                       ygroup = as.factor(ygroup))
> dataset$ObsNo <- 1:n
> dataset <- dataset[order(dataset$x,decreasing=F),]
> 
> xopt <- c(dataset$x[1],0.95)    #Initial critical x,y values
> 
> 
> my.function <- function(xopt){
>   for (i in c(1:n)){
>     dataset$xgroup[i] <- if(dataset$x[i] < xopt[1]) 'a' else
'b'
>     dataset$ygroup[i] <- if(dataset$y[i] < xopt[2]) 'c' else
'd'
>     dataset$q.i[i] <- with(dataset, ifelse
>                            (dataset$xgroup[i]=='a' &
> dataset$ygroup[i]=='d', 1, 0))
>     dataset$q.ii[i] <- with(dataset, ifelse
>                             (dataset$xgroup[i]=='b' &
> dataset$ygroup[i]=='d', 1, 0))
>     dataset$q.iii[i] <- with(dataset, ifelse
>                              (dataset$xgroup[i]=='b' &
> dataset$ygroup[i]=='c', 1, 0))
>     dataset$q.iv[i] <- with(dataset, ifelse
>                             (dataset$xgroup[i]=='a' &
> dataset$ygroup[i]=='c', 1, 0))
> 
>     dataset$q.err[i] <- sum(dataset$q.i[i] + dataset$q.iii[i])
>   }
>   min.qerr <- sum(dataset$q.err)
>   q.I <- sum(dataset$q.i)
>   q.II <- sum(dataset$q.ii)
>   q.III <- sum(dataset$q.iii)
>   q.IV <- sum(dataset$q.iv)
>   q.Iandq.III <- sum(dataset$q.err)
>   print(c(q.I,
>           q.II,
>           q.III,
>           q.IV,
>           q.Iandq.III))
>   return(min.qerr)
> }
> 
> my.function(xopt)
> 
> #-------Algorithm
> 
> xmin=c(0,0.60)
> xmax=c(0.95,0.95)
> 
> res= hydroPSO(par=xopt,my.function,method="spso2011",
>               lower=xmin,upper=xmax,control=list(maxit=3))
> 
> 
> res= modFit(my.function,p=xopt,method="Pseudo",
>             lower=xmin,upper=xmax,control=list(numiter=3))
> 
> res= optimr(par=xopt,my.function,method="L-BFGS-B",
>             lower=xmin,upper=xmax,control=list(maxit=3,trace=T))
> 
> #Only the hjkb form dfoptim that has argument maxfeval and it works.
> 
> res=hjkb(par=xopt,my.function,
>          lower=xmin,upper=xmax,control=list(maxfeval=3,tol=2^-10,info=T))
> 
> 
> 
> 
> 
> 
> Ahmed Attia, Ph.D.
> Agronomist & Crop Modeler
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>