Skip to content

Commit

Permalink
Merge pull request #155 from rsetienne/develop
Browse files Browse the repository at this point in the history 
v4.3.3
  • Loading branch information
@rsetienne
rsetienne authored Mar 10, 2023
2 parents 2d1c832 + b373327 commit a19c458
Show file tree
Hide file tree
Showing 13 changed files with 860 additions and 144 deletions.
3 changes: 3 additions & 0 deletions .Rbuildignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,3 +16,6 @@
^\.covrignore$
^LICENSE$
^\.zenodo\.json$

^\.vscode$
^source_odeint\.R
5 changes: 2 additions & 3 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -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:
Expand Down Expand Up @@ -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,
<DOI:10.1111/ele.12461>.
NeedsCompilation: yes
SystemRequirements: C++14
Encoding: UTF-8
VignetteBuilder: knitr
URL: https://github.com/rsetienne/DAISIE
Expand Down
19 changes: 19 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand Down
248 changes: 124 additions & 124 deletions src/DAISIE_loglik_rhs_FORTRAN.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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] +
Expand All @@ -172,23 +172,23 @@ 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)) + &
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 + 1) = 0
dConc(2 * N + 1) = 0

ENDDO

Expand All @@ -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)
Expand Down Expand Up @@ -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] +
Expand All @@ -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] +
Expand All @@ -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] +
Expand All @@ -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)) * &
Expand All @@ -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)
Expand Down Expand Up @@ -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] +
Expand All @@ -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
Expand Down
Loading

0 comments on commit a19c458

Please sign in to comment.