R/DAISIE_loglik_IW0.R

In xieshu95/DAISIE_new: Dynamical Assembly of Islands by Speciation, Immigration and Extinction

kmini0 <- function(dec2binmatk, lxm, lxe, sysdim = dim(dec2binmatk)[1]) {
   posc <- Matrix::rowSums(dec2binmatk)
   negc <- log2(sysdim) - posc
   kmin <- rep(negc, each = lxm * lxe)
   dim(kmin) <- c(lxm, lxe, sysdim)
   ki <- kimat(dec2binmatk)
   res <- list(kmin = kmin, ki = ki)
   return(res)
}

nndivdep0 <- function(lxm, lxe, sysdim, Kprime, M, k) {
  nnm <- c(0, 0:(lxm + 1))
  nne <- c(0, 0, 0:(lxe + 1))
  lnnm <- length(nnm)
  lnne <- length(nne)
  nn <- matrix(nnm, nrow = lnnm, ncol = lnne, byrow = F) +
       matrix(nne, nrow = lnnm, ncol = lnne, byrow = T)
  nn <- replicate(sysdim, nn)
  nil2lxm <- 2:(lxm + 1)
  nil2lxe <- 3:(lxe + 2)
  nile <- rep(1, lxe)
  allc <- 1:sysdim
  divdepfac <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)),
                    1 - (nn + k) / Kprime)
  divdepfacmin1 <- pmax(array(0, dim = c(lxm + 3, lxe + 4, sysdim)),
                        1 - (nn + k - 1) / Kprime)
  divdepfac <- divdepfac[nil2lxm, nil2lxe, allc]
  divdepfacmin1 <- divdepfacmin1[nil2lxm, nil2lxe, allc]
  Mminm <- M - nn[nil2lxm, nile, allc]
  res <- list(nn = nn,
              divdepfac = divdepfac,
              divdepfacmin1 = divdepfacmin1,
              Mminm = Mminm)
  return(res)
}

DAISIE_loglik_rhs_IW0 <- function(t, x, pars) {
  lac <- pars[[1]][1]
  mu <- pars[[1]][2]
  Kprime <- pars[[1]][3]
  gam <- pars[[1]][4]
  laa <- pars[[1]][5]
  M <- pars[[1]][6]
  kk <- pars[[2]]
  ddep <- pars[[3]]
  lxm <- pars[[4]]$lxm
  lxe <- pars[[4]]$lxe
  sysdim <- pars[[4]]$sysdim
  kmin <- pars[[5]]$kmin
  kplus <- kk - kmin
  ki <- pars[[5]]$ki
  nn <- pars[[6]]$nn
  divdepfac <- pars[[6]]$divdepfac
  divdepfacmin1 <- pars[[6]]$divdepfacmin1
  Mminm <- pars[[6]]$Mminm

  dim(x) <- c(lxm, lxe, sysdim)
  xx <- array(0, dim = c(lxm + 2, lxe + 3, sysdim))
  xx[2:(lxm + 1), 3:(lxe + 2), 1:sysdim] <- x
  nil2lxm <- 2:(lxm + 1)
  nil2lxe <- 3:(lxe + 2)
  allc <- 1:sysdim
  nile <- rep(1, lxe)
  nilm <- rep(1, lxm)
  Mminm2 <- Mminm + 1 - kmin
  Mminm2[Mminm2 < 0] <- 0
  Mminm[Mminm < 0] <- 0
  if (sysdim == 1) {
    dim(Mminm) <- c(lxm, lxe)
  }
  dx <-  gam * divdepfacmin1 * Mminm2 * xx[nil2lxm - 1, nil2lxe, allc] + #immigration
        mu * nn[nil2lxm + 1, nile, allc] * xx[nil2lxm + 1, nil2lxe, allc] + #extinction non-endemics
        mu * nn[nilm, nil2lxe + 1, allc] * xx[nil2lxm, nil2lxe + 1, allc] + #extinction endemics
        lac * divdepfacmin1 * nn[nil2lxm + 1, nile, allc] * xx[nil2lxm + 1, nil2lxe - 2, allc] + #cladogenesis non-endemics
        lac * divdepfacmin1 * nn[nilm, nil2lxe - 1, allc] * xx[nil2lxm, nil2lxe - 1, allc] + #cladogenesis endemics
        2 * kplus * lac * divdepfacmin1 * xx[nil2lxm, nil2lxe - 1, allc] + #cladogenesis species in tree
        laa * nn[nil2lxm + 1, nile, allc] * xx[nil2lxm + 1, nil2lxe - 1, allc] + #anagenesis non-endemics
       -(laa * (nn[nil2lxm, nile, allc] + kmin) + (gam * divdepfac * Mminm) +
        (lac * divdepfac + mu) * (nn[nil2lxm, nil2lxe, allc] + kk)) * xx[nil2lxm, nil2lxe, allc]
  if (sysdim > 1) {
      dx <- dx +
      laa * tensor::tensor(xx[nil2lxm, nil2lxe, allc], ki, 3, 2) + # anagenesis in colonizing lineage
      2 * lac * divdepfacmin1 * tensor::tensor(xx[nil2lxm, nil2lxe - 1, allc], ki, 3, 2) # cladogenesis in colonizing lineage
  }
  dim(dx) <- c(sysdim * lxm * lxe, 1)
  return(list(dx))
}


DAISIE_loglik_IW0 <- function(
  pars1,
  pars2,
  datalist,
  methode = "ode45",
  abstolint = 1E-16,
  reltolint = 1E-14,
  verbose = verbose
  ) {
  # pars1 = model parameters
  # - pars1[1] = lac = (initial) cladogenesis rate
  # - pars1[2] = mu = extinction rate
  # - pars1[3] = K = maximum number of species possible in the clade
  # - pars1[4] = gam = (initial) immigration rate
  # - pars1[5] = laa = (initial) anagenesis rate
  # - pars1[6] = M = number of mainland species
  # pars2 = model settings
  # - pars2[1] = lx = length of ODE variable x
  # - pars2[2] = ddep = diversity-dependent model,mode of diversity-dependence
  #  . ddep == 0 : no diversity-dependence
  #  . ddep == 1 : linear dependence in speciation rate (anagenesis and cladogenesis)
  #  . ddep == 11 : linear dependence in speciation rate and immigration rate
  #  . ddep == 3 : linear dependence in extinction rate
  # - pars2[3] = cond = conditioning
  #  . cond == 0 : no conditioning
  #  . cond == 1 : conditioning on presence on the island (not used in this single loglikelihood)
  # - pars2[4] = parameters and likelihood should be printed (1) or not (0)
  
  brts <- c(-abs(datalist[[1]]$brts_table[, 1]), 0)
  clade <- datalist[[1]]$brts_table[, 2]
  event <- datalist[[1]]$brts_table[, 3]
  pars1 <- as.numeric(pars1)

  ddep <- pars2[2]
  cond <- pars2[3]
  
  lac <- pars1[1]
  mu <- pars1[2]
  Kprime <- pars1[3]
  if (ddep == 0)
  {
    Kprime <- Inf
  }
  gam <- pars1[4]
  laa <- pars1[5]
  M <- pars1[6]
  
  if (min(pars1) < 0) {
    cat("One or more parameters are negative.\n")
    loglik <- -Inf
    return(loglik)
  }
  if ((ddep == 1 | ddep == 11) & ceiling(Kprime) < length(brts)) {
    if (verbose) {
      cat("The proposed value of K is incompatible with the number of species
          in the clade. Likelihood for this parameter set
          will be set to -Inf. \n")
    }
    loglik <- -Inf
    return(loglik)
  }
  if (ddep == 1 | ddep == 11) {
    lx <- min(1 + ceiling(Kprime), DDD::roundn(pars2[1]) )
  } else {
    lx <- DDD::roundn(pars2[1])
  }
  lxm <- min(lx, M + 1)
  lxe <- lx
  sysdimchange <- 1
  sysdim <- 1
  totdim <- lxm * lxe * sysdim
  probs <- rep(0, totdim)
  probs[1] <- 1
  loglik <- 0
  expandvec <- NULL
  for (k in 0:(length(brts) - 2)) {
    if (pars2[4] == 2) {
       print(paste("k = ", k, ", sysdim = ", sysdim, sep = ""))
       utils::flush.console()
    }
    dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim)
    if (sysdimchange == 1) {
      if (sysdim > 1) {
        dec2binmatk <- dec2binmat(log2(sysdim))
        kmi <- kmini0(dec2binmatk, lxm, lxe, sysdim)
      } else if (sysdim == 1) {
        kmi <- list(kmin = 0, ki = NULL)
      }
      sysdimchange <- 0
    }
    nndd <- nndivdep0(lxm, lxe, sysdim, Kprime, M, k)
    parslist <- list(pars = pars1, kk = k, ddep = ddep, dime = dime, kmi = kmi, nndd = nndd)
    y <- deSolve::ode(y = probs, times = brts[(k + 1):(k + 2)], func = DAISIE_loglik_rhs_IW, parms = parslist, rtol = reltolint, atol = abstolint, method = methode)
    probs <- y[2, 2:(totdim + 1)]
    cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]]
    dim(probs) <- c(lxm, lxe, sysdim)
    
    if (k < (length(brts) - 2)) {
      divdepfac <- nndd$divdepfac
      if (event[k + 2] == 1) {
        #Mminmminl = nndd$Mminm - kmi$kmin
        Mminmminl <- nndd$Mminm - clade[k + 1]
        Mminmminl[Mminmminl < 0] <- 0
        probs <- gam * divdepfac * Mminmminl * probs[ , , 1:sysdim]
        probs <- c(probs, rep(0, totdim))

        sysdim <- sysdim * 2
        expandvec <- c(expandvec, clade[k + 2])
        sysdimchange <- 1
      } else {
        probs <- lac * divdepfac * probs[ , , 1:sysdim]
        if (event[k + 2] == 2) {
          tocollapse <- which(expandvec == clade[k + 2])
          sr <- selectrows(sysdim, 2 ^ (tocollapse - 1))
          probs <- probs[ , , sr[, 1]] + probs[ , , sr[, 2]]
          sysdim <- sysdim / 2
          dim(probs) <- c(lxm, lxe, sysdim)
          expandvec <- expandvec[-tocollapse]
          sysdimchange <- 1
        }
      }
      cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]]
      totdim <- lxm * lxe * sysdim
      dim(probs) <- c(totdim, 1)
      #print(head(probs,n = 5))
    }
  }
  dim(probs) <- c(lxm, lxe, sysdim)
  expandedclades <- which(pracma::histc(clade, 1:length(clade))$cnt == 1)
  status <- rep(0, lexpandedclades <- length(expandedclades))
  if (lexpandedclades > 0) {
    for (i in lexpandedclades:1) {
      if (datalist[[1 + expandedclades[i]]]$stac == 2) {
        status[i] <- 1
      }
    }
  }
  endemic <- 0
  nonendemic <- 0
  for (i in 2:length(datalist)) {
      endemic <- endemic + (datalist[[i]]$stac == 5)
      nonendemic <- nonendemic + (datalist[[i]]$stac == 1) + (datalist[[i]]$stac == 3)
  }
  if (length(status) > 0) {
    decstatus <- bin2dec(status)
  } else {
    decstatus <- 0
  }
  print(loglik + log(probs))
  loglik <- loglik + log(probs[1 + nonendemic, 1 + endemic, 1 + decstatus])

  if (cond > 0) {
    sysdim <- 1
    totdim <- lxm * lxe * sysdim
    dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim)
    probs <- rep(0, totdim)
    probs[1] <- 1
    kmi <- list(kmin = 0, ki = NULL)
    nndd <- nndivdep0(lxm, lxe, sysdim, Kprime, M, k = 0)
    parslist <- list(pars = pars1, kk = k, ddep = ddep, dime = dime, kmi = kmi, nndd = nndd)
    y <- deSolve::ode(y = probs, times = brts[(k + 1):(k + 2)], func = DAISIE_loglik_rhs_IW0, parms = parslist, rtol = reltolint, atol = abstolint, method = methode)
    probs <- y[2, 2:(totdim + 1)]
    dim(probs) <- c(lxm, lxe, sysdim)
    logcond <- log(1 - probs[1, 1, 1])
    loglik <- loglik - logcond
  }
  
  if (pars2[4] > 0) {
    s1 <- sprintf("Parameters: %f %f %f %f %f %d", pars1[1], pars1[2], pars1[3], pars1[4], pars1[5], pars1[6])
    s2 <- sprintf(", Loglikelihood: %f", loglik)
    cat(s1, s2, "\n", sep = "")
    utils::flush.console()
  }
  return(as.numeric(loglik))
}

DAISIE_loglik_IW_M1 <- function(
  pars1,
  pars2,
  brts,
  stac,
  missnumspec,
  methode = "ode45",
  abstolint = 1E-16,
  reltolint = 1E-14,
  verbose
  ) {
  # pars1 = model parameters
  # - pars1[1] = lac = (initial) cladogenesis rate
  # - pars1[2] = mu = extinction rate
  # - pars1[3] = K = maximum number of species possible in the clade
  # - pars1[4] = gam = (initial) immigration rate
  # - pars1[5] = laa = (initial) anagenesis rate
  # pars2 = model settings
  # - pars2[1] = lx = length of ODE variable x
  # - pars2[2] = ddep = diversity-dependent model,mode of diversity-dependence
  #  . ddep == 0 : no diversity-dependence
  #  . ddep == 1 : linear dependence in speciation rate (anagenesis and cladogenesis)
  #  . ddep == 11 : linear dependence in speciation rate and immigration rate
  #  . ddep == 3 : linear dependence in extinction rate
  # - pars2[3] = cond = conditioning
  #  . cond == 0 : no conditioning
  #  . cond == 1 : conditioning on presence on the island (not used in this single loglikelihood)
  # - pars2[4] = parameters and likelihood should be printed (1) or not (0)
  
  if (is.na(pars2[4])) {
    pars2[4] <- 0
  }
  ddep <- pars2[2]
  cond <- pars2[3]
  brts <- c(-abs(brts), 0)
  pars1 <- as.numeric(pars1)
  lac <- pars1[1]
  mu <- pars1[2]
  Kprime <- pars1[3]
  if (ddep == 0) {
    Kprime <- Inf
  }
  gam <- pars1[4]
  laa <- pars1[5]
  pars1[6] <- 1
  M <- pars1[6]
  
  if (min(pars1) < 0) {
    cat("One or more parameters are negative.\n")
    loglik <- -Inf
    return(loglik)
  }
  if (ddep == 1 | ddep == 11) {
    lx <- min(1 + ceiling(Kprime), round(pars2[1]) )
  } else {
    lx <- DDD::roundn(pars2[1])
  }
  lxm <- min(lx, M + 1)
  lxe <- lx
  sysdimchange <- 1
  sysdim <- 1
  totdim <- lxm * lxe * sysdim
  probs <- rep(0, totdim)
  probs[1] <- 1
  loglik <- 0
  expandvec <- NULL
  for (k in 0:(length(brts) - 2)) {
    if (pars2[4] == 2) {
      print(paste("k = ", k, ", sysdim = ", sysdim, sep = ""))
      utils::flush.console()
    }
    dime <- list(lxm = lxm, lxe = lxe, sysdim = sysdim)
    if (sysdimchange == 1) {
      if (sysdim > 1) {
        dec2binmatk <- dec2binmat(log2(sysdim))
        kmi <- kmini0(dec2binmatk, lxm, lxe, sysdim)
      } else if (sysdim == 1) {
        kmi <- list(kmin = 0, ki = NULL)
      }
      sysdimchange <- 0
    }
    nndd <- nndivdep0(lxm, lxe, sysdim, Kprime, M, k)
    parslist <- list(pars = pars1, kk = k, ddep = ddep, dime = dime, kmi = kmi, nndd = nndd)
    y <- deSolve::ode(y = probs, times = brts[(k + 1):(k + 2)], func = DAISIE_loglik_rhs_IW0, parms = parslist, rtol = reltolint, atol = abstolint, method = methode)
    probs <- y[2, 2:(totdim + 1)]
    cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]]      
    dim(probs) <- c(lxm, lxe, sysdim)
    
    if (k < (length(brts) - 2)) {
      divdepfac <- nndd$divdepfac
      if(k == 0) {
        #Mminmminl = nndd$Mminm - kmi$kmin
        Mminmminl <- nndd$Mminm
        Mminmminl[Mminmminl < 0] <- 0
        probs <- gam * divdepfac * Mminmminl * probs[ , , 1:sysdim]
        probs <- c(probs, rep(0, totdim))
        sysdim <- 2
        sysdimchange <- 1
      } else {
        probs <- lac * divdepfac * probs[ , , 1:sysdim]
        if (k == 1) {
          probs <- probs[ , , 1] + probs[ , , 2]
          sysdim <- 1
          dim(probs) <- c(lxm, lxe, sysdim)
          sysdimchange <- 1
        }
      }
      cp <- checkprobs2(NA, loglik, probs, verbose); loglik = cp[[1]]; probs = cp[[2]]
      totdim <- lxm * lxe * sysdim
      dim(probs) <- c(totdim, 1)
    }
  }
  dim(probs) <- c(lxm, lxe, sysdim)
  endemic <- (stac == 5)
  nonendemic <- (stac == 1) + (stac == 3)
  decstatus <- (stac == 2) * (sysdim > 1) #when stac = 4, decstatus = 0
  #print(probs)
  loglik <- loglik + log(probs[1 + nonendemic, 1 + endemic, 1 + decstatus])
  
  if (pars2[4] > 0) {
    s1 <- sprintf("Parameters: %f %f %f %f %f %d", pars1[1], pars1[2], pars1[3], pars1[4], pars1[5], pars1[6])
    s2 <- sprintf(", Loglikelihood: %f", loglik)
    cat(s1, s2, "\n", sep = "")
    utils::flush.console()
  }
  return(as.numeric(loglik))
}

DAISIE_loglik_IW0_choosepar <- function(
  trparsopt,
  trparsfix,
  idparsopt,
  idparsfix,
  pars2,
  datalist,
  methode,
  abstolint,
  reltolint
) {
  trpars1 <- rep(0, 6)
  trpars1[idparsopt] <- trparsopt
  if (length(idparsfix) != 0) {
    trpars1[idparsfix] <- trparsfix
  }
  if (max(trpars1) > 1 | min(trpars1) < 0) {
    loglik <- -Inf
  } else {
    pars1 <- trpars1 / (1 - trpars1)
    if (min(pars1) < 0) {
      loglik <- -Inf
    } else {
      loglik <- DAISIE_loglik_IW0(pars1, pars2, datalist, methode, abstolint, reltolint)
    }
    if (is.nan(loglik) || is.na(loglik)) {
      cat("There are parameter values used which cause numerical problems.\n")
      loglik <- -Inf
    }
  }
  return(loglik)
}