Skip to content

[R-sig-dyn-mod] Problem with simecol and ode - crash from R

8 messages · Thomas Petzoldt, Rainer M Krug, Soetaert, Karline

Original
Thomas Petzoldt
#
Rainer M Krug wrote:
Hi Rainer,

as you see yourself, this is not a problem of package simecol, it is a
problem of lsoda from the deSolve package. The reason is, that with the
parameters given, your simulation runs to minus infinity.

As a result, former versions of deSolve/lsoda produced a buffer overflow
and this resulted in a "stack smashing" error.

The "stack smashing" bug was solved in version 1.5 of deSolve, available
from CRAN since last weekend. Instead of breaking R lsoda now issues an
error message. I admit that this error message is not very helpful, but
there are several possibilities to diagnose the problem:

One option is to use a fixed step solver (rk4, euler) as you did, or a
different variable time step solver (rk alias ode45, vode, lsode or
daspk). Given your model and parameters, try:

solver(lgH) <- "vode"
plot(sim(lgH))

lsoda seems the only one that fails completely, but even that case can
be solved (in the diagnostic sense, i.e. by potentially loosing
accuracy) if you limit hmin, e.g.:

solver(lgH) <- "lsoda"
plot(sim(lgH, hmin = 1e-12))


Hope it helps and please update deSolve.

Thomas Petzoldt
Thomas Petzoldt
#
Rainer M Krug wrote:
Rainer,

sorry, but I cannot reproduce this crash. Not on R-2.9.0 or R-2.11-devel 
on Windows and not on R 2.7.1 / Debian Linux.

We had this problem in the past and I was sure that it was solved. Are 
you sure that you installed deSolve 1.5, i.e. was your CRAN mirror up to 
date?

Nevertheless, using another solver (e.g. vode) should be be a reasonable 
workaround for you and we keep it in mind for the next release.

Thank you for your report

Thomas P.
Soetaert, Karline
# 
Rainer,
 
I cannot recreate your crash; I get a lot of warnings and finally an error from lsoda, which terminates the simulation (but not R).  
The warnings are there because, between 43 and 44 simulation units the state variable becomes -Inf.  This is due to the 0-order consumption rate (h), so that "N" continues being consumed even if there is nothing to consume. This generates negative densities, which are also being consumed etc..., and N decreases exponentially unitl infinity.
 
 
A more appropriate model would be:
 
 
NoCrash <- function (time, N, parms) {
            with(as.list(parms),
                 {
                   dn <-  r * N * (1 - (N/K)) - h*(N>0)
                   if (! is.numeric(dn)) dn <- 0
                   list(dn)
                 }
                 )
          }

require(deSolve)
parms  = c(r = 0.1, K = 10000, h=200)
init   = c(N = 2000)

out <- ode(func=NoCrash, y = init, parms=parms, times=0:50, method="lsoda")
plot(out)
 
Due to the (N>0) there is only harvesting as long as N>0.
 
note that the plot(out) only works from deSolve version 1.5
 
 
 
Karline Soetaert

________________________________

From: r-sig-dynamic-models-bounces at r-project.org on behalf of Rainer M Krug
Sent: Wed 10/21/2009 12:14 PM
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] Problem with simecol and ode - crash from R



Hi

I have a problem with simecol, particularly with the logistic growth model
with harvesting: with the parameters as given below, it crashes R.

When I change h to values below 100, it works, changing the solver to
"euler" or "rk4", it works

Any idea what is happening, or am I doing something wrong?


Thanks,

Rainer

library(simecol)
lgH <- new(
          "odeModel",
          main = function (time, init, parms) {
            with(
                 as.list(c(init, parms)),
                 {
                   dn <-  r * N * (1 - (N/K)) - h
                   list(dn)
                 }
                 )
          },
          parms  = c(r = 0.1, K = 10000, h=200),
          times  = c(from = 0, to = 100, by = 0.5),
          init   = c(N = 2000),
          solver = "lsoda"
)
lgH <-  sim(lgH)
*** stack smashing detected ***: /usr/lib/R/bin/exec/R terminated
======= Backtrace: =========
/lib/tls/i686/cmov/libc.so.6(__fortify_fail+0x48)[0xb7c69138]
/lib/tls/i686/cmov/libc.so.6(__fortify_fail+0x0)[0xb7c690f0]
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/deSolve/libs/deSolve.so[0xb74822b4]
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/deSolve/libs/deSolve.so(dlsoda_+0x1279)[0xb746c329]
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/deSolve/libs/deSolve.so(call_lsoda+0xa2d)[0xb743628d]
/usr/lib/R/lib/libR.so[0xb7d6999c]
/usr/lib/R/lib/libR.so(Rf_eval+0x714)[0xb7d92fb4]
/usr/lib/R/lib/libR.so[0xb7d962fe]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so[0xb7d94530]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so(Rf_applyClosure+0x2ac)[0xb7d96c6c]
/usr/lib/R/lib/libR.so(Rf_eval+0x349)[0xb7d92be9]
/usr/lib/R/lib/libR.so[0xb7d23ce6]
/usr/lib/R/lib/libR.so[0xb7dd2d7a]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so[0xb7d94530]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so(Rf_applyClosure+0x2ac)[0xb7d96c6c]
/usr/lib/R/lib/libR.so(Rf_eval+0x349)[0xb7d92be9]
/usr/lib/R/lib/libR.so[0xb7d962fe]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so[0xb7d94530]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so[0xb7d96487]
/usr/lib/R/lib/libR.so(R_execMethod+0x2c1)[0xb7d96901]
/usr/lib/R/library/methods/libs/methods.so[0xb764b66a]
/usr/lib/R/lib/libR.so[0xb7dd47ab]
/usr/lib/R/lib/libR.so(Rf_eval+0x598)[0xb7d92e38]
/usr/lib/R/lib/libR.so(Rf_applyClosure+0x2ac)[0xb7d96c6c]
/usr/lib/R/lib/libR.so(Rf_eval+0x349)[0xb7d92be9]
/usr/lib/R/lib/libR.so[0xb7d962fe]
/usr/lib/R/lib/libR.so(Rf_eval+0x451)[0xb7d92cf1]
/usr/lib/R/lib/libR.so(Rf_ReplIteration+0x1c5)[0xb7dc0685]
/usr/lib/R/lib/libR.so(run_Rmainloop+0x102)[0xb7dc0a62]
/usr/lib/R/lib/libR.so(Rf_mainloop+0x1c)[0xb7dc0abc]
/usr/lib/R/bin/exec/R(main+0x46)[0x8048776]
/lib/tls/i686/cmov/libc.so.6(__libc_start_main+0xe0)[0xb7b92450]
/usr/lib/R/bin/exec/R[0x8048691]
======= Memory map: ========
08048000-08049000 r-xp 00000000 08:02 811281     /usr/lib/R/bin/exec/R
08049000-0804a000 rw-p 00000000 08:02 811281     /usr/lib/R/bin/exec/R
0804a000-08c8b000 rw-p 0804a000 00:00 0          [heap]
b742f000-b7485000 r-xp 00000000 08:03 10371078
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/deSolve/libs/deSolve.so
b7485000-b7486000 rw-p 00055000 08:03 10371078
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/deSolve/libs/deSolve.so
b7486000-b7487000 rw-p b7486000 00:00 0
b7487000-b74e9000 r-xp 00000000 08:02 1444606
/usr/lib/R/library/stats/libs/stats.so
b74e9000-b74eb000 rw-p 00062000 08:02 1444606
/usr/lib/R/library/stats/libs/stats.so
b74eb000-b750c000 r-xp 00000000 08:02 1810475
/usr/lib/R/library/grDevices/libs/grDevices.so
b750c000-b750d000 rw-p 00020000 08:02 1810475
/usr/lib/R/library/grDevices/libs/grDevices.so
b7555000-b755f000 r-xp 00000000 08:02 1146976    /lib/libgcc_s.so.1
b755f000-b7560000 rw-p 0000a000 08:02 1146976    /lib/libgcc_s.so.1
b7560000-b7562000 r-xp 00000000 08:03 4572127
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/simecol/libs/simecol.so
b7562000-b7563000 rw-p 00001000 08:03 4572127
/home/rkrug/R/i486-pc-linux-gnu-library/2.9/simecol/libs/simecol.so
b7563000-b75d5000 rw-p b7563000 00:00 0
b75d5000-b75d6000 r-xp 00000000 08:02 1226930    /usr/lib/gconv/ISO8859-1.so
b75d6000-b75d8000 rw-p 00000000 08:02 1226930    /usr/lib/gconv/ISO8859-1.so
b75d8000-b7648000 rw-p b75d8000 00:00 0
b7648000-b764e000 r-xp 00000000 08:02 1794882
/usr/lib/R/library/methods/libs/methods.so
b764e000-b764f000 rw-p 00005000 08:02 1794882
/usr/lib/R/library/methods/libs/methods.so
b764f000-b767b000 rw-p b764f000 00:00 0
b767b000-b7684000 r-xp 00000000 08:02 467611     /lib/tls/i686/cmov/
libnss_files-2.7.so
b7684000-b7686000 rw-p 00008000 08:02 467611     /lib/tls/i686/cmov/
libnss_files-2.7.so
b7686000-b768e000 r-xp 00000000 08:02 467613     /lib/tls/i686/cmov/
libnss_nis-2.7.so
b768e000-b7690000 rw-p 00007000 08:02 467613     /lib/tls/i686/cmov/
libnss_nis-2.7.so
b7690000-b76a4000 r-xp 00000000 08:02 467608     /lib/tls/i686/cmov/li
Process R aborted at Wed Oct 21 12:05:59 2009

--
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology,
UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Natural Sciences Building
Office Suite 2039
Stellenbosch University
Main Campus, Merriman Avenue
Stellenbosch
South Africa

Cell:           +27 - (0)83 9479 042
Fax:            +27 - (0)86 516 2782
Fax:            +49 - (0)721 151 334 888
email:          Rainer at krugs.de

Skype:          RMkrug
Google:         R.M.Krug at gmail.com


_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models


-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/ms-tnef
Size: 11212 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-dynamic-models/attachments/20091021/9face222/attachment.bin>
Thomas Petzoldt
#
Soetaert, Karline wrote:
[...]

 > A more appropriate model would be:
 >
 >
 > NoCrash <- function (time, N, parms) {
 >             with(as.list(parms),
 >                  {
 >                    dn <-  r * N * (1 - (N/K)) - h*(N>0)
 >                    if (! is.numeric(dn)) dn <- 0
 >                    list(dn)
 >                  }
 >                  )
 >           }
 >
 > require(deSolve)
 > parms  = c(r = 0.1, K = 10000, h=200)
 > init   = c(N = 2000)

Well, this works in practice, but allow me to two additional points:

First, the line "if (! is.numeric(dn)) dn <- 0" is redundant, probably a 
residual of testing?

Second, instead of a hard limit (that is of course a good test case for 
the solvers), I usually prefer a more smooth limit, e.g:

NoCrash <- function (time, N, parms) {
              with(as.list(parms), {
                    dn <-  r * N * (1 - (N/K)) - h * N / (N + hlim)
                    list(dn)
              })
            }

parms  <- c(r = 0.1, K = 10000, h = 200, hlim = 1)
init   <- c(N = 2000)

out <- ode(func=NoCrash, y = init, parms = parms,
   times = seq(0,50,.01), method = "lsoda")

plot(out, type = "l")

You can of course try different values of hlim, e.g. hlim=100.

Thomas P.