Skip to content

Commit

Permalink
Drop all fortran code
Browse files Browse the repository at this point in the history
  • Loading branch information
@richfitz
richfitz committed Dec 20, 2023
1 parent e9c1346 commit c7242dd
Show file tree
Hide file tree
Showing 25 changed files with 16 additions and 3,561 deletions.
68 changes: 0 additions & 68 deletions R/expm.R
View file
Original file line number Diff line number Diff line change
@@ -1,68 +0,0 @@
## This is significantly less efficient than proper checkpointing
## within ddexpmv -> ddexpmvi. However, the interface should be OK.
expmv.expokit.dense <- function(Q, t, v) {
n <- as.integer(length(v))
nt <- length(t)
if ( nt != 1 ) {
ret <- matrix(0, length(v), nt)
for ( i in seq_len(nt) ) {
res <- .Fortran(f_ddexpmv, Q, n, v, t[i],
out = numeric(n), iflag = numeric(1))
if ( res$iflag != 0 )
stop("Expokit failed with code ", res$iflag)
ret[,i] <- res$out
}
} else {
res <- .Fortran(f_ddexpmv, Q, n, v, t, out = numeric(n), iflag = numeric(1))
if ( res$iflag != 0 )
stop("Expokit failed with code ", res$iflag)
ret <- res$out
}
ret
}

expmv.expokit.sparse <- function(Q.sparse, t, v, tol=1e-8) {
if ( is.matrix(Q.sparse) )
Q.sparse <- expm.expokit.sparse.pars(Q.sparse)
Q <- Q.sparse$Q
iq <- Q.sparse$iq
jq <- Q.sparse$jq
qnorm <- Q.sparse$qnorm

n <- as.integer(length(v))
lt <- as.integer(length(t))

res <- .Fortran(f_dsexpmvi, Q, n, iq, jq, length(iq),
qnorm, v, t, lt, tol, out = numeric(n*lt),
iflag = integer(1))
if ( res$iflag != 0 )
stop("Expokit failed with code ", res$iflag)

matrix(res$out, n, lt)
}

## Convert a sparse Q matrix into COO format in preparation for
## running expmv(sparse) on it.
expm.expokit.sparse.pars <- function(Q) {
idx <- which(Q != 0, TRUE)
nz <- as.integer(nrow(idx))
qq <- Q[idx]
list(Q=qq, iq=idx[,1], jq=idx[,2], nz=nrow(idx), qnorm=max(abs(qq)))
}

expm.dense <- function(Q, t) {
if ( !is.matrix(Q) )
stop("Q must be a matrix")
n <- as.integer(nrow(Q))
if ( n != ncol(Q) )
stop("Q must be a square matrix")
if ( n < 1 )
stop("Q must have positive dimensions")
t <- as.numeric(t)
if ( length(t) != 1 )
stop("t must be a scalar")
if ( !is.finite(t) )
stop("t must be finite")

matrix(.Fortran(f_dexpmf, Q, n, t, out = numeric(n*n), numeric(1))$out, n, n)
}
63 changes: 2 additions & 61 deletions R/model-bisse-unresolved.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,32 +15,7 @@
## faster, sometimes slower. Play around.
bucexp <- function(nt, la0, la1, mu0, mu1, q01, q10, t, scal=1,
tol=1e-10, m=15) {
stopifnot(all(c(la0, la1, mu0, mu1, q01, q10) >= 0))
lt <- length(t)
stopifnot(lt > 0)
stopifnot(scal >= 1)
stopifnot(m > 1)
n <- (nt*(nt+1)/2+1)
res <- .Fortran(f_bucexp,
nt = as.integer(nt),
la0 = as.numeric(la0),
la1 = as.numeric(la1),
mu0 = as.numeric(mu0),
mu1 = as.numeric(mu1),
q01 = as.numeric(q01),
q10 = as.numeric(q10),
t = as.numeric(t),
lt = as.integer(lt),
scal = as.numeric(scal),
tol = as.numeric(tol),
m = as.integer(m),
w = numeric(2*n*lt),
flag = integer(1))
if ( res$flag < 0 )
stop("Failure in the Fortran code")
else if ( res$flag > 0 )
cat(sprintf("Flag: %d\n", res$flag))
array(res$w, c(n, lt, 2))
stop("Unresolved clade calculations no longer supported (since 0.10.0)")
}

## This is more useful; given model parameters, and vectors of times
Expand All @@ -52,41 +27,7 @@ bucexp <- function(nt, la0, la1, mu0, mu1, q01, q10, t, scal=1,
## probabilities of extinction.
bucexpl <- function(nt, la0, la1, mu0, mu1, q01, q10, t, Nc, nsc, k,
scal=1, tol=1e-10, m=15) {
tt <- sort(unique(t))
ti <- as.integer(factor(t))

lt <- length(tt)
lc <- length(Nc)
stopifnot(lc > 0, lt > 0, max(Nc) < nt,
length(Nc) == lc, length(nsc) == lc, length(k) == lc,
all(nsc <= Nc), all(nsc >= 0), all(k <= nsc),
scal >= 1, m > 1)

res <- .Fortran(f_bucexpl,
nt = as.integer(nt),
la0 = as.numeric(la0),
la1 = as.numeric(la1),
mu0 = as.numeric(mu0),
mu1 = as.numeric(mu1),
q01 = as.numeric(q01),
q10 = as.numeric(q10),
t = as.numeric(tt),
lt = as.integer(lt),
ti = as.integer(ti),
Nc = as.integer(Nc),
nsc = as.integer(nsc),
k = as.integer(k),
lc = as.integer(lc),
scal = as.numeric(scal),
tol = as.numeric(tol),
m = as.integer(m),
ans = numeric(4*lc),
flag = integer(1))
if ( res$flag < 0 )
stop("Failure in the Fortran code")
else if ( res$flag > 0 )
cat(sprintf("Flag: %d\n", res$flag))
matrix(res$ans, ncol=4)
stop("Unresolved clade calculations no longer supported (since 0.10.0)")
}

## Compute the PDF of the hypergeometric distribution, using the same
Expand Down
71 changes: 2 additions & 69 deletions R/model-bisseness-unresolved.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,36 +15,7 @@
## faster, sometimes slower. Play around.
nucexp <- function(nt, la0, la1, mu0, mu1, q01, q10, p0c, p0a,
p1c, p1a, t, scal=1, tol=1e-10, m=15) {
stopifnot(all(c(la0, la1, mu0, mu1, q01, q10, p0c, p0a, p1c, p1a) >= 0))
lt <- length(t)
stopifnot(lt > 0)
stopifnot(scal >= 1)
stopifnot(m > 1)
n <- (nt*(nt+1)/2+1)
res <- .Fortran(f_nucexp,
nt = as.integer(nt),
la0 = as.numeric(la0),
la1 = as.numeric(la1),
mu0 = as.numeric(mu0),
mu1 = as.numeric(mu1),
q01 = as.numeric(q01),
q10 = as.numeric(q10),
p0c = as.numeric(p0c),
p0a = as.numeric(p0a),
p1c = as.numeric(p1c),
p1a = as.numeric(p1a),
t = as.numeric(t),
lt = as.integer(lt),
scal = as.numeric(scal),
tol = as.numeric(tol),
m = as.integer(m),
w = numeric(2*n*lt),
flag = integer(1))
if ( res$flag < 0 )
stop("Failure in the Fortran code")
else if ( res$flag > 0 )
cat(sprintf("Flag: %d\n", res$flag))
array(res$w, c(n, lt, 2))
stop("Unresolved clade calculations no longer supported (since 0.10.0)")
}

## This is more useful; given model parameters, and vectors of times
Expand All @@ -56,45 +27,7 @@ nucexp <- function(nt, la0, la1, mu0, mu1, q01, q10, p0c, p0a,
## probabilities of extinction.
nucexpl <- function(nt, la0, la1, mu0, mu1, q01, q10, p0c, p0a,
p1c, p1a, t, Nc, nsc, k, scal=1, tol=1e-10, m=15) {
tt <- sort(unique(t))
ti <- as.integer(factor(t))

lt <- length(tt)
lc <- length(Nc)
stopifnot(lc > 0, lt > 0, max(Nc) < nt,
length(Nc) == lc, length(nsc) == lc, length(k) == lc,
all(nsc <= Nc), all(nsc >= 0), all(k <= nsc),
scal >= 1, m > 1)

res <- .Fortran(f_nucexpl,
nt = as.integer(nt),
la0 = as.numeric(la0),
la1 = as.numeric(la1),
mu0 = as.numeric(mu0),
mu1 = as.numeric(mu1),
q01 = as.numeric(q01),
q10 = as.numeric(q10),
p0c = as.numeric(p0c),
p0a = as.numeric(p0a),
p1c = as.numeric(p1c),
p1a = as.numeric(p1a),
t = as.numeric(tt),
lt = as.integer(lt),
ti = as.integer(ti),
Nc = as.integer(Nc),
nsc = as.integer(nsc),
k = as.integer(k),
lc = as.integer(lc),
scal = as.numeric(scal),
tol = as.numeric(tol),
m = as.integer(m),
ans = numeric(4*lc),
flag = integer(1))
if ( res$flag < 0 )
stop("Failure in the Fortran code")
else if ( res$flag > 0 )
cat(sprintf("Flag: %d\n", res$flag))
matrix(res$ans, ncol=4)
stop("Unresolved clade calculations no longer supported (since 0.10.0)")
}

## Construct the transition matrix. Again, non-R style, as this was
Expand Down
31 changes: 0 additions & 31 deletions R/model-mkn-expokit.R
View file

This file was deleted.

2 changes: 1 addition & 1 deletion R/model-mkn.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ rootfunc.mkn <- function(res, pars, root, root.p, intermediates) {
make.all_branches.mkn <- function(cache, control) {
if ( control$method == "ode" ) {
if ( !is.null(control$backend) && control$backend == "expokit" )
make.all_branches.mkn.expokit(cache, control)
stop("expokit no longer suppported")
else
make.all_branches.dtlik(cache, control, initial.conditions.mkn)
} else { # method == "pij"
Expand Down
4 changes: 4 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,10 @@ If fftw is not found, installation will continue, but the (relatively)
fast C based QuaSSE integration will not be available. The R based
fft integrator and the method-of-lines integrator will be available.

## Unresolved clades

As of version 0.10.0, diversitree can no longer work with unresolved clades (FitzJohn, Maddison and Otto 2009's method), due to the package being long in retirement from development and difficulties adapting the Fortran code to meet CRAN's requirements. Users can install an older package from github (e.g., with `remotes::install_github("mrc-ide/diversitree@e587755")` to install the last released version that contained this code). This will require a working Fortran compiler. Alternatively, a suitably motivated person could restore the code (reverting changes contained in `aaaa`) and patch the code to work with the current version of `flang` as used by CRAN.

## Branches

The "master" branch contains the bleeding edge version of diversitree.
Expand Down
48 changes: 0 additions & 48 deletions inst/tests/test-expm.R
View file

This file was deleted.

3 changes: 3 additions & 0 deletions man/make.bisse.Rd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,9 @@ starting.point.bisse(tree, q.div=5, yule=FALSE)
}
\section{Unresolved clade information}{
Since 0.10.10 this is no longer supported. See the package README for
more information.
This must be a \code{\link{data.frame}} with at least the four columns
\itemize{
\item \code{tip.label}, giving the name of the tip to which the data
Expand Down
17 changes: 0 additions & 17 deletions man/make.bisse.split.Rd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,23 +115,6 @@ lik.s(pars4) - lik.b(pars)
## You can use the labelled nodes rather than indices:
lik.s2 <- make.bisse.split(phy, phy$tip.state, nodes)
identical(lik.s(pars4), lik.s2(pars4))

## This also works where some tips are unresolved clades. Here are a
## few:
unresolved <-
data.frame(tip.label=c("sp12", "sp32", "sp9", "sp22", "sp11"),
Nc=c(2,5,3,2,5), n0=c(1, 4, 3, 2, 4), n1=c(1, 1, 0, 0, 1))

## Plain BiSSE with unresolved clades:
lik.u.b <- make.bisse(phy, phy$tip.state, unresolved=unresolved)
lik.u.b(pars) # -139.3688

## Split BiSSE with unresolved clades:
lik.u.s <- make.bisse.split(phy, phy$tip.state, nodes,
unresolved=unresolved)
lik.u.s(pars4) # -139.3688

lik.u.b(pars) - lik.u.s(pars4) # numerical error only
}

\author{Richard G. FitzJohn}
Expand Down
Loading

0 comments on commit c7242dd

Please sign in to comment.