Bug fix add_brt_table and proposal for new function to replace it

R/DAISIE_format_IW.R

@@ -175,7 +175,7 @@ add_brt_table <- function(island, full_table = FALSE) {
           for (k in 1:length(island[[i]]$all_colonisations)) {
             the_brts <- island[[i]]$all_colonisations[[k]]$event_times[-1]
             pos2 <- pos1 + length(the_brts)
-            ff <- matrix(ncol = 5, nrow = pos2 - pos1 + 1)
+            ff <- matrix(ncol = 5, nrow = pos2 - pos1)
             ff[1:(pos2 - pos1), 1] <- the_brts
             ff[1:(pos2 - pos1), 2] <- j
             ff[1:(pos2 - pos1), 3] <- seq(1, length(the_brts))
@@ -203,3 +203,89 @@ add_brt_table <- function(island, full_table = FALSE) {
   colnames(island[[1]]$brts_table) <- c("brt", "clade", "event", "endemic", "col")
   return(island)
 }
+
+#The function below makes more sense, but DAISIE_loglik_IW needs to be adjusted to it
+add_brt_table2 <- function(island, full_table = FALSE) {
+  island_age <- island[[1]]$island_age
+  island_top <- island[[1]]
+  if (length(island) == 1) {
+    brts_table <- matrix(ncol = 5, nrow = 1)
+    brts_table[1, ] <-  c(island_age, 0, 0, NA, NA)
+    island[[1]]$brts_table <- brts_table
+  } else {
+    island_top <- island[[1]]
+    island[[1]] <- NULL
+    btimes <- list()
+    for (i in 1:length(island)) {
+      btimes[[i]] <- island[[i]]$branching_times[-1]
+    }
+    island <- island[rev(order(sapply(btimes, "[", 1)))]
+    il <- unlist(island)
+    stac1s <- which(il[which(names(il) == "stac")] == "1")
+    stac5s <- which(il[which(names(il) == "stac")] == "5")
+    stac1_5s <- sort(c(stac1s, stac5s))
+    if (length(stac1_5s) != 0) {
+      if (length(stac1_5s) == length(island)) {
+        brts_table <- matrix(ncol = 5, nrow = 1)
+        brts_table[1, ] <- c(island_age, 0, 0, NA, NA)
+        island_no_stac1or5 <- NULL
+      } else {
+        island_no_stac1or5 <- island[-stac1_5s]
+      }
+    }
+    if (length(stac1_5s) == 0) {
+      island_no_stac1or5 <- island
+    }
+    if (length(island_no_stac1or5) != 0) {
+      btimes <- list()
+      for (i in 1:length(island_no_stac1or5)) {
+        btimes[[i]] <- island_no_stac1or5[[i]]$branching_times[-1]
+      }
+      brts <- rev(sort(unlist(btimes)))
+      brts_IWK <- NULL
+      pos1 <- 0
+      for (i in 1:length(btimes)) {
+        the_stac <- island_no_stac1or5[[i]]$stac
+        if(!is.null(island[[i]]$all_colonisations) & full_table == TRUE) {
+          if(length(island[[i]]$all_colonisations) > 0) print(i)
+          for (k in 1:length(island[[i]]$all_colonisations)) {
+            the_brts <- island[[i]]$all_colonisations[[k]]$event_times[-1]
+            pos2 <- pos1 + length(the_brts)
+            ff <- matrix(ncol = 5, nrow = pos2 - pos1)
+            ff[1:(pos2 - pos1), 1] <- the_brts
+            ff[1:(pos2 - pos1), 2] <- i
+            ff[1:(pos2 - pos1), 3] <- seq(1, length(the_brts))
+            ff[1:(pos2 - pos1), 4] <- (the_stac == 2) + (the_stac == 3)
+            ff[1:(pos2 - pos1), 5] <- k
+            brts_IWK <- rbind(brts_IWK,ff)
+            pos1 <- pos2
+          }
+        } else {
+          the_brts <- btimes[[i]]
+          pos2 <- pos1 + length(the_brts)
+          ff <- matrix(ncol = 5, nrow = pos2 - pos1)
+          ff[1:(pos2 - pos1), 1] <- the_brts
+          ff[1:(pos2 - pos1), 2] <- i
+          ff[1:(pos2 - pos1), 3] <- seq(1, length(the_brts))
+          ff[1:(pos2 - pos1), 4] <- (the_stac == 2) + (the_stac == 3)
+          ff[1:(pos2 - pos1), 5] <- 1
+          brts_IWK <- rbind(brts_IWK,ff)
+          pos1 <- pos2
+        }
+      }
+      brts_table <- brts_IWK[rev(order(brts_IWK[, 1])), ]
+      brts_table <- rbind(c(island_age, 0, 0, NA, NA), brts_table)
+    }
+    island_top$brts_table <- brts_table
+    if (length(stac1_5s) != 0) {
+      for (i in 1:length(stac1_5s)) {
+        island[[length(island) + 1]] <- island[[stac1_5s[i]]]
+        island[[stac1_5s[i]]] <- NULL
+        stac1_5s <- stac1_5s - 1
+      }
+    }
+    island <- append(list(island_top), island)
+  }
+  colnames(island[[1]]$brts_table) <- c("brt", "clade", "event", "endemic", "col")
+  return(island)
+}