diff --git a/.Rbuildignore b/.Rbuildignore index 5d79eddc..5954eef4 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -16,3 +16,6 @@ ^\.covrignore$ ^LICENSE$ ^\.zenodo\.json$ + +^\.vscode$ +^source_odeint\.R diff --git a/DESCRIPTION b/DESCRIPTION index fdb6c78d..813d99dc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: DAISIE Type: Package Title: Dynamical Assembly of Islands by Speciation, Immigration and Extinction -Version: 4.3.2 -Date: 2023-02-24 +Version: 4.3.3 +Date: 2023-03-10 Depends: R (>= 4.1.0) biocViews: Imports: @@ -111,7 +111,6 @@ Description: Simulates and computes the (maximum) likelihood of a dynamical extinction. See e.g. Valente et al. 2015. Ecology Letters 18: 844-852, . NeedsCompilation: yes -SystemRequirements: C++14 Encoding: UTF-8 VignetteBuilder: knitr URL: https://github.com/rsetienne/DAISIE diff --git a/NEWS.md b/NEWS.md index fd38c444..eaf98bc1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,22 @@ +# DAISIE 4.3.3 + +* Address problem detected by valgrind: unitialized member variable +bulirsch_stoer<>::m_dt_last. + * Patched version of `boost/numeric/odeint/stepper/bulirsch_stoer.hpp`. This + is done by including the patched header file `src/patched_bulrisch_stoer.h` + before `boost/numeric/odeint` to shadow + `boost/numeric/odeint/stepper/bulrisch_stoer.hpp` + The issue is *not* fixed in BOOST_VERSION == 1.81.0. + Must check for fixes in upcomming boost (BH) releases. +* Fix tab spacing in `src/DAISIE_loglik_rhs_FORTRAN.f95`. +* Re-implement clang16 `-D_HAS_AUTO_PTR_ETC` flag fix via Makevars[.win] to +comply with CRAN requests and to force C++ standard C++14 without using +SystemRequirements line in DESCRIPTION, at CRAN's request. +* Skip plotting tests that cause problems in headless systems. These should +be run manually. +* Skip tests of as-of-now experimental steppers from ODEINT rosenbrock4 and +adams bashfort moulton 1 because they are too slow. + # DAISIE 4.3.2 * Apply CRAN suggested fixes to clang16 issues with deprecated C++ functions included the Boost library, which are used in some of the stepper functions. diff --git a/src/DAISIE_loglik_rhs_FORTRAN.f95 b/src/DAISIE_loglik_rhs_FORTRAN.f95 index 2bb49994..5c121830 100644 --- a/src/DAISIE_loglik_rhs_FORTRAN.f95 +++ b/src/DAISIE_loglik_rhs_FORTRAN.f95 @@ -82,8 +82,8 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) INTEGER :: neq, ip(*), i, ii DOUBLE PRECISION :: t, Conc(2 * N + 1), dConc(2 * N + 1), yout(*) DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3 - INTEGER :: il1(N), il2(N), il3in3(N), il4(N) - INTEGER :: in1(N), in2ix2(N) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N) INTEGER :: ix1(N), ix3(N), ix4(N) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) @@ -131,30 +131,30 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) xx1(1) = 0 xx1(2) = 0 xx2(1) = 0 - xx2(2) = 0 + xx2(2) = 0 DO I = 3, N + 2 - xx1(I) = Conc(I - 2) - xx2(I) = Conc(N + I - 2) - il1(I - 2) = I + kk - 1 - il2(I - 2) = I + kk + 1 - il3in3(I - 2) = I + kk - il4(I - 2) = I + kk - 2 - in1(I - 2) = I + 2 * kk - 1 - in2ix2(I - 2) = I + 1 - ix1(I - 2) = I - 1 - ix3(I - 2) = I - ix4(I - 2) = I - 2 - ENDDO - xx1(N + 3) = 0 - xx2(N + 3) = 0 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 xx3 = Conc(2 * N + 1) DO I = 1, N + 4 + 2 * kk - laavec(I) = P(I) - lacvec(I) = P(I + N + 4 + 2 * kk) - muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO ! dx1 = laavec[il1 + 1] * xx2[ix1] + @@ -172,15 +172,15 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) ! dx3 <- 0 DO I = 1, N - dConc(I) = & - laavec(il1(I) + 1) * xx2(ix1(I)) + & - lacvec(il4(I) + 1) * xx2(ix4(I)) + & - muvec(il2(I) + 1) * xx2(ix3(I)) + & - lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & - (muvec(il3in3(I)) + lacvec(il3in3(I))) * & - nn(il3in3(I)) * xx1(ix3(I)) - & - gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(I) = & + laavec(il1(I) + 1) * xx2(ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I) + 1) * xx2(ix3(I)) + & + lacvec(il1(I)) * nn(in1(I)) * xx1(ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) dConc(N + I) = & gamvec(il3in3(I)) * xx1(ix3(I)) + & lacvec(il1(I) + 1) * nn(in1(I)) * xx2(ix1(I)) + & @@ -188,7 +188,7 @@ SUBROUTINE daisie_runmod (neq, t, Conc, dConc, yout, ip) (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & nn(il3in3(I) + 1) * xx2(ix3(I)) - & laavec(il3in3(I) + 1) * xx2(ix3(I)) - dConc(2 * N + 1) = 0 + dConc(2 * N + 1) = 0 ENDDO @@ -209,8 +209,8 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) INTEGER :: neq, ip(*), i, ii DOUBLE PRECISION :: t, Conc(4 * N), dConc(4 * N), yout(*) DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3(N + 3), xx4(N + 3) - INTEGER :: il1(N), il2(N), il3in3(N), il4(N) - INTEGER :: in1(N), in2ix2(N), in4ix1(N) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N), in4ix1(N) INTEGER :: ix3(N), ix4(N) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) @@ -262,37 +262,37 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) xx1(1) = 0 xx1(2) = 0 xx2(1) = 0 - xx2(2) = 0 + xx2(2) = 0 xx3(1) = 0 - xx3(2) = 0 + xx3(2) = 0 xx4(1) = 0 - xx4(2) = 0 + xx4(2) = 0 DO I = 3, N + 2 - xx1(I) = Conc(I - 2) - xx2(I) = Conc(N + I - 2) - xx3(I) = Conc(2 * N + I - 2) - xx4(I) = Conc(3 * N + I - 2) - il1(I - 2) = I + kk - 1 - il2(I - 2) = I + kk + 1 - il3in3(I - 2) = I + kk - il4(I - 2) = I + kk - 2 - in1(I - 2) = I + 2 * kk - 1 - in2ix2(I - 2) = I + 1 - in4ix1(I - 2) = I - 1 - ix3(I - 2) = I - ix4(I - 2) = I - 2 - ENDDO - xx1(N + 3) = 0 - xx2(N + 3) = 0 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + xx3(I) = Conc(2 * N + I - 2) + xx4(I) = Conc(3 * N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + in4ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 xx3(N + 3) = 0 xx4(N + 3) = 0 DO I = 1, N + 4 + 2 * kk - laavec(I) = P(I) - lacvec(I) = P(I + N + 4 + 2 * kk) - muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO ! dx1 <- lacvec[il1] * xx1[ix1] + @@ -304,15 +304,15 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) ! -gamvec[il3] * xx1[ix3] DO I = 1, N - dConc(I) = & - lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & - laavec(il1(I) + 1) * xx2(in4ix1(I)) + & - lacvec(il4(I) + 1) * xx2(ix4(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) + & - muvec(il3in3(I) + 1) * xx2(ix3(I)) - & - (muvec(il3in3(I)) + lacvec(il3in3(I))) * & - nn(il3in3(I)) * xx1(ix3(I)) - & - gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(I) = & + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & + laavec(il1(I) + 1) * xx2(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) + & + muvec(il3in3(I) + 1) * xx2(ix3(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) ! dx2 <- gamvec[il3] * xx1[ix3] + ! gamvec[il3] * xx3[ix3] + @@ -326,11 +326,11 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) gamvec(il3in3(I)) * xx1(ix3(I)) + & gamvec(il3in3(I)) * xx3(ix3(I)) + & gamvec(il3in3(I) + 1) * xx4(ix3(I)) + & - lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & - muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & - (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & - nn(il3in3(I) + 1) * xx2(ix3(I)) - & - laavec(il3in3(I) + 1) * xx2(ix3(I)) + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) ! dx3 <- lacvec[il1] * nn[in1] * xx3[ix1] + ! laavec[il1 + 1] * xx4[ix1] + @@ -355,7 +355,7 @@ SUBROUTINE daisie_runmod1 (neq, t, Conc, dConc, yout, ip) ! -(lacvec[il3 + 1] + muvec[il3 + 1]) * nn[in3 + 1] * xx4[ix3] + ! -gamvec[il3 + 1] * xx4[ix3] - dConc(3 * N + I) = & + dConc(3 * N + I) = & lacvec(il1(I) + 1) * nn(in1(I)) * xx4(in4ix1(I)) + & muvec(il2(I) + 1) * nn(in2ix2(I)) * xx4(in2ix2(I)) - & (lacvec(il3in3(I) + 1) + muvec(il3in3(I) + 1)) * & @@ -382,8 +382,8 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) INTEGER :: neq, ip(*), i, ii DOUBLE PRECISION :: t, Conc(3 * N), dConc(3 * N), yout(*) DOUBLE PRECISION :: xx1(N + 3), xx2(N + 3), xx3(N + 3) - INTEGER :: il1(N), il2(N), il3in3(N), il4(N) - INTEGER :: in1(N), in2ix2(N), in4ix1(N) + INTEGER :: il1(N), il2(N), il3in3(N), il4(N) + INTEGER :: in1(N), in2ix2(N), in4ix1(N) INTEGER :: ix3(N), ix4(N) DOUBLE PRECISION :: laavec(N + 4 + 2 * kk),lacvec(N + 4 + 2 * kk) DOUBLE PRECISION :: muvec(N + 4 + 2 * kk),gamvec(N + 4 + 2 * kk) @@ -433,33 +433,33 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) xx1(1) = 0 xx1(2) = 0 xx2(1) = 0 - xx2(2) = 0 + xx2(2) = 0 xx3(1) = 0 - xx3(2) = 0 + xx3(2) = 0 DO I = 3, N + 2 - xx1(I) = Conc(I - 2) - xx2(I) = Conc(N + I - 2) - xx3(I) = Conc(2 * N + I - 2) - il1(I - 2) = I + kk - 1 - il2(I - 2) = I + kk + 1 - il3in3(I - 2) = I + kk - il4(I - 2) = I + kk - 2 - in1(I - 2) = I + 2 * kk - 1 - in2ix2(I - 2) = I + 1 - in4ix1(I - 2) = I - 1 - ix3(I - 2) = I - ix4(I - 2) = I - 2 - ENDDO - xx1(N + 3) = 0 - xx2(N + 3) = 0 + xx1(I) = Conc(I - 2) + xx2(I) = Conc(N + I - 2) + xx3(I) = Conc(2 * N + I - 2) + il1(I - 2) = I + kk - 1 + il2(I - 2) = I + kk + 1 + il3in3(I - 2) = I + kk + il4(I - 2) = I + kk - 2 + in1(I - 2) = I + 2 * kk - 1 + in2ix2(I - 2) = I + 1 + in4ix1(I - 2) = I - 1 + ix3(I - 2) = I + ix4(I - 2) = I - 2 + ENDDO + xx1(N + 3) = 0 + xx2(N + 3) = 0 xx3(N + 3) = 0 DO I = 1, N + 4 + 2 * kk - laavec(I) = P(I) - lacvec(I) = P(I + N + 4 + 2 * kk) - muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) - gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) - nn(I) = P(I + 4 * (N + 4 + 2 * kk)) + laavec(I) = P(I) + lacvec(I) = P(I + N + 4 + 2 * kk) + muvec(I) = P(I + 2 * (N + 4 + 2 * kk)) + gamvec(I) = P(I + 3 * (N + 4 + 2 * kk)) + nn(I) = P(I + 4 * (N + 4 + 2 * kk)) ENDDO ! dx1 = (laavec[il3] * xx3[ix3] + @@ -484,33 +484,33 @@ SUBROUTINE daisie_runmod2 (neq, t, Conc, dConc, yout, ip) ! -(laavec[il3] + gamvec[il3]) * xx3[ix3] DO I = 1, N - dConc(I) = 0 - IF(kk .EQ. 1) THEN - dConc(I) = laavec(il3in3(I)) * xx3(ix3(I)) + & - 2 * lacvec(il1(I)) * xx3(in4ix1(I)) - ENDIF - dConc(I) = dConc(I) + & - laavec(il1(I) + 1) * xx2(in4ix1(I)) + & - lacvec(il4(I) + 1) * xx2(ix4(I)) + & - muvec(il2(I) + 1) * xx2(ix3(I)) + & - lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & - (muvec(il3in3(I)) + lacvec(il3in3(I))) * & - nn(il3in3(I)) * xx1(ix3(I)) - & - gamvec(il3in3(I)) * xx1(ix3(I)) - dConc(N + I) = & - gamvec(il3in3(I)) * xx1(ix3(I)) + & - lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & - muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & - (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & - nn(il3in3(I) + 1) * xx2(ix3(I)) - & - laavec(il3in3(I) + 1) * xx2(ix3(I)) - dConc(2 * N + I) = & - lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) + & - muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - & - (lacvec(il3in3(I)) + muvec(il3in3(I))) * & - nn(il3in3(I)) * xx3(ix3(I)) - & - (laavec(il3in3(I)) + gamvec(il3in3(I))) * xx3(ix3(I)) + dConc(I) = 0 + IF(kk .EQ. 1) THEN + dConc(I) = laavec(il3in3(I)) * xx3(ix3(I)) + & + 2 * lacvec(il1(I)) * xx3(in4ix1(I)) + ENDIF + dConc(I) = dConc(I) + & + laavec(il1(I) + 1) * xx2(in4ix1(I)) + & + lacvec(il4(I) + 1) * xx2(ix4(I)) + & + muvec(il2(I) + 1) * xx2(ix3(I)) + & + lacvec(il1(I)) * nn(in1(I)) * xx1(in4ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx1(in2ix2(I)) - & + (muvec(il3in3(I)) + lacvec(il3in3(I))) * & + nn(il3in3(I)) * xx1(ix3(I)) - & + gamvec(il3in3(I)) * xx1(ix3(I)) + dConc(N + I) = & + gamvec(il3in3(I)) * xx1(ix3(I)) + & + lacvec(il1(I) + 1) * nn(in1(I)) * xx2(in4ix1(I)) + & + muvec(il2(I) + 1) * nn(in2ix2(I)) * xx2(in2ix2(I)) - & + (muvec(il3in3(I) + 1) + lacvec(il3in3(I) + 1)) * & + nn(il3in3(I) + 1) * xx2(ix3(I)) - & + laavec(il3in3(I) + 1) * xx2(ix3(I)) + dConc(2 * N + I) = & + lacvec(il1(I)) * nn(in4ix1(I)) * xx3(in4ix1(I)) + & + muvec(il2(I)) * nn(in2ix2(I)) * xx3(in2ix2(I)) - & + (lacvec(il3in3(I)) + muvec(il3in3(I))) * & + nn(il3in3(I)) * xx3(ix3(I)) - & + (laavec(il3in3(I)) + gamvec(il3in3(I))) * xx3(ix3(I)) ENDDO END SUBROUTINE daisie_runmod2 diff --git a/src/DAISIE_odeint.h b/src/DAISIE_odeint.h index bb84dd62..a2d531a5 100644 --- a/src/DAISIE_odeint.h +++ b/src/DAISIE_odeint.h @@ -6,6 +6,7 @@ #include "config.h" #include "ublas_types.h" +#include "patched_bulrisch_stoer.h" // shadow buggy boost header #include #include #include diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 00000000..1fada8d1 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +CXX_STD = CXX14 +PKG_CPPFLAGS = -D_HAS_AUTO_PTR_ETC=0 diff --git a/src/Makevars.win b/src/Makevars.win new file mode 100644 index 00000000..1fada8d1 --- /dev/null +++ b/src/Makevars.win @@ -0,0 +1,2 @@ +CXX_STD = CXX14 +PKG_CPPFLAGS = -D_HAS_AUTO_PTR_ETC=0 diff --git a/src/config.h b/src/config.h index 348a5311..32cfa85b 100644 --- a/src/config.h +++ b/src/config.h @@ -5,18 +5,24 @@ // Special case to make use of some steppers that would include // boost/functional.hpp +// moved to Makevars[.win] +/* #if __cplusplus >= 201703L -#ifdef _HAS_AUTO_PTR_ETC -#undef _HAS_AUTO_PTR_ETC -#endif -#define _HAS_AUTO_PTR_ETC 0 +# ifdef _HAS_AUTO_PTR_ETC +# undef _HAS_AUTO_PTR_ETC +# endif +# define _HAS_AUTO_PTR_ETC 0 #endif + */ // Special case to make use of some steppers that would include // boost/get_pointer.hpp -#ifdef BOOST_NO_AUTO_PTR -#undef BOOST_NO_AUTO_PTR +#ifndef BOOST_NO_AUTO_PTR +# define BOOST_NO_AUTO_PTR #endif -#define BOOST_NO_AUTO_PTR + +// uncomment if unitialized member variable bulirsch_stoer::m_dt_last +// is fixed in boost (BH) +#define USE_BULRISCH_STOER_PATCH #endif diff --git a/src/patched_bulrisch_stoer.h b/src/patched_bulrisch_stoer.h new file mode 100644 index 00000000..cad4f521 --- /dev/null +++ b/src/patched_bulrisch_stoer.h @@ -0,0 +1,661 @@ +/* + Patched version of + boost/numeric/odeint/stepper/bulirsch_stoer.hpp + + Addresses unitialized member variable bulirsch_stoer<>::m_dt_last. + + Include this header before boost/numeric/odeint to shadow + boost/numeriuc/odeint/stepper/bulrisch_stoer.hpp + + The issue is *not* fixed in BOOST_VERSION == 1.81.0. + We need to check for fixes in upcomming boost (BH) releases. + + Hanno Hildenbrandt 2023 +*/ + +/* + [auto_generated] + boost/numeric/odeint/stepper/bulirsch_stoer.hpp + + [begin_description] + Implementation of the Burlish-Stoer method. As described in + Ernst Hairer, Syvert Paul Norsett, Gerhard Wanner + Solving Ordinary Differential Equations I. Nonstiff Problems. + Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Second revised edition 1993. + [end_description] + + Copyright 2011-2013 Mario Mulansky + Copyright 2011-2013 Karsten Ahnert + Copyright 2012 Christoph Koke + + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or + copy at http://www.boost.org/LICENSE_1_0.txt) +*/ + +#ifdef USE_BULRISCH_STOER_PATCH + +#ifndef BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED +#define BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED + + +#include + +#include + +#include // for min/max guidelines + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace boost { +namespace numeric { +namespace odeint { + +template< + class State , + class Value = double , + class Deriv = State , + class Time = Value , + class Algebra = typename algebra_dispatcher< State >::algebra_type , + class Operations = typename operations_dispatcher< State >::operations_type , + class Resizer = initially_resizer + > +class bulirsch_stoer { + +public: + + typedef State state_type; + typedef Value value_type; + typedef Deriv deriv_type; + typedef Time time_type; + typedef Algebra algebra_type; + typedef Operations operations_type; + typedef Resizer resizer_type; +#ifndef DOXYGEN_SKIP + typedef state_wrapper< state_type > wrapped_state_type; + typedef state_wrapper< deriv_type > wrapped_deriv_type; + typedef controlled_stepper_tag stepper_category; + + typedef bulirsch_stoer< State , Value , Deriv , Time , Algebra , Operations , Resizer > controlled_error_bs_type; + + typedef typename inverse_time< time_type >::type inv_time_type; + + typedef std::vector< value_type > value_vector; + typedef std::vector< time_type > time_vector; + typedef std::vector< inv_time_type > inv_time_vector; //should be 1/time_type for boost.units + typedef std::vector< value_vector > value_matrix; + typedef std::vector< size_t > int_vector; + typedef std::vector< wrapped_state_type > state_table_type; +#endif //DOXYGEN_SKIP + const static size_t m_k_max = 8; + + bulirsch_stoer( + value_type eps_abs = 1E-6 , value_type eps_rel = 1E-6 , + value_type factor_x = 1.0 , value_type factor_dxdt = 1.0 , + time_type max_dt = static_cast(0)) + : m_error_checker( eps_abs , eps_rel , factor_x, factor_dxdt ) , m_midpoint() , + m_last_step_rejected( false ) , m_first( true ) , + m_dt_last(static_cast(-1)) , // !!! patched !!! + m_max_dt(max_dt) , + m_interval_sequence( m_k_max+1 ) , + m_coeff( m_k_max+1 ) , + m_cost( m_k_max+1 ) , + m_facmin_table( m_k_max+1 ) , + m_table( m_k_max ) , + STEPFAC1( 0.65 ) , STEPFAC2( 0.94 ) , STEPFAC3( 0.02 ) , STEPFAC4( 4.0 ) , KFAC1( 0.8 ) , KFAC2( 0.9 ) + { + BOOST_USING_STD_MIN(); + BOOST_USING_STD_MAX(); + /* initialize sequence of stage numbers and work */ + for( unsigned short i = 0; i < m_k_max+1; i++ ) + { + m_interval_sequence[i] = 2 * (i+1); + if( i == 0 ) + m_cost[i] = m_interval_sequence[i]; + else + m_cost[i] = m_cost[i-1] + m_interval_sequence[i]; + m_coeff[i].resize(i); + m_facmin_table[i] = pow BOOST_PREVENT_MACRO_SUBSTITUTION( STEPFAC3 , static_cast< value_type >(1) / static_cast< value_type >( 2*i+1 ) ); + for( size_t k = 0 ; k < i ; ++k ) + { + const value_type r = static_cast< value_type >( m_interval_sequence[i] ) / static_cast< value_type >( m_interval_sequence[k] ); + m_coeff[i][k] = 1.0 / ( r*r - static_cast< value_type >( 1.0 ) ); // coefficients for extrapolation + } + } + reset(); + } + + + /* + * Version 1 : try_step( sys , x , t , dt ) + * + * The overloads are needed to solve the forwarding problem + */ + template< class System , class StateInOut > + controlled_step_result try_step( System system , StateInOut &x , time_type &t , time_type &dt ) + { + return try_step_v1( system , x , t, dt ); + } + + /** + * \brief Second version to solve the forwarding problem, can be used with Boost.Range as StateInOut. + */ + template< class System , class StateInOut > + controlled_step_result try_step( System system , const StateInOut &x , time_type &t , time_type &dt ) + { + return try_step_v1( system , x , t, dt ); + } + + /* + * Version 2 : try_step( sys , x , dxdt , t , dt ) + * + * this version does not solve the forwarding problem, boost.range can not be used + */ + template< class System , class StateInOut , class DerivIn > + controlled_step_result try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) + { + m_xnew_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_xnew< StateInOut > , detail::ref( *this ) , detail::_1 ) ); + controlled_step_result res = try_step( system , x , dxdt , t , m_xnew.m_v , dt ); + if( res == success ) + { + boost::numeric::odeint::copy( m_xnew.m_v , x ); + } + return res; + } + + /* + * Version 3 : try_step( sys , in , t , out , dt ) + * + * this version does not solve the forwarding problem, boost.range can not be used + */ + template< class System , class StateIn , class StateOut > + typename boost::disable_if< boost::is_same< StateIn , time_type > , controlled_step_result >::type + try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) + { + typename odeint::unwrap_reference< System >::type &sys = system; + m_dxdt_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateIn > , detail::ref( *this ) , detail::_1 ) ); + sys( in , m_dxdt.m_v , t ); + return try_step( system , in , m_dxdt.m_v , t , out , dt ); + } + + + /* + * Full version : try_step( sys , in , dxdt_in , t , out , dt ) + * + * contains the actual implementation + */ + template< class System , class StateIn , class DerivIn , class StateOut > + controlled_step_result try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) + { + if( m_max_dt != static_cast(0) && detail::less_with_sign(m_max_dt, dt, dt) ) + { + // given step size is bigger then max_dt + // set limit and return fail + dt = m_max_dt; + return fail; + } + + BOOST_USING_STD_MIN(); + BOOST_USING_STD_MAX(); + + static const value_type val1( 1.0 ); + + if( m_resizer.adjust_size( in , detail::bind( &controlled_error_bs_type::template resize_impl< StateIn > , detail::ref( *this ) , detail::_1 ) ) ) + { + reset(); // system resized -> reset + } + + if( dt != m_dt_last ) + { + reset(); // step size changed from outside -> reset + } + + bool reject( true ); + + time_vector h_opt( m_k_max+1 ); + inv_time_vector work( m_k_max+1 ); + + time_type new_h = dt; + + /* m_current_k_opt is the estimated current optimal stage number */ + for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ ) + { + /* the stage counts are stored in m_interval_sequence */ + m_midpoint.set_steps( m_interval_sequence[k] ); + if( k == 0 ) + { + m_midpoint.do_step( system , in , dxdt , t , out , dt ); + /* the first step, nothing more to do */ + } + else + { + m_midpoint.do_step( system , in , dxdt , t , m_table[k-1].m_v , dt ); + extrapolate( k , m_table , m_coeff , out ); + // get error estimate + m_algebra.for_each3( m_err.m_v , out , m_table[0].m_v , + typename operations_type::template scale_sum2< value_type , value_type >( val1 , -val1 ) ); + const value_type error = m_error_checker.error( m_algebra , in , dxdt , m_err.m_v , dt ); + h_opt[k] = calc_h_opt( dt , error , k ); + work[k] = static_cast( m_cost[k] ) / h_opt[k]; + + if( (k == m_current_k_opt-1) || m_first ) + { // convergence before k_opt ? + if( error < 1.0 ) + { + //convergence + reject = false; + if( (work[k] < KFAC2*work[k-1]) || (m_current_k_opt <= 2) ) + { + // leave order as is (except we were in first round) + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(k)+1 ) ); + new_h = h_opt[k]; + new_h *= static_cast( m_cost[k+1] ) / static_cast( m_cost[k] ); + } else { + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max)-1 , max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(k) ) ); + new_h = h_opt[k]; + } + break; + } + else if( should_reject( error , k ) && !m_first ) + { + reject = true; + new_h = h_opt[k]; + break; + } + } + if( k == m_current_k_opt ) + { // convergence at k_opt ? + if( error < 1.0 ) + { + //convergence + reject = false; + if( (work[k-1] < KFAC2*work[k]) ) + { + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(m_current_k_opt)-1 ); + new_h = h_opt[m_current_k_opt]; + } + else if( (work[k] < KFAC2*work[k-1]) && !m_last_step_rejected ) + { + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max-1) , static_cast(m_current_k_opt)+1 ); + new_h = h_opt[k]; + new_h *= static_cast(m_cost[m_current_k_opt])/static_cast(m_cost[k]); + } else + new_h = h_opt[m_current_k_opt]; + break; + } + else if( should_reject( error , k ) ) + { + reject = true; + new_h = h_opt[m_current_k_opt]; + break; + } + } + if( k == m_current_k_opt+1 ) + { // convergence at k_opt+1 ? + if( error < 1.0 ) + { //convergence + reject = false; + if( work[k-2] < KFAC2*work[k-1] ) + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( 2 , static_cast(m_current_k_opt)-1 ); + if( (work[k] < KFAC2*work[m_current_k_opt]) && !m_last_step_rejected ) + m_current_k_opt = min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(m_k_max)-1 , static_cast(k) ); + new_h = h_opt[m_current_k_opt]; + } else + { + reject = true; + new_h = h_opt[m_current_k_opt]; + } + break; + } + } + } + + if( !reject ) + { + t += dt; + } + + if( !m_last_step_rejected || boost::numeric::odeint::detail::less_with_sign(new_h, dt, dt) ) + { + // limit step size + if( m_max_dt != static_cast(0) ) + { + new_h = detail::min_abs(m_max_dt, new_h); + } + m_dt_last = new_h; + dt = new_h; + } + + m_last_step_rejected = reject; + m_first = false; + + if( reject ) + return fail; + else + return success; + } + + /** \brief Resets the internal state of the stepper */ + void reset() + { + m_first = true; + m_last_step_rejected = false; + // crude estimate of optimal order + m_current_k_opt = 4; + /* no calculation because log10 might not exist for value_type! + const value_type logfact( -log10( max BOOST_PREVENT_MACRO_SUBSTITUTION( eps_rel , static_cast< value_type >(1.0E-12) ) ) * 0.6 + 0.5 ); + m_current_k_opt = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( 1 ) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast( m_k_max-1 ) , logfact )); + */ + } + + + /* Resizer methods */ + + template< class StateIn > + void adjust_size( const StateIn &x ) + { + resize_m_dxdt( x ); + resize_m_xnew( x ); + resize_impl( x ); + m_midpoint.adjust_size( x ); + } + + +private: + + template< class StateIn > + bool resize_m_dxdt( const StateIn &x ) + { + return adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable::type() ); + } + + template< class StateIn > + bool resize_m_xnew( const StateIn &x ) + { + return adjust_size_by_resizeability( m_xnew , x , typename is_resizeable::type() ); + } + + template< class StateIn > + bool resize_impl( const StateIn &x ) + { + bool resized( false ); + for( size_t i = 0 ; i < m_k_max ; ++i ) + resized |= adjust_size_by_resizeability( m_table[i] , x , typename is_resizeable::type() ); + resized |= adjust_size_by_resizeability( m_err , x , typename is_resizeable::type() ); + return resized; + } + + + template< class System , class StateInOut > + controlled_step_result try_step_v1( System system , StateInOut &x , time_type &t , time_type &dt ) + { + typename odeint::unwrap_reference< System >::type &sys = system; + m_dxdt_resizer.adjust_size( x , detail::bind( &controlled_error_bs_type::template resize_m_dxdt< StateInOut > , detail::ref( *this ) , detail::_1 ) ); + sys( x , m_dxdt.m_v ,t ); + return try_step( system , x , m_dxdt.m_v , t , dt ); + } + + + template< class StateInOut > + void extrapolate( size_t k , state_table_type &table , const value_matrix &coeff , StateInOut &xest ) + /* polynomial extrapolation, see http://www.nr.com/webnotes/nr3web21.pdf + uses the obtained intermediate results to extrapolate to dt->0 + */ + { + static const value_type val1 = static_cast< value_type >( 1.0 ); + for( int j=k-1 ; j>0 ; --j ) + { + m_algebra.for_each3( table[j-1].m_v , table[j].m_v , table[j-1].m_v , + typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][j] , -coeff[k][j] ) ); + } + m_algebra.for_each3( xest , table[0].m_v , xest , + typename operations_type::template scale_sum2< value_type , value_type >( val1 + coeff[k][0] , -coeff[k][0]) ); + } + + time_type calc_h_opt( time_type h , value_type error , size_t k ) const + /* calculates the optimal step size for a given error and stage number */ + { + BOOST_USING_STD_MIN(); + BOOST_USING_STD_MAX(); + using std::pow; + value_type expo( 1.0/(2*k+1) ); + value_type facmin = m_facmin_table[k]; + value_type fac; + if (error == 0.0) + fac=1.0/facmin; + else + { + fac = STEPFAC2 / pow BOOST_PREVENT_MACRO_SUBSTITUTION( error / STEPFAC1 , expo ); + fac = max BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(facmin/STEPFAC4) , min BOOST_PREVENT_MACRO_SUBSTITUTION( static_cast(1.0/facmin) , fac ) ); + } + return h*fac; + } + + controlled_step_result set_k_opt( size_t k , const inv_time_vector &work , const time_vector &h_opt , time_type &dt ) + /* calculates the optimal stage number */ + { + if( k == 1 ) + { + m_current_k_opt = 2; + return success; + } + if( (work[k-1] < KFAC1*work[k]) || (k == m_k_max) ) + { // order decrease + m_current_k_opt = k-1; + dt = h_opt[ m_current_k_opt ]; + return success; + } + else if( (work[k] < KFAC2*work[k-1]) || m_last_step_rejected || (k == m_k_max-1) ) + { // same order - also do this if last step got rejected + m_current_k_opt = k; + dt = h_opt[ m_current_k_opt ]; + return success; + } + else + { // order increase - only if last step was not rejected + m_current_k_opt = k+1; + dt = h_opt[ m_current_k_opt-1 ] * m_cost[ m_current_k_opt ] / m_cost[ m_current_k_opt-1 ] ; + return success; + } + } + + bool in_convergence_window( size_t k ) const + { + if( (k == m_current_k_opt-1) && !m_last_step_rejected ) + return true; // decrease stepsize only if last step was not rejected + return ( (k == m_current_k_opt) || (k == m_current_k_opt+1) ); + } + + bool should_reject( value_type error , size_t k ) const + { + if( k == m_current_k_opt-1 ) + { + const value_type d = m_interval_sequence[m_current_k_opt] * m_interval_sequence[m_current_k_opt+1] / + (m_interval_sequence[0]*m_interval_sequence[0]); + //step will fail, criterion 17.3.17 in NR + return ( error > d*d ); + } + else if( k == m_current_k_opt ) + { + const value_type d = m_interval_sequence[m_current_k_opt] / m_interval_sequence[0]; + return ( error > d*d ); + } else + return error > 1.0; + } + + default_error_checker< value_type, algebra_type , operations_type > m_error_checker; + modified_midpoint< state_type , value_type , deriv_type , time_type , algebra_type , operations_type , resizer_type > m_midpoint; + + bool m_last_step_rejected; + bool m_first; + + time_type m_dt_last; + time_type m_t_last; + time_type m_max_dt; + + size_t m_current_k_opt; + + algebra_type m_algebra; + + resizer_type m_dxdt_resizer; + resizer_type m_xnew_resizer; + resizer_type m_resizer; + + wrapped_state_type m_xnew; + wrapped_state_type m_err; + wrapped_deriv_type m_dxdt; + + int_vector m_interval_sequence; // stores the successive interval counts + value_matrix m_coeff; + int_vector m_cost; // costs for interval count + value_vector m_facmin_table; // for precomputed facmin to save pow calls + + state_table_type m_table; // sequence of states for extrapolation + + value_type STEPFAC1 , STEPFAC2 , STEPFAC3 , STEPFAC4 , KFAC1 , KFAC2; +}; + + +/******** DOXYGEN ********/ +/** + * \class bulirsch_stoer + * \brief The Bulirsch-Stoer algorithm. + * + * The Bulirsch-Stoer is a controlled stepper that adjusts both step size + * and order of the method. The algorithm uses the modified midpoint and + * a polynomial extrapolation compute the solution. + * + * \tparam State The state type. + * \tparam Value The value type. + * \tparam Deriv The type representing the time derivative of the state. + * \tparam Time The time representing the independent variable - the time. + * \tparam Algebra The algebra type. + * \tparam Operations The operations type. + * \tparam Resizer The resizer policy type. + */ + + /** + * \fn bulirsch_stoer::bulirsch_stoer( value_type eps_abs , value_type eps_rel , value_type factor_x , value_type factor_dxdt ) + * \brief Constructs the bulirsch_stoer class, including initialization of + * the error bounds. + * + * \param eps_abs Absolute tolerance level. + * \param eps_rel Relative tolerance level. + * \param factor_x Factor for the weight of the state. + * \param factor_dxdt Factor for the weight of the derivative. + */ + + /** + * \fn bulirsch_stoer::try_step( System system , StateInOut &x , time_type &t , time_type &dt ) + * \brief Tries to perform one step. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param x The state of the ODE which should be solved. Overwritten if + * the step is successful. + * \param t The value of the time. Updated if the step is successful. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + /** + * \fn bulirsch_stoer::try_step( System system , StateInOut &x , const DerivIn &dxdt , time_type &t , time_type &dt ) + * \brief Tries to perform one step. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param x The state of the ODE which should be solved. Overwritten if + * the step is successful. + * \param dxdt The derivative of state. + * \param t The value of the time. Updated if the step is successful. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + /** + * \fn bulirsch_stoer::try_step( System system , const StateIn &in , time_type &t , StateOut &out , time_type &dt ) + * \brief Tries to perform one step. + * + * \note This method is disabled if state_type=time_type to avoid ambiguity. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param in The state of the ODE which should be solved. + * \param t The value of the time. Updated if the step is successful. + * \param out Used to store the result of the step. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + + /** + * \fn bulirsch_stoer::try_step( System system , const StateIn &in , const DerivIn &dxdt , time_type &t , StateOut &out , time_type &dt ) + * \brief Tries to perform one step. + * + * This method tries to do one step with step size dt. If the error estimate + * is to large, the step is rejected and the method returns fail and the + * step size dt is reduced. If the error estimate is acceptably small, the + * step is performed, success is returned and dt might be increased to make + * the steps as large as possible. This method also updates t if a step is + * performed. Also, the internal order of the stepper is adjusted if required. + * + * \param system The system function to solve, hence the r.h.s. of the ODE. + * It must fulfill the Simple System concept. + * \param in The state of the ODE which should be solved. + * \param dxdt The derivative of state. + * \param t The value of the time. Updated if the step is successful. + * \param out Used to store the result of the step. + * \param dt The step size. Updated. + * \return success if the step was accepted, fail otherwise. + */ + + + /** + * \fn bulirsch_stoer::adjust_size( const StateIn &x ) + * \brief Adjust the size of all temporaries in the stepper manually. + * \param x A state from which the size of the temporaries to be resized is deduced. + */ + +} +} +} + +#endif // BOOST_NUMERIC_ODEINT_STEPPER_BULIRSCH_STOER_HPP_INCLUDED + +#endif // USE_BULRISCH_STOER_PATCH diff --git a/tests/testthat/test-DAISIE_plot_sims.R b/tests/testthat/test-DAISIE_plot_sims.R index c6da7367..93417224 100644 --- a/tests/testthat/test-DAISIE_plot_sims.R +++ b/tests/testthat/test-DAISIE_plot_sims.R @@ -1,5 +1,5 @@ test_that("Example 1", { - + skip("Plots: Run and manually inspect output") data(islands_1type_1000reps) expect_silent( DAISIE_plot_sims( @@ -9,6 +9,7 @@ test_that("Example 1", { }) test_that("Example 2", { + skip("Plots: Run and manually inspect output") data(islands_2types_1000reps) expect_silent( DAISIE_plot_sims( @@ -18,7 +19,7 @@ test_that("Example 2", { }) test_that("Plot plus one", { - + skip("Plots: Run and manually inspect output") data(islands_1type_1000reps) expect_silent( DAISIE_plot_sims( diff --git a/tests/testthat/test-DAISIE_plot_stt.R b/tests/testthat/test-DAISIE_plot_stt.R index b73b1ff0..45a9faea 100644 --- a/tests/testthat/test-DAISIE_plot_stt.R +++ b/tests/testthat/test-DAISIE_plot_stt.R @@ -1,4 +1,5 @@ test_that("use", { + skip("Plots: Run and manually inspect output") utils::data(islands_1type_1000reps, package = "DAISIE") plot_lists <- DAISIE:::DAISIE_convert_to_classic_plot( simulation_outputs = islands_1type_1000reps diff --git a/tests/testthat/test_odeint.R b/tests/testthat/test-odeint.R similarity index 85% rename from tests/testthat/test_odeint.R rename to tests/testthat/test-odeint.R index 6a47f4a9..eb03992a 100644 --- a/tests/testthat/test_odeint.R +++ b/tests/testthat/test-odeint.R @@ -33,18 +33,11 @@ test_that("odeint solvers give the same result as deSolve solvers", { loglik_rkd5 <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) methode <- 'odeint::bulirsch_stoer' loglik_bs <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) - methode <- 'odeint::rosenbrock4' - loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) - methode <- 'odeint::adams_bashforth_moulton_1' - DAISIE_CS_max_steps(100000000) - DAISIE_abm_factor(0.000001) - loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) expect_equal(loglik_lsodes, loglik_rkck54) expect_equal(loglik_lsodes, loglik_rkf78) expect_equal(loglik_lsodes, loglik_rkd5) expect_equal(loglik_lsodes, loglik_bs) - expect_equal(loglik_lsodes, loglik_rb) - expect_equal(loglik_lsodes, loglik_abm, tol = 1E-6) + pars1a <- pars1 pars1a[6] <- Inf diff --git a/tests/testthat/test-rosenbrock_abm1_steppers.R b/tests/testthat/test-rosenbrock_abm1_steppers.R new file mode 100644 index 00000000..6037871c --- /dev/null +++ b/tests/testthat/test-rosenbrock_abm1_steppers.R @@ -0,0 +1,28 @@ +test_that("rosenbrock4 and adams bashforth moulton1 work", { + skip("Too slow to run") + utils::data(Galapagos_datalist_2types) + pars1 <- c( + 0.195442017, + 0.087959583, + Inf, + 0.002247364, + 0.873605049, + 3755.202241, + 8.909285094, + 14.99999923, + 0.002247364, + 0.873605049, + 0.163 + ) + pars2 <- c(40, 11, 0, 1) + methode <- 'deSolve_R::lsodes' + loglik_lsodes_R <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::rosenbrock4' + loglik_rb <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + methode <- 'odeint::adams_bashforth_moulton_1' + DAISIE_CS_max_steps(100000000) + DAISIE_abm_factor(0.000001) + loglik_abm <- DAISIE_loglik_all(pars1, pars2, Galapagos_datalist_2types, methode = methode) + expect_equal(loglik_lsodes, loglik_rb) + expect_equal(loglik_lsodes, loglik_abm, tol = 1E-6) +})